Multi-omics analysis identifies essential regulators of mitochondrial stress response in two wild-type C. elegans strains.

The mitochondrial unfolded protein response (UPRmt) is a promising pharmacological target for aging and age-related diseases. However, the integrative analysis of the impact of UPRmt activation on different signaling layers in animals with different genetic backgrounds is lacking. Here, we applied systems approaches to investigate the effect of UPRmt induced by doxycycline (Dox) on transcriptome, proteome, and lipidome in two genetically divergent worm strains, named N2 and CB4856. From the integrated omics datasets, we found that Dox prolongs lifespan of both worm strains through shared and strain-specific mechanisms. Specifically, Dox strongly impacts mitochondria, upregulates defense response, and lipid metabolism, while decreasing triglycerides. We further validated that lipid genes acs-2/20 and fat-7/6 were required for Dox-induced UPRmt and longevity in N2 and CB4856 worms, respectively. Our data have translational value as they indicate that the beneficial effects of Dox-induced UPRmt on lifespan are consistent across different genetic backgrounds through different regulators.

Introduction

Mitochondria are essential organelles for numerous processes, such as energy harvesting, intermediate metabolism, autophagy, and immune response (Nunnari and Suomalainen, 2012; Quiros et al., 2016; West and Shadel, 2017). Changes in mitochondrial number, morphology, and functions not only impact cellular metabolism but also influence whole body metabolism, health, and lifespan (Andreux et al., 2013; Nunnari and Suomalainen, 2012; Vafai and Mootha, 2012). There has been increasing evidence that mitochondrial dysfunction accumulates upon aging and correlates with the development of many age-associated diseases (Sun et al., 2016). Age-associated mitochondrial impairments include decreased efficiency of oxidative phosphorylation, increases in oxidative damage, aggregation of mitochondrial proteins, alteration of mitochondrial quality control (e.g., mitophagy), accumulation of mtDNA mutations, as well as dysregulation of many aspects of mitochondrial metabolism (D'Amico et al., 2019; Jang et al., 2018; Sun et al., 2016).

Given the central role of mitochondria in health span and lifespan, mitochondria evolved an elaborate quality control system directing pleiotropic mitochondrial stress response (MSR) pathways to ensure optimal mitochondrial function and promote cell survival upon stress and aging. One of the MSR pathways is the mitochondrial unfolded protein response (UPRmt), which is well-characterized in C. elegans (Quiros et al., 2016; Shpilka and Haynes, 2018). The UPRmt is a proteotoxic stress response that senses protein-folding perturbations, which overload the capacities of the mitochondrial quality control network (Jovaisaite et al., 2014). The prototypical UPRmt is best characterized in C. elegans; when unfolded or misfolded proteins accumulate in the mitochondrial matrix in the worm, CLPP-1, a mitochondrial protease, cleaves these proteins into small peptides, which are then exported into the cytosol by the mitochondrial inner membrane peptide transporter HAF-1. On the other hand, those released peptides can impede the import of proteins into the mitochondria (Haynes et al., 2010; Naresh and Haynes, 2019; Yoneda et al., 2004). The transcription factor ATFS-1 therefore cannot be imported into the mitochondrial matrix where it is degraded by the LON protease (Nargund et al., 2012). As a result, ATFS-1 translocases and accumulates in the nucleus, forming a complex with the small ubiquitin-like protein UBL-5 and the transcription factor DVE-1. Together, this transcription complex activates the expression of nuclear-encoded protein quality components including the heat-shock proteins hsp-6 and hsp-60 to reestablish mitochondrial homeostasis (Benedetti et al., 2006; Jovaisaite et al., 2014; Nargund et al., 2012).

Paradoxically, induction of mild mitochondrial perturbations that activate the UPRmt in a moderate fashion have been shown to lead to the extension of lifespan in worms (Dillin et al., 2002; Durieux et al., 2011; Houtkooper et al., 2013), flies (Copeland et al., 2009), and mice (Houtkooper et al., 2013), suggesting that approaches to enhance UPRmt activation may potentially be useful to manage certain age-related diseases. Tetracyclines, such as doxycycline (Dox), are antibiotics that inhibit both bacterial and mitochondrial translation (Houtkooper et al., 2013). Dox has therefore been applied as a pharmacological inducer of the UPRmt (Moullan et al., 2015). Despite evidence that the UPRmt has a beneficial effect on health and aging, the molecular mechanisms that couple the Dox-mediated UPRmt with health span and lifespan extension remain to be elucidated. In addition, these effects of Dox have only been demonstrated in the reference Bristol N2 strain, and the impact of the genetic background on the protective effect of Dox has never been characterized. In this study, we have pharmacologically induced the UPRmt by Dox administration in two genetically divergent worm strains, the N2 and the Hawaii CB4856 strains (Figure 1A). By obtaining accurate quantitative data, we captured a holistic view of different molecular layers with integration of multi-omics data profiles that have been powerful in generating novel biological insights (Overmyer et al., 2021a, 2021b; Stefely et al., 2016). We show that the Dox-related protective effects are likely mediated at multiple levels of biological regulations, including transcription, translation, and lipid levels. Activation of mitochondrial stress response by Dox prolongs lifespan of both worm strains, coupling with some shared and strain-specific features at different layers of regulation. Our findings highlight a strong impact of Dox on mitochondrial functions in both worm strains. Beyond that, Dox administration also upregulated stress response and lipid metabolism, while lowering triglycerides in both strains.

Figure 1

Doxycycline (Dox) activates the mitochondrial stress response (MSR) and prolongs lifespan in two wild type C. elegans strains, N2 (Bristol), and CB4856 (Hawaii)

(A) Flowchart of the strategy to identify shared and strain-specific MSR regulators.

(B) Survival analysis of worms fed with control bacteria (E. coli HT115) culture on plates with or without Dox (15 μg/mL) by merging nine independent lifespan experiments (n = 2400). The shadow area represents the 95% confidence intervals. p values represent comparison with the controls calculated using log-rank test (∗∗: p <0.01; ∗∗∗: p <0.001; ∗∗∗∗: p <0.0001).

(C) Dox reduced oxygen consumption rate (OCR) in both worm strains. Error bars denote SEM. Statistical analysis was performed by one-way ANOVA followed by Tukey post-hoc test (∗p <0.05; ∗∗p <0.01; ∗∗∗p <0.001; ∗∗∗∗p <0.0001; N.S., not significant). Values in the figure are mean ± SEM.

(D) Principal component analysis (PCA) of transcripts/mRNAs, proteins, and lipids measured in N2 and CB4856 worms treated with or without Dox. Related to Figure S1, Tables S1, S2, and S5.

Results

Doxycycline prolongs lifespan of both N2 and CB4856 strains

To determine whether doxycycline (Dox) has a strain-dependent effect on lifespan, we cultured both N2 and CB4856 worms on plates with or w/o Dox (15 μg/mL) and measured their lifespan in nine independent experiments (Figure 1B). In basal conditions, CB4856 worms had a shorter lifespan compared to N2 worms. Upon Dox exposure, both worm strains showed an increased lifespan compared to the controls. Because Dox was known to block mitochondrial translation and attenuate respiration (Houtkooper et al., 2013), we measured the oxygen consumption rate (OCR) in N2 and CB4856 worms fed with Dox (Figure 1C). Respiration was decreased upon Dox addition in both strains. These data suggest that Dox has a universal effect on lifespan and mitochondrial function in both worm strains.

To further investigate the difference in the regulation of Dox-induced longevity in N2 and CB4856 strains, we collected worm samples either exposed to Dox or not, and extracted total mRNAs, proteins, and lipids for multi-omic analysis. Before the in-depth multi-omics analysis, we ascertained that the slightly different duration of Dox treatment applied in our lifespan (from maternal L4) and omics studies (from L1) did not impact the results of lifespan, OCR, and the UPRmt activation (Figures S1A–S1C). Next, we assessed the strain and Dox effects on each of the omic profiles (Figure 1D). Interestingly, a primary separation by Dox exposure and a second separation by strain were observed at the transcript and lipid layers (Figure 1D). However at the protein layer, the first component separated the strains and the second component separated the treatments. These data suggest that Dox showed a stronger effect on the transcriptomic and lipidomic layers and the genetic background seems to play a larger role in determining the proteomic layer.

Next, we questioned whether these alterations are commonly shared between the two wild type strains upon Dox (Figure S1D). In line with the PCA analysis (Figure 1D), up to 80% of transcripts showed similar Dox-induced changes in both worm strains, in which 2,414 and 2,021 genes were significantly upregulated and downregulated, respectively. Consistent with the PCA of the protein profiles, we detected fewer overlapping changes between the two strains upon Dox exposure, in which 127 and 205 proteins were significantly upregulated and downregulated, respectively. In particular, the majority of the altered proteins in CB4856 worms were not found to be altered in N2 worms upon Dox (Figure S1D; Table S2). In addition, we also detected a group of genes that showed reciprocal regulation between the two strains, including the 24 at transcript and 65 at the protein level. Lipid levels altered similarly in both N2 and CB4856 strains during mitochondrial stress. We then asked how many genes were significantly altered at both transcript and the corresponding protein levels (Figure S1E; Tables S1 and S2). In N2 worms, 271 transcript-protein pairs were similarly altered in the transcriptomics and proteomics profiles, in which 188 and 83 pairs were upregulated and downregulated, respectively, and 66 transcript-protein pairs were reciprocally regulated by Dox. In CB4856 worms, 196 and 181 transcript-protein pairs were unregulated and downregulated, respectively, whereas 172 pairs of mRNAs/proteins were reciprocally regulated by Dox. However, in both worm strains, the genes/proteins that did not change at both transcript and protein levels account for a substantial part of the total altered genes. This reveals that mitochondrial stress induced by Dox may have a global impact on gene expression as well as on posttranscriptional regulation. Collectively, these results suggest that mitochondrial stress induced by Dox prolongs lifespan of both worm strains and the underlying changes at the molecular levels exhibit shared and strain-specific features.

Defense response and oxidative reduction are affected by genetic variants in CB4856 compared with N2 worms

Because N2 and CB4856 at the basal levels have distinct lifespans as well as different OCR (Figures 1B and 1C), we cataloged the genetic differences and their consequences, and explored gene and protein expression differences in basal conditions (Figure 2). We compiled high impact variants present in the CB4856 strain (retrieved from the Caenorhabditis elegans Natural Diversity Resource: https://elegansvariation.org), relative to the N2 reference strain. Among these, we filtered from all the detected variants and selected the major types of consequences from homozygous variants that are of high impact, including frame-shift, start codon lost, stop codon lost, and stop codon gained for further investigations. Across the six chromosomes, most of the detected variants were distributed at the ends of each chromosome (Figure 2A; Table S3). Of these protein-coding genes, 2,755 had a frameshift (Table S3), 32 lost a start codon, 52 lost a stop codon, and 305 had a stop-gain in CB4856 relative to N2 (Figure 2C).

Figure 2

Variant analysis and gene set enrichment analyses (GSEA) determined the differences in the genetic background and gene expression between the N2 and CB4856 strains

(A) distribution of high impact variants in CB4856 per chromosome compared to the N2 reference strain. CB4856 VCF (Variant Call Format) file was retrieved from https://elegansvariation.org.

(B) Main classes of genes with high impact homozygous variants detected in the CB4856 genome are listed in the table.

(C) As an example, a circos plot of high impact variants (hard-filtered) present in the CB4856 strain compared the N2 reference strain. Analyzed variants are homozygous and with one of the following consequences: Brown: genes with start codon lost; Dark pink: genes with stop codon lost; Blue: genes with stop codon gained. The full list of variants is shown Table S3.

(D) GSEA reveals the differences between the N2 and CB4856 strains at both mRNA and protein expression levels under basal conditions. Gene sets used for enrichment analysis were retrieved from https://bio.tools/go2msig for C. elegans. Q-value: false discovery rate adjusted p values, gray dots: nonsignificant (n.s.); gene ratio: ratio of found genes within a given gene set. NES: normalized enrichment score. Related to Figure S2 and Table S3.

These genes with codon variants could be pseudogenes in CB4856 as a rapid means of adaptation (Olson, 1999); however, the persistence of these regions may suggest the presence of functional and essential genes. The genes altered by the variants are composed of members of large gene classes, including 107 fbx (F-box), 19 math (MATH), 47 clec (C-lectin), 16 cyp (cytochrome P450 family), 12 bath (BTB-MATH) genes, 28 nhr (nuclear hormone receptor), and 190 sr/str (serpentine receptor superfamily) genes (Figure 2B). F-box family, MATH, and BTB-MATH gene families encode ubiquitin-dependent proteasome adapters that were largely involved in targeting foreign proteins for proteolysis during pathogen defense (Thomas, 2006). The serpentine receptor superfamily belongs to the G-protein-coupled receptors (GPCRs) that are involved in signal transmission and are playing a vital role in controlling innate immunity against bacterial infections (Kaur and Aballay, 2020; Nagarathnam et al., 2012). The variants detected in these genes may lead to defects in pathogen defense in CB4856 worms, which partly could contribute to their shorter lifespan at basal condition compared with the N2 strain. Mammalian C-lectins are carbohydrate-binding proteins that have very narrow ligand specificity, and many C-lectins are involved in innate immune response. Although there are 264 genes encoding C-lectins in C. elegans, functions and roles of most C-lectins remain unclear (including the four detected clec genes), except for the few of them that were annotated as innate immune genes (Pees et al., 2016). The cytochrome P450 family is composed of enzymes that catalyze oxidative reactions and their substrates include lipids, exogenous, and xenobiotic chemicals (Harlow et al., 2018; Herholz et al., 2019). Nuclear hormone receptors are a family of transcription factors that often influence lipid metabolism in C. elegans (Ashrafi et al., 2003; Chinetti et al., 2000; Taubert et al., 2006). Variants found in the stop codons of these genes could impair protein function and thus affect lipid oxidation or other oxidation-reduction processes. Overall, differences in these gene families could partially explain the lifespan difference between the two wild type strains under control conditions.

To further explore the impact of the genetic background on the transcript and protein levels, we examined the top gene sets enriched in N2 or CB4856 at the control condition (Figures 2D and S2). Interestingly, we detected overall more pronounced differences between the two strains at the protein level, as compared to those detected at the transcript level. At the transcript level, the top enriched gene sets were mainly in CB4856 worms, including those involved in cell fate and transcription factor activity. In contrast, the majority of the top detected gene sets at the protein level were more pronouncedly enriched in N2 worms, including lipid metabolism, carboxylic acid metabolism, ribonucleoprotein granule, metabolic process, and biosynthesis of small molecules, fatty acids, organic acids, and monocarboxylic acids. In CB4856 worms, proteins involved in chromatin organization, histone modification, extracellular region, anchoring junction, and taxis were more enriched, compared to those detected in N2 worms.

Next, we generated enrichment maps for each worm strain to broaden the view and explore other interesting gene sets (Figure S2). Of note, we observed a presence of various gene sets involved in the nucleotide metabolism, and lipid metabolism in the N2 worms at the transcript level, compared to those in CB4856 (Figure S2A). Additionally, a large gene set cluster enriched for lipid metabolism and a cluster of genes enriched for defense response was detected at the protein level in the N2 worms (Figure S2C). Hence, in line with our variant analysis results (Figure 2), genes involved in defense response, and lipid metabolism were more enriched in N2 worms in control conditions compared with those in CB4856 worms (Figures S2A and S2C). In the CB4856 worms, we detected a trend of a very large gene set for nucleic acid metabolism, especially those associated with RNA metabolism at the mRNA level (Figure S2B). At protein level, we detected a cluster of gene sets enriched for neuron development (Figure S2D), in addition to the gene sets already mentioned above (Figure 2D). Altogether, these data provide details on strain specificity of the gene set enrichment at both transcriptomic and proteomic level and suggest that some of these differences could be attributed to genetic variants present in the CB4856 genome.

Dox upregulates lipid metabolism, proteolysis, and stress response in both strains

After observing a considerable impact of genetic variants in the CB4856 strain, we wondered whether the molecular changes upon mitochondrial stress will be affected by the different genetic backgrounds as well. We first analyzed differentially expressed genes using gene ontology (GO) analysis on the 2,414 upregulated transcripts and 127 upregulated proteins in both worm strains upon Dox exposure (Figures S1D, 3A and 3B). At the transcript level, factors involved in proteolysis, lipid metabolism and transport, stress response, innate immune response, and muscle contraction were enriched in both worm strains upon Dox (Figure 3A). At the protein level, we found fewer GO-terms compared to those at the transcript level, and most of them were related to defense response, including response to stimulus, immune response, response to oxidative stress, and response to osmotic stress (Figure 3B). Additionally, we also detected a group of GO-terms related to metabolic processes of flavonoids, organic acids, and monocarboxylic acids.

Figure 3

GO-term enrichment and GSEA revealed shared key players and strain-specific regulators of Dox-mediated longevity in N2 and CB4856

(A) GO term enrichment analysis (biological process) on the upregulated transcripts shared between N2 and CB4856 upon Dox exposure identified via David and ReviGO. The size of the dots indicates the frequency of the GO term in the underlying Gene Ontology Annotation database; the plots are color coded according to significance (Log10-transformed); level of significance increases from red to blue. GO terms belonging to the same cluster were grouped and circled in dark gray dashed lines.

(B) Go term enrichment analysis on upregulated proteins shared by N2 and CB4856 worms upon Dox.

(C and D) GSEA showed top 18 enriched gene sets in N2 (C) and CB4856 (D) worms upon Dox at both transcript and protein level. The shared gene sets between the two strains are highlighted in red. Gene sets used for enrichment analysis were retrieved from https://bio.tools/go2msig for C. elegans. Q-value: false discovery rate adjusted p-values; gene ratio: ratio of found genes within a given gene set. NES: normalized enrichment score.

(E and F) Scatter plots of differentially expressed genes (p <0.05, Log2(fold-change)>1 or <-1)) at mRNA (x-axis) and protein level (y-axis) of N2 (E) and CB4856 (F) upon Dox. Go-term enrichment analysis was then performed on each of the 8 different categories (defined based on directionality change at mRNA & protein levels). Significantly enriched GO terms are presented as a circle (radius determined by gene ratio) and edges to the related genes are shown in gray. Related to Figure S3.

Similar to the number of transcripts exclusively altered in either N2 or CB4856 worms upon Dox (Figure S1D), we also detected fewer N2-specific and more CB4856-specific enriched GO terms (Figures S3A and S3B). In the N2 strain, transcripts related to proteolysis, defense response, carbohydrate metabolism, cuticle development, and oxidation-reduction process were upregulated exclusively upon Dox (Figure S3A). In CB4856, the majority of the GO-term clusters have already been identified in the shared GO-term profile (Figure 3A), suggesting that CB4856 worms required expression of more genes in order to cope with the mitochondrial stress induced by Dox.

In addition, we sought to examine whether other effects were conferred by Dox treatment in both worm strains. We compared the transcripts of the genes involved in multiple forms of mitochondrial stress, including the integrated stress response (ISR), mitophagy, and mitochondrial dynamics (Figures S3C and S3D). Because the ISR is dispensable for the induction of the UPRmt in worms (Baker et al., 2012), and acts in complement to the UPRmt, it was not surprising to see mostly similar alterations in the expression of these ISR genes between N2 and CB4856 worms upon Dox. Of note, there was a downregulation of atf-6 (ER stress (Burkewitz et al., 2020)) and an upregulation of atf-5, a transcription factor required for regulating cytosolic translation upon UPRmt (Molenaars et al., 2020) in both strains (Figure S3C). Furthermore, genes that are involved in mitophagy and mitochondrial dynamics were mostly downregulated in both worm strains upon Dox exposure, which was different from the induction in UPRmt genes (Figure S3D).

Dox induces various responses in N2 and CB4856 worms at transcript and protein levels

Next, we determined the top enriched gene sets at mRNA and protein levels separately to assess the effects of Dox on N2 and CB4856 strains. Most of the top enriched gene sets at the mRNA level in N2 worms were upregulated, whereas the majority of the top detected gene sets at protein level were downregulated (Figure 3C). In CB4856 worms, top enriched genes were primarily downregulated at both transcript and protein levels (Figure 3D). Moreover, immune responses were upregulated in both worm strains upon Dox at the protein level, and gene sets annotated as extracellular region and defense response were upregulated in N2 and CB4856, respectively. Downregulated proteins upon Dox showed strain-specific gene set enrichments, in which proteins related to ribonucleoprotein complex and reproduction were only decreased in N2, whereas those involved in cell cycle and chromosome organization were exclusively decreased in CB4856. Although Dox is known to inhibit mitochondrial translation (Houtkooper et al., 2013), our data confirm that Dox may also attenuate cytosolic translation and this in turn leads to a downregulation of protein sets in different worm strains (D'Amico et al., 2017; Molenaars et al., 2020).

To better understand whether the changes detected at transcriptomic and proteomic levels were representative of one another, we compared the transcriptomic and proteomic profiles of N2 and CB4856 upon Dox (Figures 3E and 3F). We divided significantly altered transcripts and proteins (adjusted p-value <0.05 & an abs. fold-change>1) into eight categories based on their differences in a co-regulation analysis. An overrepresentation analysis on the detected genes was performed on each of the eight categories to determine the major altered gene sets. In N2 worms, gene sets related to molting cycle, structure organization, and cuticle development were co-upregulated at both mRNA and protein level upon Dox (Figure 3E, pink box). Genes enriched for immune response and defense response were primarily upregulated at the mRNA level (yellow box). In contrast, genes related to eggshell formation and amino sugar metabolism were co-downregulated at both mRNA and protein levels (blue box). In CB4856 worms, those that encode factors for immune response, defense response, and stress response were co-upregulated upon Dox (Figure 3F, pink box). Gene sets involved in transmembrane transport and muscle system process were only upregulated at transcriptional level (yellow box). As expected, a number of gene sets were exclusively downregulated at the protein level, including those enriched for chromosome organization, cell cycle, nuclear division, organelle fission, chromatin organization, transcription, biosynthesis of RNA, and macromolecules (purple box). Lastly, genes involved in DNA replication, cell cycle, and chromosome organization were co-downregulated at both levels (blue box). These data suggest that shared and strain-specific regulation is present at both transcriptomic and proteomic levels upon mitochondrial stress.

Dox affects mitochondrial gene expressions of both worm strains in a similar fashion

As mitochondria are the target of Dox, we then sought to investigate whether expression of specific groups of mitochondrial genes were altered by Dox in the two strains. About 138 mitochondria-related genes were significantly changed in at least one worm strain upon Dox. We first annotated the altered mRNAs based on their associated functions to assess the effects on the mitochondria (Figures 4A and S4B; Table S4). Most genes that encode the mitochondrial complex I, II, and IV, were downregulated upon Dox exposure in both worm strains, whereas those that encode complex III protein remained unchanged. Moreover, we also detected Dox-induced downregulation in transcripts involved in energy consuming processes, such as DNA replication, DNA repair, apoptosis, mitochondrial translation, and oxidation-reduction process. Many genes involved in transport, such as mitochondrial ATP translocase, ant-1.3 and ant-1.4 (Hoshino et al., 2019), and mitochondrial ion transport, sfxn-1.1, sfxn-1.2, sfxn-1.3, sfxn-1.4, and sfxn-5, were upregulated in both worm strains. For pathways that are involved in energy production, although transcripts involved in the TCA cycles were not changed in the same trend, we detected an overall upregulation in genes involved in mitochondrial fatty acid oxidation (FAO). These data suggest an enhanced energy support from FAO upon Dox.

Figure 4

Dox induced significant alterations in mitochondrial genes in both N2 and CB4856 worms

(A) Mitochondrial genes associated with different functions are significantly altered in both worm strains upon Dox at the transcript level. Outer ring: N2; inner ring: CB4856. Orange bars: upregulated genes by Dox; blue bars: downregulated genes by Dox.

(B) Changes in the expression of mitochondrial transcripts encoding for proteins localizing at different mitochondrial compartments after Dox.

(C and D) KEGG pathway enrichment analysis of upregulated and downregulated mitochondrial genes upon Dox in N2 (C) and CB4856 (D). Orange bars: upregulated genes; blue bars: downregulated genes. Pathway in red: shared between the two strains.

(E) GO enrichment analysis of the Dox-affected mitochondrial genes in both worm strains. Gene sets used for GSEA analysis were retrieved from http://ge-lab.org/gskb/for C. elegans. Related to Figure S4 and Table S4.

When we plotted the altered mitochondrial transcripts based on their localization within the mitochondria, we noticed that genes that encode proteins/enzymes located in the inter membrane space (IMS) and for membrane proteins were downregulated, whereas those at outer mitochondrial membrane (OMM) were mainly upregulated in both worm strains (Figure 4B). Genes that encode for proteins functioning in the mitochondrial matrix and inner mitochondrial membrane (IMM) were partly upregulated and downregulated upon Dox exposure. Because Dox mainly inhibits mitochondrial translation (Houtkooper et al., 2013), we then sought to assess the alterations of mitochondrial proteins upon Dox (Figures S4A and S4B). However, because of technical challenges in proteomics measurements, we were not able to collect the comprehensive changes at the protein level for these mitochondrial genes upon Dox (only 59 of the detected mitochondrial proteins were significantly altered) (Figure S4B).

To determine pathways that are significantly affected by Dox, we performed KEGG pathway and GO-term enrichment analysis on these transcripts (Figures 4C–4E). In N2 worms, degradation of BCAAs and fatty acids were significantly upregulated whereas the FOXO signaling pathway and those associated with peroxisome functions were downregulated by Dox (Figure 4C). In CB4856 worms, metabolic processes involved in BCAAs, propionate, and carbons were upregulated and those involved in base excision repair and sulfur metabolism were downregulated by Dox (Figure 4D). We then wondered if there were more similarities in the mitochondrial gene sets shared between the two strains besides BCAA degradation. Under Dox, we detected a number of upregulated GO-terms in both worm strains, including those involved in catalytic activity, oxidation-reduction process, carbon transmembrane transporter activity, cation transport, and mitochondrial inner membrane (Figure 4E). The downregulated mitochondrial genes shared between the two worm strains upon Dox were mainly enriched for DNA repair, DNA recombination, and transferase activity transferring acyl-groups (Figure 4E). These GO terms were similar to the ones that we detected in co-downregulated gene sets (Figures 3E and 3F), suggesting that the majority of downregulated proteins were related to mitochondria, which could be a consequence of inhibited mitochondrial translation by Dox.

Dox induces triglyceride degradation by lysosomal lipase to sustain energy production by fatty acid oxidation

In the abovementioned analysis performed on transcript and protein levels, we detected a large number of genes directly or indirectly linked to lipid metabolism, such as lipid transport, lipid catabolism, and mitochondrial fatty acid metabolism. To expand on the characterization of lipid-related features, we profiled lipids for N2 and CB4856 worms exposed to Dox (Figures 5 and S5; Table S5). As the primary component of lipidomic PCA was separated by treatment and the second component was separated by strains (Figure 1D), we expected to detect more shared and less strain-specific lipidomic changes in the two worm strains. In total we detected 1,572 lipids that belong to 38 lipid classes, which could be categorized further into seven main lipid categories consisting of monoacylglycerol lipids, diacylglycerol lipids, fatty acids/esters, glycerophospholipids, sphingolipids, sterol lipids, and triacylglycerol lipids (Figure S5A).

Figure 5

Lipids were significantly altered in Dox-treated worms

(A) Lipid class over-representation analysis was performed to identify top changes in the lipidomic layer in N2 and CB4856 worms upon Dox exposure. The asterisk (∗): indicates the glycerol position sn-2.

(B and C) Percentages of increased and decreased lipids in different sphingolipids, including ceramides, hydroxy ceramides, sphingomyelins, and hydroxy-sphingomyelins in N2 (B) and CB4856 (C) worms upon Dox exposure (p <0.05). A Benjamini-Hochberg corrected p <0.05 was applied to determine statistical significance.

(D and E) Changes in different lipid classes, including triglyceride (TG), diglyceride (DG), phosphatidic acid (PA), phosphatidylcholine (PC), phosphatidylethanolamine (PE), and phosphatidylglycerol (PG) in Dox-treated N2 worms compared to N2 controls (D) and in Dox-treated CB4856 compared to CB4856 worms (E).

(F) Lysosomal lipase genes are upregulated at transcript level upon Dox in both N2 and CB4856 worms.

(G)The majority of acylcarnitines was significantly increased in Dox-treated worms in a variant of their genetic background. Y-axis: number of double bonds; x-axis: number of carbon-chain length. Related to Figure S5 and Table S5.

We then determined the percentage of lipids showing significant changes in each worm strain upon Dox (Figures S5B and S5C). Although the overall lipid changes were similar between the two strains, we detected more decreased lipid species exclusively in the N2 worms (126 reduced lipids), compared to those only lowered in CB4856 worms (31 reduced lipids) upon Dox (Figure S1D). In addition, about half of the lowered lipids by Dox are TGs, of which 143 and 105 TGs were decreased in N2 and CB4856 worms, respectively (Table S5). To better understand the overall impact on the lipidomic profiles upon Dox, we performed overrepresentation analysis on the lipid classes to identify the top changed lipid species upon Dox in each worm strain (Figure 5A). Among the top altered lipid classes, most of them increased in N2 worms upon Dox, including LPE, PC, PA, LPA, DG, PLC, PE, CE, PS, DG, SM, CL, and PI. Different from those detected in N2 worms, fewer classes of lipids were increased in CB4856 worms upon Dox, including PE, PA, LPE, CE, PC, LPA, PS, and PI (Figure 5A). Intriguingly, three classes of lipids, Cer, TG, and PG were lowered in both worms upon Dox (Figure 5A).

Notably, the sphingolipids, such as ceramides, hydroxy ceramides, sphingomyelins, and hydroxy-sphingomyelins were significantly altered upon Dox; many of them decreased in both strains except few hydroxy species which were increased exclusively in CB4856 worms (Figures 5B and 5C). Low levels of sphingolipids, especially low ceramides, have been shown to couple with increased lifespan in C. elegans (Cutler et al., 2014; Huang et al., 2014; Kim et al., 2016) and with caloric restriction in mouse liver (Green et al., 2017). It is therefore not surprising that these sphingolipids were decreased in both worm strains upon Dox exposure, which universally increased lifespan.

In addition, we also observed some shared features in other lipid species, including increased levels of DG, PA, PC, and PE, and decreased levels of TG and PG (Figure 5A). Although some TGs with three double bonds slightly increased in both worm strains and alkyl-diacylglycerol (TG[O]) increased in N2, the majority of the significantly altered TGs were decreased upon Dox (Figures 5D and 5E). In line with this, we also detected a marked increase of DGs with >2 double bonds and decrease of those with ≤1 double bond, irrespective of the carbon chain lengths in Dox-treated worms. This could be because of either enhanced lipolysis of TGs or reduced TG synthesis. To assess this hypothesis, we checked the changes in the expression of genes encode the enzymes that are involved in both TG synthesis and breakdown pathways. The transcripts of genes involved in TG synthesis, i.e, diglycerol acyltransferase dgat-2 and mboa-2 remained unchanged in both worms upon Dox (Table S1), whereas transcripts of genes involved in TG breakdown, belonging to both adipocyte triglyceride lipase ATGL family (Srinivasan, 2015) and lysosomal lipase family, lipl-1, was upregulated in both worms upon Dox (Figure 5F). The transcript of another gene that belongs to the ATGL family, lipl-3, was upregulated exclusively in CB4856 worms upon Dox. In addition, transcripts of many other lysosomal lipase genes, including lipl-2, lipl-4, lipl-5, and lipl-6, were higher in CB4856 worms, and lipl-4, and lipl-6 in N2 worms upon Dox (Figure 5F). However, in the proteomics analysis, lipl-2, lipl-5, and lipl-7 abundances were not significantly changed upon Dox (Table S2), suggesting that the regulation of the ATGL family enzymes occurred mainly at the transcriptional level. The link between genome and lipid profiles is still largely unknown and merits further investigation (Linke et al., 2020). In addition, we also observed a significant increase of PA, PC, and PE with >2 double bonds, irrespective of the carbon chain lengths in Dox-treated worms (Figures 5D and 5E), suggesting an upregulation of either lipolysis or polyunsaturated fatty acid (PUFA) biosynthesis. We then analyzed the fatty acid and acylcarnitine profiles (Figures 5G S5D and S5E). In both N2 and CB4856 worms, fatty acids remained largely unchanged upon Dox, except for increases of two PUFAs, C18:3 and C20:4, and a decrease of saturated fatty acid (SFA) C22:0 in N2, and increase of C18:0, C18:1, and C15:1 in CB4856 worms (Figure S5E). Strikingly, we detected a major change in the acylcarnitines upon Dox in both strains (Figure 5G), with levels of medium-chain and long-chain acylcarnitines being largely increased in response to mitochondrial stress. Additionally, we determined fat storage by lipid staining of both worm strains upon Dox treatment. As expected, we detected a tendency towards decreased TG storage in the lipid droplets, confirming the overall decreased TG levels measured in Dox treated worms by lipidomics analysis (Figures S5F and S5G). In line with the detection of upregulated FAO transcripts in both worms upon Dox (Figures 4A–4D), we therefore speculated that FAO was activated during mitochondrial stress induced by Dox, to produce energy using fatty acids that were released from TGs by lipolysis. Consistent with this, CPT-1 protein was significantly upregulated in both strains with Dox treatment (Table S2).

Lipid genes acs-2/20 and PUFA synthesis genes fat-7/6 are essential in mediating Dox-induced longevity in N2 and CB4856 worms, respectively

To validate the biological relevance of the candidate regulators identified from the multi-Omics analysis, we selected a group of genes that are involved in different branches of lipid metabolism and performed lifespan measurement. As we have speculated the potential involvement of lipolysis, FAO, and PUFA synthesis, we selected lysosomal acid lipase lipl-4 (Folick et al., 2015) and FAO gene acs-2 (acyl-CoA synthase) (Ramachandran et al., 2019), PUFA synthesis genes fat-2, fat-6, and fat-7 (Gao et al., 2017; Watts and Browse, 2002), as the top candidates. In addition, as we found cuticle development as the most essential gene/protein sets only in Dox-treated N2 (Figures 3E and S3A), we also included acs-20, a gene encodes acyl-CoA synthetase that plays an essential role in cuticle structures and surface barrier function (Kage-Nakadai et al., 2010). Interestingly, RNAi of acs-2 and acs-20 exclusively attenuated Dox-mediated longevity in N2 worms, whereas RNAi of fat-7 and fat-6 attenuated the lifespan extension in CB4856 worms (Figures 6A and 6B). Next, we took advantage of the UPRmt reporter hsp-6p::gfp strain (on an N2 background), to test the role of the candidate genes in UPRmt activation upon Dox exposure. We found that RNAi of acs-2 and acs-20 (but not fat-7), strongly suppressed UPRmt activation in these N2 worms (Figure 6C). To verify whether fat-7 and fat-6 are required for UPRmt activation in CB4856 worms, we collected worm samples and examined changes in UPRmt genes by qRT-PCR. Of note, although RNAi of fat-7/6 did not affect the expression of UPRmt genes involved in mitochondrial function and metabolic processes (hsp-6 and gpd-2), RNAi of these two lipid genes significantly suppressed the expression of immune response genes lys-2 and clec-4 (Figure 6D). This result confirms our omics analysis showing that the immune response was one of the most essential gene sets regulated in the CB4856 strain upon Dox treatment. Taken together, Dox-induced longevity is hence most likely regulated in a strain-specific manner.

Figure 6

Investigation of lipid genes involved in Dox-mediated longevity

(A) Lifespan measurement in N2 (A) and CB4856 (B) worms with RNAi of fat-7/6, fat-2 (PUFA synthesis), acs-2 (FAO), acs-20 (acyl-CoA synthetase) with or without Dox, respectively.

(C) RNAi of acs-2 and acs-20, but not fat-7, attenuates the GFP expression of hsp-6p::gfp worms (in the background of N2) induced by Dox. Scale bar: 10 μm.

(D) RNAi of fat-7/6 suppressed the expression of immune response genes in CB4856. Statistical analysis of lifespan assays was performed by log-rank test (∗∗∗∗p <0.0001; ∗∗∗p <0.001; ∗∗p <0.01; ∗p <0.05; n.s., not significant). Statistical analysis of qPCR results was performed by one-way ANOVA followed by Tukey post-hoc test (∗∗∗∗p <0.0001; ∗∗∗p <0.001; ∗∗p <0.01; ∗p <0.05; n.s., not significant). Values in the figure are mean ± SEM. Related to Figure S6.

Discussion

In this study, we employed a comprehensive systems approach to explore effects of Dox-induced mitochondrial stress on multiple layers of biology and elucidate the molecular changes in two wild type C. elegans strains, i.e., N2 (Bristol) and CB4856 (Hawaii). This included the integrated analyses of transcripts, proteins, and lipids collected in these two worm strains either not exposed or exposed to Dox. Through this approach, we uncovered shared and strain-specific effects of Dox, and important pathways/regulators that are involved in Dox-induced mitochondrial stress response and lifespan extension. Previously, we have shown the robust activation of UPRmt by administering Dox to N2 worms (Houtkooper et al., 2013; Li et al., 2021). It remains, however, unknown whether Dox could also activate the UPRmt and prolong lifespan in other C. elegans strains. Although CB4856 worms are genetically divergent and have a shorter lifespan compared with the N2 strain under normal conditions (Banse et al., 2019; Dingley et al., 2014; Doroszuk et al., 2009; Thompson et al., 2015), we observed a comparable effect of Dox on the lifespan and respiratory rate in both worm strains. These data suggest a general and beneficial effect of Dox in terms of lifespan extension. Given the inherent complexity of certain phenotypes, such as lifespan, we assessed whether Dox exhibited shared or strain specific effects at different molecular levels. Based on the PCA of each molecular layer, we noticed that both strain and treatment effects could explain the majority of the difference at each layer of biological regulation.

To better understand the influence of the genetic background in response to mitochondrial stress, we first explored the genetic variants between the two strains in the control condition. In total, we cataloged 3,658 variants with potential disruptive impact on the protein translation that are distributed on six chromosomes of CB4856 compared to the N2 reference. This number is comparable to the scope of variations detected between these two strains in prior studies (Kim et al., 2019; Thompson et al., 2015). Genes with these homozygous variants mostly belonged to the F-box, MATH, BTB-MATH, C-lectins, Cytochrome P450, serpentine receptor, and nuclear hormone receptor families. Variants that were detected in these families of genes may lead to impaired defense response, oxidative reduction, and lipid metabolism in CB4856 worms. This may partly contribute to their shorter lifespan compared to N2 worms because all these pathways are closely involved in mediation of lifespan (Goudeau et al., 2011; Kaur and Aballay, 2020). At the protein levels, we detected more enriched gene sets in N2 worms for lipid metabolism and defense response, supporting our findings in the variant analysis. In addition, we detected more significantly changed mRNAs/proteins in our study compared with a previous study (Kamkina et al., 2016).

Differences in the genetic backgrounds of the two worm strains not only influenced molecular changes at the unstressed condition but also played a vital role upon mitochondrial stress. When analyzing either transcriptomic or proteomic data, we observed a number of shared and strain-specific regulations of different pathways and gene sets. Genes involved in lipid metabolism, proteolysis, and stress response were upregulated in both worms upon Dox at the transcript level, which is consistent with the results of previous studies activating UPRmt via other inducers in N2 worms (Li et al., 2021; Liu et al., 2014, 2020; Pellegrino et al., 2014; Sorrentino et al., 2017). Of note, after we integrated transcriptome and proteome profiles, we observed more strain-specific alterations upon Dox. In N2 worms, immune and defense response-associated genes were exclusively upregulated at transcriptomic level, whereas these genes were upregulated at both transcriptomic and proteomic levels in CB4856 worms. Upregulation of these genes at both transcriptional and translational levels might be a part of a compensatory mechanism balancing the defects in gene expression caused by the large number of disruptive variants present in the CB4856 compared to the N2 worms. In addition, downregulated proteins also showed to be strain-specific, in which ribonucleoprotein complex and reproduction-related proteins were exclusively decreased in N2, and cell cycle and chromosome organization associated proteins were exclusively decreased in CB4856. Our results validated that blocking mitochondrial translation reduces cytosolic translation and leads to a detection of reduced levels of proteins in both N2 and CB4856 worms as a consequence (D'Amico et al., 2017; Molenaars et al., 2020).

As the primary molecular target of Dox is the mitochondrial ribosome and mitochondrial translation (Houtkooper et al., 2013), we took a deeper look at the differentially expressed mitochondrial genes at both mRNA level and protein level according to their function and localization within the mitochondria. At the transcript level, we observed a common effect of Dox on genes involved in various mitochondrial functions. Overall energy consuming pathways were attenuated, including mitochondrial translation (Houtkooper et al., 2013), consistent with reported observations in worms with mitochondrial stress induced by other stressors (Li et al., 2021; Melber and Haynes, 2018; Nargund et al., 2015). Moreover, we observed a metabolic remodeling upon Dox in both worm strains, in which FAO genes were significantly upregulated in both worm strains, suggesting a switch to lipid breakdown to fuel energy production. This speculation was further supported by the upregulated transcripts that belong to the ATGL and lysosomal lipase families and by the lowered TGs as well as the increased acylcarnitines levels in both worm strains upon Dox. Evidence for enhanced FAO in mitochondrial-stressed worms has been documented by others (Baruah et al., 2014; Liu et al., 2020; Tharyan et al., 2020). Moreover, increased expression of fatty acid β-oxidation gene acs-2 (fatty acyl-CoA synthetase (Van Gilst et al., 2005)) in mitochondrial stressed worms was previously shown to be atfs-1-dependent (Wu et al., 2018), suggesting a specific requirement of FAO upon mitochondrial stress. Although upregulated lipid metabolism at the transcript level was not mirrored at the proteomic level, the lipidomic layers still underwent significant alterations in worms exposed to Dox. Besides TGs, sphingolipids, and PGs were significantly decreased in both strains. This was accompanied with an increase of DG, PA, PC, and PE.

Next, we validated the involvement of lipid genes in Dox-mediated longevity. We focused on lipid genes involved in PUFA synthesis (fat-2, fat-6, and fat-7), mitochondrial FAO (acs-2), cuticle development (acs-20), and lysosomal lipase (lipl-4). Except cuticle development (acs-20), which was found to be the most essential gene/protein exclusive to Dox-treated N2, the other lipid pathways were affected in both worm strains upon Dox exposure. We hence performed lifespan measurements to determine the roles of these candidates in each worm strain. Interestingly, knockdown of acs-2 and acs-20 exclusively suppressed Dox-mediated lifespan extension in N2 worms, confirming an essential role of FAO and cuticle development (epidermal innate immune response (Kage-Nakadai et al., 2010; Taffoni and Pujol, 2015)). In contrast, knockdown of fat-7 and fat-6 attenuated the expression of immune response genes as well as lifespan extension in CB4856 worms, confirming that stearoyl-CoA desaturase is necessary for the Dox-induced activation of immune response (Anderson et al., 2019) and lifespan extension in CB4856 worms.

Collectively, our comprehensive multi-omics systems approach revealed shared and strain-specific pathways that are important for mitochondrial stress-induced longevity in two wild type worm strains. Effects of Dox-induced mitochondrial stress at different omics layers in the two different worm strains showed more shared features, suggesting universal benefits of Dox. Our data provide new mechanistic insights in Dox-induced mitochondrial stress, which is accompanied by a metabolic rewiring that favors lipid metabolism and TG degradation. This highlights the importance of lipid metabolism in mitochondrial stress-mediated longevity in genetically divergent worm strains.

Limitations of the study

The current study has used only two genetically divergent strains to study the effects of Dox on biological and molecular layers, precluding an in-depth analysis of genetic factors that contribute to the differences. A study in a C. elegans genetic reference population (GRP), such as the recombinant inbred lines (RILs) or the recombinant intercross advanced inbred lines (RIAILs) will hence be better suited to identify the alleles underlying the subtle variations in physiological and molecular traits linked with the UPRmt. We also were unable to identify the tissue contributing to the differences between the two strains.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resource and reagents should be directed to and will be fulfilled by the lead contact, Johan Auwerx (admin.auwerx@epfl.ch).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

C. elegans and bacterial feeding strains

The Bristol strain (N2) and Hawaii strain (CB4856) were used as the wild-type strains, and SJ4100 [zcIs13(hsp-6p::GFP)] were obtained from the Caenorhabditis Genetics Center (CGC; Minneapolis, MN). E. coli OP50 and HT115 strains were also obtained from the CGC. RNAi clones against fat-2, fat-6, fat-7, acs-2, acs-20, and lipl-4 were obtained from the Ahringer and Vidal libraries and verified by sequencing before use. Worms were cultured and maintained at 20°C and fed with E. coli OP50 on Nematode Growth Media (NGM) plates unless otherwise indicated.

Method details

Lifespan measurements

Worm lifespan was performed as described previously (Mouchiroud et al., 2011). In brief, 5-10 L4 worms of each strain were transferred onto RNAi plates (containing 2 mM IPTG and 25 mg/mL carbenicillin) or RNAi plates containing 15 μg/mL Doxycycline (Dox, Cat. D9891, Sigma) seeded with either E. coli HT115 bacteria or RNAi bacteria. After the progenies reached the last larval stage L4, worms were then transferred onto RNAi (or RNAi +15 μg/mL Dox) plates containing 10 μM 5FU. Approximately 80 worms were used for each condition and scored every other day. Nine independent lifespan experiments were performed.

Oxygen consumption rate analysis

Oxygen consumption rate (OCR) was measured with the Seahorse XF96 (Seahorse Bioscience) as described (Koopman et al., 2016). In brief, ∼100 N2 and CB4856 worms were cultured on plates with or without Dox, and washed off with sterile M9 buffer after worms reached day 1 adulthood. After washing with M9 buffer for three times, worms were then transferred in a 96-well Seahorse plate and OCR was measured for six times. Mitochondrial OCR was measured for each condition.

Imaging of worms with UPRmt activation by Dox

RNAi bacteria were cultured overnight in lysogeny broth (LB) medium containing 100 μg/mL ampicillin at 37°C. Next, the bacteria were seeded onto either RNAi plates or RNAi plates supplemented with Dox. L4/young adult worms were picked onto the RNAi bacteria-seeded plates and cultured at 20°C until their progenies reached young adult stage. 6 - 10 progenies were then randomly picked in a drop of 10 mM tetramisole (Cat. T1512, Sigma) and then aligned on an empty NGM plate. Fluorescent images, with the same exposure time for each condition, were captured using a Nikon SMZ1000 microscope.

Worm fixation and lipid content staining

Worms were cultured and maintained at 20°C and fed with E. coli OP50 on NGM plates. Eggs were obtained by alkaline hypochlorite treatment of gravid adults and ∼200 eggs were transferred onto the experimental plates. At L4 larval stage, worms were collected, washed twice with 1 x PBS and then suspended in 120 μL of PBS to which an equal volume of 2X MRWB buffer (160 mM KCl, 40 mM NaCl, 4 mM EGTA, 30 mM PIPES pH 7.4, 1 mM spermidine, 0.4 mM spermine, 2% paraformaldehyde, 0.2% beta-mercaptoethanol) was added for 1 h. After fixation, the worms were freeze-thawed for three times with dry ice and in a 37°C water bath, followed by centrifuging for 1 min at 14,000 x g and then washing once in PBS buffer. Prior to Oil Red O staining, worms were re-suspended and dehydrated in 60% isopropanol for 15 min 250 μL of 60% Oil Red O stain was added to each sample, and samples were incubated overnight at 4°C. Worms were washed twice in 60% isopropanol solution after Oil Red O staining (Li et al., 2018). Worm images were acquired with the Leica DM5500 Upright Microscope. The lipid staining been quantified with the open-source image analysis software package FIJI/ImageJ from https://imagej.net/Fiji.

Sample collection for RNA-seq, proteomics, and lipidomics analyses

Worms of each strain were cultured on plates seeded with E. coli OP50, then eggs were obtained by alkaline hypochlorite treatment of gravid adults. A synchronized L1 population was obtained by incubating the egg suspension in M9 butter overnight at room temperature. Approximately 2000 L1 worms were transferred onto plates with or without 15 μg/mL Dox seeded with E. coli HT115. L4 worms were harvested after 50 h by three times of washing with M9 buffer. Tubes containing worm pellets were immediately submerged in liquid nitrogen for snap freezing and stored at −80°C until use. Three biological replicates were collected for each condition.

RNA extraction, RNA-seq data analysis and qRT-PCR

On the day of the extraction, 1 mL of TriPure Isolation Reagent was added to each tube. The samples were then frozen and thawed quickly eight times with liquid nitrogen and 37°C water bath to rupture cell membranes. RNA was then extracted by using a column-based kit from Macherey-Nagel. RNA-seq was performed by BGI with the BGISEQ-500 platform. FastQC was used to verify the quality of the mapping (de Sena Brandine and Smith, 2019). No low-quality reads were present and no trimming was needed. Alignment was performed against worm genome (WBcel235 ce11 primary assembly and Ensembl release 89 annotation) following the STAR (version 2.73a) manual guidelines (Dobin et al., 2013). The obtained STAR gene-counts for each alignment were analyzed for differentially expressed genes using the R package edgeR (version 3.24.3) using a generalized linear model (Robinson et al., 2010). The trimmed mean of M values (TMM) method was chosen to normalize the counts and the Cox-Reid common dispersion method for computing the dispersion parameter (tag wise dispersion was also computed). A principal component analysis was also generated to explore the primary variation in the data (Lê et al., 2008; Risso et al., 2014).

For qRT-PCR, worms were cultured and total RNA was extracted as described for the RNA-seq sample preparation. Synthesis of cDNA was conducted from total RNA using the Reverse Transcription Kit (Qiagen, 205314). qPCR was carried out using the Light Cycler 480 SYBR green I Master kit (Roche, Cat. 04887352001). The primers used for qPCR are listed in the key resources table. Primers for worm pmp-3 were used as housekeeping controls.

Protein extraction and proteomics analysis

Proteins were denatured, extracted with vigorously shaking and digested with LysC and trypsin. Peptides were analyzed on our liquid chromatography setup with the Fusion Lumos Mass Spectrometer (Hebert et al., 2014). In brief, pellets containing ∼2,000 worms were resuspended in 900 μL lysis buffer A (6M Guanidine hydrochloride, 100 mM Tris pH 8.0), and one large stainless-steel bead was added to each 2 mL tube. Tubes were shaken in a bead miller (Retsch) at 30 Hz continuously for 10 min. Then the lysates were placed in a sand bath and heated at 95°C for 10 min. Protein concentrations of lysates were determined using protein BCA assay (Pierce). To extract proteins, 100 μL of lysate were added to 900 μL 100% methanol. Tubes were then spun for 7 min at 15,000 rpm, supernatants were discarded, and protein pellets were allowed to briefly air dry. 50 μL lysis buffer B (8 M urea, 100 mM Tris, 40 mM TCEP, and 10 mM 2-chloroacetamide) were added to the pellets, and the tubes were vigorously vortexed for 10 min. LysC (Wako) was added in 1:50 enzyme: protein ratio, and samples were digested overnight at ambient temperature on a rocker (Fisher Scientific). On the next day, digestion mixtures were diluted down to 1.5 M urea with addition of 50 mM Tris, pH 8.0. Sequencing grade trypsin (Promega) was added in 1:50 enzyme: protein ratio for further 3 h digestion. Samples were then acidified to pH of ∼2 by addition of 10% TFA and desalted using Strata SPE columns (Phenomenex). Desalted peptides were lyophilized to dryness in a SpeeVac (Thermo), then resuspended in 50 μL 0.2% formic acid. Final peptide concentrations were determined using Peptide Colorimetric Assay (Pierce). 1 μg peptides were injected onto an in-house high pressure packed capillary column (Shishkova et al., 2018), housed at 55°C, and separated using Dionex Ultimate 3000 nano HPLC system (Thermo Fisher) over 120 min gradient at 315 nL/min. Mobile phase A consisted of 0.2% formic acid in water, and mobile phase B consisted of 0.2% formic acid in 70% acetonitrile. Eluting peptides were electro-sprayed into and analyzed on Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher). Orbitrap MS1 scans were collected at resolution of 240,000 at 200 m/z with an AGC target of 1 × 106 ions and maximum injection time set to 50 ms. Precursors were isolated in a quadrupole with the isolation window of 0.7 Th. Tandem MS scans were collected in the ion trap using rapid scan rate, AGC target of 1 × 104 ions, HCD fragmentation with NCE of 25, and dynamic exclusion of 20 s. RAW files were searched using MaxQuant (Cox et al., 2014) against Uniprot database of C. elegans, containing canonical sequences and isoforms. Unless specified, default settings were used. Protein abundances were quantified using label-free quantification with count of 1 and match-between-runs enabled. MS2 tolerance was set to 0.27 Da. Similar differential analysis was performed as for the transcriptome layer using the R package edgeR (version 3.24.3) with a generalized linear model, similar normalization and principal component analysis methods (Robinson et al., 2010). The proteomic data files were shared using the in-house web-based data visualization tool Argonaut (Brademan et al., 2020).

Lipid extraction and lipidomics data analyses

Around ∼2500 worms were collected in a 2 mL tube and the following amounts of internal standards dissolved in 1:1 (v/v) methanol:chloroform were added to each sample: BMP(14:0)2 (0.2 nmol), CL(14:0)4 (0.1 nmol), LPA(14:0) (0.1 nmol), LPC(14:0) (0.5 nmol), LPE(14:0) (0.1 nmol), LPG(14:0) (0.02 nmol), PA(14:0)2 (0.5 nmol), PC(14:0)2 (0.2 nmol), PE(14:0)2 (0.5 nmol), PG(14:0)2 (0.1 nmol), PS(14:0)2 (5 nmol), ceramide phosphocholine SM(d18:1/12:0) (2 nmol) (Avanti Polar Lipid). A steel bead and 1:1 (v/v) methanol:chloroform was added to each sample to a volume of 1.5 mL. Samples were homogenized using a TissueLyser II (Qiagen) for 5 min at 30 Hz and centrifuged for 10 min at 20,000 g. The supernatant was transferred to a 4 mL glass vial and evaporated under a stream of nitrogen at 45°C. The lipids were dissolved in 150 μL of 1:1 (v/v) chloroform:methanol and lipidomics analysis was performed as described (Molenaars et al., 2021).

In brief, the HPLC system consisted of a Thermo Fisher Scientific Ultimate 3000 binary UPLC coupled to a Q Exactive Plus Orbitrap mass spectrometer using Nitrogen as the nebulizing gas. The column temperature was maintained at 25°C. For normal-phase separation, 2 μL lipid extract was injected onto a Phenomenex® LUNA silica, 250 × 2 mm, 5 μm 100 Åand for reverse phase 5 μL of each sample was injected onto a Waters HSS T3 column (150 × 2.1 mm, 1.8 μm particle size). A Q Exactive Plus Orbitrap (Thermo Scientific) mass spectrometer was used in the negative and positive electrospray ionization mode. In both ionization modes, mass spectra of the lipid species were acquired by continuous scanning over the range m/z 150-2000 with a resolution of 280,000 full width at half maximum (FWHM). For quantification of the data, lipids were normalized to corresponding internal standards according to lipid class, as well as total protein content in samples, determined using a PierceTM BCA Protein Assay Kit (Herzog et al., 2016).

Similar differential analysis was performed as for the transcriptome and proteome, using the Bioconductor package limma version 3.42.2 (Ritchie et al., 2015), with a generalized linear model, similar normalization and principal component analysis methods. Results of the statistical tests were corrected for multiple testing using the Benjamini-Hochberg method.

Gene set enrichment analysis (GSEA)

We used the clusterProfiler R package (Yu et al., 2012) to conduct GSEA analysis (Figures 2D, and 3C–3D). Gene sets used for enrichment analysis were retrieved from https://bio.tools/go2msig for C. elegans. We used a minimum gene set size of 10, a maximum gene set size of 500, and performed 10,000 permutations. For the transcriptome data, we used a gene list that is ordered by log2(Fold-changes) from the differential expression analysis. In the proteomic analysis, the gene name for each protein were retrieved and the gene list was then ordered based on fold change values of the protein quantity obtained from the differential analysis of the proteomic level.

Transcriptomic/proteomic GO enrichment analysis

All significantly changing genes/proteins (adjusted p value <0.05 and an absolute fold change >1) were split into 8 groups based on the combination of direction of the fold change, includes genes up & down regulated significantly in a single analysis (i.e. transcriptome OR proteomics). An overrepresentation analysis using GO gene sets was performed on each of the eight groups to determine the main – if any – gene sets changing. GO term enrichment analyses shown in Figures 3 and S3 were performed with David and ReviGo (Huang da et al., 2009; Supek et al., 2011).

Variant analysis

Variants for the CB4856 strain were retrieved from the “Caenorhabditis elegans Natural Diversity Resource” (https://elegansvariation.org/). Soft filtered variant file was retrieved from http://storage.googleapis.com/elegansvariation.org/releases/20200815/variation/WI.20200815.soft-filter.vcf.gz. Hard filtered variant file was retrieved from http://storage.googleapis.com/elegansvariation.org/releases/20200815/variation/WI.20200815.hard-filter.vcf.gz. These variants were then filtered for the CB4856 strain keeping only 1/1 variants with high impact consequence.

Figure generation

Circos genomics tool (http://circos.ca/) was used to generate the circular representations (Figures 2C, 4A, 4B, and S4A). ReviGO was used to generate clustering of enrichment analysis (Figures 3A, 3B, S3A, and S3B). Results from ReviGO were downloaded as R scripts and images were further edited in Illustrator. The emaplots (Figures S2A–S2D) were generated using the R package clusterProfiler's emapplot function. Enrichment analysis was performed using the "EnrichGO" function of the ClusterProfiler R package.GO term enrichment analyses shown in Figures 3 and S3 were performed with David and ReviGo (Huang da et al., 2009; Supek et al., 2011).

Quantification and statistical analysis

No statistical methods were applied to pre-determine worm sample size. Comparison between more than two groups were assessed by using one-way ANOVA test. Prism 8 (GraphPad Software) was used for statistical analysis of all the bench experiments. Variability in panels is given as the s.e.m. All p <0.05 were considered to be significant. (∗∗∗∗p <0.0001; ∗∗∗p <0.001; ∗∗p <0.01; ∗p <0.05; n.s., not significant. For lifespan, imaging, and OCR measurement in worms, sample size was determined based on the known variability of the experiments. All experiments were done non-blinded and repeated at least twice.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Bacterial and virus strains
Escherichia coli: OP50	Caenorhabditis Genetics Center	RRID:WB-STRAIN:OP50
Escherichia coli: HT115 (DE3)	Caenorhabditis Genetics Center	RRID:WB-STRAIN:HT115(DE3)
fat-2 RNAi	Vidal library	10012-D5
fat-6 RNAi	Ahringer	IV-6G22
fat-7 RNAi	Vidal	11061-F8
acs-2 RNAi	Ahringer	V-10A16
acs-20 RNAi	Ahringer	IV-6H23
lipl-4 RNAi	Vidal	11069-E7
Chemicals, peptides, and recombinant proteins
Doxycycline	Sigma-Aldrich	Cat# D9891
5-Fluorouracil (5-FU)	Sigma-Aldrich	Cat# F6627
Ampicillin sodium salt	Sigma-Aldrich	Cat# A9518
Carbenicillin disodium salt	Sigma-Aldrich	Cat# C1389
IPTG	AppliChem	Cat# A1008,0005
TriPure Isolation Reagent	Roche	Cat# 11667165001
EGTA	Sigma-Aldrich	Cat# E-3889
PIPES	Sigma-Aldrich	Cat# 80635
KCl	Fisher Chemicals	Cat# P/4240/60
Paraformaldehyde	Acros Organics	Cat# 416785000
Beta-mercaptoethanol	Acros Organics	Cat# 125470010
PBS	Gibco	Cat# 10010-015
Oil red O	Sigma-Aldrich	Cat# O0625
Critical commercial assays
NucleoSpin RNA, Mini kit for RNA purification	Macherey-Nagel	Cat# 740955.250
Seahorse XFe96 Extracellular Flux Assay kit	Agilent	Cat# 102416-100
RNA using the Reverse Transcription Kit	Qiagen	Cat# 205314
LightCycler 480 SYBR Green I Master kit	Roche	Cat# 04887352001
Deposited data
C. elegans RNA-seq data	This paper	Table S1, GEO: GSE159228
C. elegans proteomics data	This paper	Table S2
C. elegans lipidomics data	This paper	Table S5
Experimental models: Organisms/strains
C. elegans: N2 Bristol	Caenorhabditis Genetics Center(CGC); https://cbs.umn.edu/cgc/home	CGC:10570RRID:WB-STRAIN:N2_(ancestral)
C. elegans: CB4856 Hawaii	Caenorhabditis Genetics Center(CGC); https://cbs.umn.edu/cgc/home	CGC:7525RRID:WB-STRAIN:CB4856
C elegans: SJ4100 [zcIs13(hsp-6p::GFP)]	Caenorhabditis Genetics Center(CGC); https://cbs.umn.edu/cgc/home	CGC: 23223RRID:WB-STRAIN:SJ4100
Oligonucleotides
hsp-6Fw: AGAGCCAAGTTCGAGCAGAT	Sigma-Aldrich	N/A
hsp-6Rv: TCTTGAACAGTGGCTTGCAC	Sigma-Aldrich	N/A
lys-2Fw: GCTGCGATTGCAAATACCGA	Sigma-Aldrich	N/A
lys-2Rv: TGGCAGTTGGATTGTTTGGC	Sigma-Aldrich	N/A
clec-4Fw: GGGATATGGAGCGACACTGG	Sigma-Aldrich	N/A
clec-4Rv: TCGCAAATCTTCTGGCCCTT	Sigma-Aldrich	N/A
gpd-2Fw: AAGGCCAACGCTCACTTG AA	Sigma-Aldrich	N/A
gpd-2Rv: GGTTGACTCCGACGACGA AC	Sigma-Aldrich	N/A
pmp-3Fw: GTTCCCGTGTTCATCACTCAT	Sigma-Aldrich	N/A
pmp-3Rv: ACACCGTCGAGAAGCTGTAGA	Sigma-Aldrich	N/A
Software and algorithms
GraphPad Prism v8	GraphPad Software, Inc.	https://www.graphpad.com/scientificsoftware/prism/
FIJI/ImageJ	FIJI/ImageJ	https://imagej.net/Fiji
David	Huang da et al., 2009	https://david.ncifcrf.gov/tools.jsp
REVIGO	Supek et al., 2011	http://revigo.irb.hr/
R	The R Foundation	https://www.r-project.org/
Variant analysis	Caenorhabditis elegans Natural Diversity Resource	https://elegansvariation.orgSoft filtered variant file was retrieved from: http://storage.googleapis.com/elegansvariation.org/releases/20200815/variation/WI.20200815.soft-filter.vcf.gz.Hard filtered variant file was retrieved from:http://storage.googleapis.com/elegansvariation.org/releases/20200815/variation/WI.20200815.hard-filter.vcf.gz.
Adobe Illustrator 2021	Adobe	https://www.adobe.com/products/illustrator.html
Caenorhabditis elegans GSEA GO Gene Set	GO2MSIG	https://bio.tools/go2msig
Circos	Circos	Krzywinski, Martin, et al. "Circos: an information aesthetic for comparative genomics." Genome research 19.9 (2009): 1639-1645.doi:10.1101/gr.092759.109http://circos.ca/
reshape2	reshape2	https://cran.r-project.org/web/packages/reshape2/index.html
edgeR	edgeR	https://bioconductor.org/packages/release/bioc/html/edgeR.html
refGenome	refGenome	https://cran.r-hub.io/web/packages/refGenome/index.html
Fastqcr	Fastqcr	https://cran.r-hub.io/web/packages/fastqcr/index.html
Limma	Limma	https://bioconductor.org/packages/release/bioc/html/limma.html
FactoMineR	FactoMineR	https://cran.r-project.org/web/packages/FactoMineR/index.html
Plotly	Plotly	https://cran.r-project.org/web/packages/plotly/index.html
Cowplot	Cowplot	https://cran.r-project.org/web/packages/cowplot/index.html
RColorBrewer	RColorBrewer	https://cran.r-project.org/web/packages/RColorBrewer/index.html
Openxlsx	Openxlsx	https://cran.r-project.org/web/packages/openxlsx/index.html
Dplyr	Dplyr	https://CRAN.R-project.org/package=dplyr
xlsx	Xlsx	https://CRAN.R-project.org/package=xlsx
clusterProfiler	clusterProfiler	https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html
enrichplot	Enrichplot	https://bioconductor.org/packages/release/bioc/html/enrichplot.html
ggrepel	Ggrepel	https://cran.r-project.org/web/packages/ggrepel/index.html
ggplot2	ggplot2	https://cran.r-project.org/web/packages/ggplot2/index.html
UpSetR	UpSetR	https://cran.r-project.org/web/packages/UpSetR/index.html
lme4	lme4	https://cran.r-project.org/web/packages/lme4/index.html
GenomicFeatures	GenomicFeatures	https://bioconductor.org/packages/release/bioc/html/GenomicFeatures.html
org.Ce.eg.db	org.Ce.eg.db	https://bioconductor.org/packages/release/data/annotation/html/org.Ce.eg.db.html
BSgenome.Celegans.UCSC.ce11	BSgenome.Celegans.UCSC.ce11	https://bioconductor.org/packages/release/data/annotation/html/BSgenome.Celegans.UCSC.ce11.html