To analyze genotype-dependent responses to CR, we took advantage of natural genetic variation across 46 diploid wild yeast isolates of S. cerevisiae (Saccharomyces cerevisiae) and six other closely related budding yeast species; Saccharomyces kluyveri, Saccharomyces bayanus, Saccharomyces kudriavzevii, Saccharomyces paradoxus, Saccharomyces castellii, and Saccharomyces mikatae. We analyzed their RLS phenotype under both control (high glucose, 2%) and CR (low glucose, 0.05%) conditions.

Within the S. cerevisiae population, we found ∼10-fold median RLS variation under the control condition, and RLS analysis under the CR condition revealed a 7-fold variation among them (Figure 1A, Supplementary Material S1). While 26 strains displayed significant (Wilcoxon rank sum test, adjusted p ≤ 0.05) CR response (either decreased or increased) in median RLS, the remaining 20 strains did not respond to CR (no significant changes in median RLS) (Figure 1A, Supplementary Material S1). Among the 26 CR responding strains, 11 of them showed various degrees of median RLS extension, while the remaining 15 displayed a wide range of decreased median RLS (Figure 1A, Supplementary Material S1). For example, under the CR condition YJM978 strain showed a decreased median RLS (50% decrease), whereas CR caused a 75% median RLS increase of the Y9 strain compared to the control condition (Figure 1B, Supplementary Material S1). Consistent with the previous report (Zou et al., 2020), the control diploid laboratory WT strains BY4743, which is a derivative of the original S288c isolate, also increased median RLS by 12% (Wilcoxon rank sum test, adjusted p = 0.009). On the other hand, the original S288c strain did not respond to CR (Figure 1B), indicating adaptive changes to laboratory conditions resulted in an alteration to CR response in this strain background.

Next, to understand whether the longevity effect of CR also varies across different species, we also analyzed the RLS of six different budding yeast species of Saccharomyces genus; S. kluyveri, S. bayanus, S. kudriavzevii, S. paradoxus, S. castellii, and S. mikatae. Interestingly, only two of them, S. mikatae and S. kluyveri showed significant (Wilcoxon rank sum test, adjusted p ≤ 0.05) lifespan (median RLS) extension when subjected to the CR (Figure 1C, Supplementary Material S1).

Overall, this data showed considerable variation in lifespan phenotype under CR conditions among the genetically diverse natural isolates of the same species and the closely related species of the same genus. This data suggests that CR at the 0.05% glucose level does not promote lifespan extension for most of the wild-derived yeast strains and species, and this applies to other non-laboratory adapted model organisms.

Gene expression variation has been suggested to play a significant role in adaptive evolution. Existing research also highlights the potential influence of gene expression levels on various phenotypic traits and its plasticity (Gilad et al., 2006; Vu et al., 2015), such as changes in lifespan (Whitaker et al., 2014; Liu et al., 2023).

Accordingly, we examined whether the comparison of transcriptome profiles of these strains, obtained under the high glucose condition, can reveal molecular signatures that can predict the CR response type. First, we examined whether evolutionary relationships based on the gene expression variation are associated with the CR response types; we constructed gene expression phylograms for S. cerevisiae strains using a distance matrix of one minus Spearman correlation coefficients (Brawand et al., 2011) and based on normalized reads (Kaya et al., 2021). We found that CR response types are mainly branched through early adaptation of strain-specific life history trajectories at gene expression level, and this adaptation is mostly acquired independently (Figure 2A). Subsequently, principal component analysis (PCA) was conducted, and the result did not reveal a distinct segregation pattern, except for a few outlier strains that were separated by PC1 and PC2. Intriguingly, this separation resembled the phylogenetic relationship, with the cumulative effect of the first three principal components accounting for approximately 35% of the total variance in gene expression (Figure 2B). We then employed a different dimensionality reduction method, partial least squares discriminant analysis (PLS-DA), particularly suited for distinguishing different groups (Please see methods). PLS-DA analyses revealed distinct clusters for positively responding strains (Figure 2C). The model evaluation metrics provided further validation, with an R2Y value of 0.926 indicating that the model accounts for 92.6% of the variability in these three groups and a Q2Y value of 0.766, suggesting a predictive accuracy of 76.6%.

Next, to identify molecular signatures and genetic regulators of CR response types, we compared the gene expression pattern between CR responding strains (increased median RLS) to non-responding strains (no significant change) and negative responding strains (decreased median RLS). Under high glucose conditions, 222 genes differentially expressed (DEGs), including 146 genes that had significantly reduced expression and 76 genes that had increased expression in the CR responding group in comparison to the non-responding group (adjusted p ≤ 0.05 and log2-fold >0.5) (Figure 3A, Supplementary Material S2). Similar analysis revealed 176 DEGs, including 50 genes with significantly reduced expression and 120 genes with increased expression in the CR responding group compared to the negative-responding group (Figure 3B, Supplementary Material S2).

To further investigate the molecular patterns associated with the CR response type, we analyzed the commonalities and differences in gene expression between two groups (positively responding vs. non-responding and positively responding vs. negatively responding). Among the downregulated DEGs identified from both comparisons, 32 genes were in common. The GO term associated with these commonly downregulated genes in positively responding strains was enriched in cellular amino acid biosynthesis (Figure 3C). Among the upregulated DEGs, there were 34 genes in common. Functional enrichments of these commonly upregulated DEGs resulted in a single GO term associated with mitochondrial translation (Figure 3D).

We then performed GO enrichment and KEGG pathway analysis for the uniquely altered genes in positively responding strains compared to the non-responding or negatively responding strains alone. The 114 unique downregulated DEGs, resulting from the comparison of positively responding versus non-responding strains, revealed that positively responding strains have decreased lipid and amino acid biosynthesis and metabolism compared to the non-responding strains. (Figures 4A, B). A similar result was also obtained for 24 downregulated DEGs from positive versus negative responding strains comparison (Figures 4C, D). The upregulated 42 unique DEGs resulted from a comparison of positively responding versus non-responding strains enriched in the ribosome (both mitochondrial and cytosolic) and mRNA surveillance pathways. Additionally, GO term enrichment indicated an increase in iron metabolism in positively responding strains (Figures 4A, B). The GO terms for 86 upregulated DEGs resulted from positively responding versus negatively responding strains comparison were enriched in mitochondrial translation, ATP synthesis coupled proton transport, ATP biosynthetic process, and mitochondrial respiratory chain complex IV assembly (Figures 4C, D). Overall, these results suggest that responding strains maintain higher translation under high glucose conditions and are uniquely adapted to regulate mitochondrial function for energy production compared to the non-responding and negatively responding strains. In addition, our data suggests that negative responding strains are characterized by high amino acid biosynthesis and metabolisms. On the other hand, non-responding strains might also be compromised for iron and copper homeostasis. Further analyses of KEGG pathways and GO term enrichment from each comparison revealed decreased fatty acid and branched-chain amino acid synthesis in positively responding strains. At the same time, these two processes are upregulated in other groups under high glucose conditions (Figures 4C, D; Supplementary Figure S1). Overall, we showed a phenotype-specific gene expression that, in comparison to the non- or negatively responding strains, specific group genes were selectively up or downregulated in positively responding strains under high glucose conditions (Figure 5 and Supplementary Material S2).

Next, to further understand the relationship between gene expression variation and CR-mediated lifespan variation, we performed, phylogenetic regression by generalized least squares (PGLS), as we described previously (Kaya et al., 2015a; Kaya et al., 2021). Our analyses further showed that expression levels of some of these genes were significantly correlated with CR response variability across all strains (Figure 6 and Supplementary Material S2). For example, we identified 118 transcripts with significant correlation with median RLS (adjusted p ≤ 0.01; 48 with positive correlation and 70 with negative correlation) (Supplementary Material S2). Among the top hits with positive correlation were a major mitochondrial D-lactate dehydrogenase (DLD1-oxidizes D-lactate to pyruvate), homeodomain-containing protein and putative transcription factor (TOS8), and SKI complex-associated protein (SKA1-involved in involved in 3′-5′ degradation of long 3′UTR-containing mRNA) (Figure 6A). The top hits with negative correlation included the genes coding for DNA end-binding protein required for nonhomologous end joining (NEJ1), GTPase-activating protein for Gpa1p (SST2), and Subunit of G protein involved in pheromone response (GPA1) (Figure 6A). The data indicates a role for GPA1 function in regulating lifespan under CR condition, which plays a role in mating-related nuclear migration and karyogamy also involved in inositol lipid-mediated signaling and regulation of MAPK export from the nucleus.

In agreement with the cross-comparison transcriptome data, the KEGG pathway analysis of these positively correlated genes revealed that strains with increased oxidative phosphorylation activity and amino acid and fatty acid degradation pathways show effective CR response (Figure 6B). Contrarily, negatively correlating genes enriched in MAPK signaling pathway, branched-chain amino acid biosynthesis, mismatched repair, and cell cycle (Figure 6B).

In conclusion, the comprehensively analyzed transcriptome data across diverse genetic backgrounds and the lifespan variation under CR conditions uncover the molecular factors associated with CR response types.

We also searched for metabolites whose abundances identified under high glucose condition are associated with CR response type. The metabolome represents a snapshot of relevant biological processes downstream of the proteome, and it has been widely used for characterizing aging-regulated metabolic pathways and biomarkers for age-associated diseases (Cheng et al., 2015; Laye et al., 2015; Ma et al., 2015). Among the 166 metabolites that we examined, none of them showed significant differences at adjusted p-value cut-off (p ≤ 0.05) in both responding versus non-responding and responding versus negatively responding comparisons (Supplementary Material S2). At the non-adjusted p-value cut-off (p ≤ 0.01), we found responding strains with decreased abundance of leucine and S-Adenosyl-homocysteine and increased abundance of inositol monophosphate in comparison to both non- and negatively responding groups. Three metabolites, nicotinic acid, mevalonate, and 1-metyhylnicotinamide were specifically characterized by decreased abundance in responding strains in comparison to the non-responding group. Oxidized glutathione and phenylalanine were specifically characterized by decreased abundance in positively responding strains in comparison to the negatively responding group (Figures 7A, B and Supplementary Material S2).

Next, we examined the correlation between metabolite abundance and CR-mediated lifespan variations. Our phylogenetic regression (PGLS) analyses identified a single metabolite, n-formylmethionine (fMet), with a significant (adjusted p = 0.041) positive correlation (Figure 7C, Supplementary Material S2). Methionine (Met) derivate fMet initiates protein synthesis in mitochondria. Our data showed that fMet abundance was not coupled to Met abundance since there was no significant correlation between Met and median RLS (Figure 7D, Supplementary Material S2). Given that both transcriptomic and metabolomic data were obtained under high glucose conditions, our regression analysis was not aimed at predicting CR response type. However, it gives an association between CR-mediated lifespan extension rates and molecule abundance which can be utilized as a predictive marker.

Next, to further investigate CR-mediated strain-specific transcriptional changes, we selected two CR responding (DBVPG1106, DBVPG1373) and non-responding (BC187, DBVPG6765) strains and collected cells that were grown under YPD-CR condition and subjected them transcriptomics analyses. We observed a distinct transcriptional response between responsive and non-responsive strains under CR conditions (Figures 8A, B). Our calculation of distances between strains based on the biological coefficient of variation (BCV) (McCarthy et al., 2012) of gene expression revealed two separate clusters. The CR regimen caused transcriptional changes that segregated the non-responding strains, mainly BCV1 and BCV2. Conversely, CR responding strains only separated along BCV2 (Figure 8B). The segregation pattern indicates unique transcriptional changes among the responding and non-responding groups. To analyze DEGs, we analyzed CR-mediated changes compared to the controls (high glucose) for each group. We found that CR caused more robust expression changes in the non-responding group. Overall, there were 2,185 DEGs (1213 downregulated and 972 upregulated) significantly altered (adjusted p ≤ 0.01, log2-fold change ≥1) in the non-responding group when subjected to CR. On the other hand, we identified only 422 DEGs (230 downregulated and 192 upregulated) significantly altered in the positively responding group when subjected to CR (Figure 8C, Supplementary Material S3). Among them, 134 downregulated and 110 upregulated genes were shared between the responding and non-responding groups (Figure 8C, Supplementary Material S3). These commonly altered genes are enriched in cytoplasmic translation for downregulated and Gluconeogenesis, TCA cycle, Biosynthesis of nucleotide sugars, Tryptophan metabolisms and Longevity regulating pathways, and Fatty Acid degradation for upregulated genes (Supplementary Material S3).

Then, we analyzed response type-specific DEGs for each group. For the non-responding group, we found that CR induces gene expression changes that upregulate the pathways associated with Carbon metabolism, TCA cycle, Oxidative phosphorylation, Propanoate metabolisms, Branched Chain amino acid and Lysine degradation, Longevity regulating pathway, Biosynthesis of cofactors, and Autophagy (Figure 9A). Among the downregulated pathways were Ribosome biogenesis, Amino acid biosynthesis and metabolism, and Nucleotide metabolisms (Figure 9B). Analysis of DEGs in CR responding strains revealed that CR caused downregulation of Ribosomal genes and Tryptophan metabolisms and upregulation of MAPK signaling, TCA cycle, Amino acid metabolic pathways, and Biosynthesis of secondary metabolites. This data suggests that CR selectively acts on specific sets of genes to alter targeted pathways to extend lifespan in responding strains (Figures 10A, B). For example, responding strains are characterized by higher mitochondrial respiration and translation under high glucose conditions (Figures 4B, D). Accordingly, the altered genes were only related to decreased ribosomal biogenesis, and there were no further alterations for genes associated with mitochondrial function in responding strains under CR conditions. However, CR-induced expression of genes alters mitochondrial function and decreases ribosomal biogenesis and translation for non-responding strains. Considering that the non-responding strains are characterized by decreased ribosomal biogenesis and translation under high glucose conditions (Figures 4B, D), nonselective repression of ribosomal biogenesis genes under CR conditions might further decrease translation. Thus, the imbalance between mitochondrial and cytosolic translation might be a confounding factor for preventing CR-mediated lifespan extension in non-responding strains, as suggested previously (Molenaars et al., 2020).

Next, to understand the regulatory network under CR conditions, we identified core Transcription Factors (TF) for up- and downregulated genes of responding and non-responding strains. In responding strains, increased expressions of genes are regulated by several transcription factors, including SKN7, RPN4, GCN4, HSF1, MSN2/4, YAP1, and MIG1 (Figure 11A). Many of these transcription factors regulated genes associated with environmental stress response genes (ESR) (Gasch et al., 2000). In fact, 53 of the 192 upregulated genes in responding strains belonged to the ESR group (Supplementary File S3). Furthermore, our data suggests that many of the ESR genes and other upregulated genes are regulated by GCN4 (Figure 11B). On the other hand, decreased expressions of genes are regulated by HMS1, CRZ1, and UME1 (Figure 11C) in the responding strains.

For the nonresponding strains, the increased expressions of genes were regulated by MET32, MSN2, HCM1, HAP4, AFT2, and AFT1 under CR condition (Supplementary Figure S2A). The data suggests that increased expressions of respiratory genes were regulated by HAP4 and HCM1 in non-responding strains. ATF1 and ATF2 regulate iron homeostasis, while MET32 is involved in transcriptional regulation of the methionine biosynthetic and other sulfur metabolic genes. Among the top TFs associated with downregulated genes were SFP1, CST6, HSF1, RAP1, and FHL1 (Supplementary Figure S2B). Many of these TFs are involved in the regulation of ribosomal gene expression and might be associated with decreased translation in non-responding strains.

We further dissected the effect of mitochondrial respiration on CR response by comparing median RLS variation between YPD-CR and YPG (3% glycerol) conditions. A respiratory growth substrate, glycerol can extend RLS (Kaeberlein et al., 2005; Burtner et al., 2009) by a switch from fermentation to respiration; however, the exact mechanism of glycerol-mediated lifespan extension is unclear. We previously observed significant differences in RLS across these strains on YPG condition (Kaya et al., 2021). Interestingly, a comparison of RLS phenotypes of these strains from YPG and YPD-CR revealed a significant positive correlation (R ² = 0.41, p = 6.18 × 10⁻⁶) (Supplementary Figure S3) that strains responding to CR also respond to YPG positively or vice versa (Figure 12, Supplementary Material S1). This data further confirms mitochondria associated with overlapping mechanisms mediating lifespan under glycerol and CR conditions.

Next, to investigate the mitochondrial role in CR-mediated lifespan extension further, we eliminated mtDNA in two responding strains to isolate respiratory deficient cells (Woo and Poyton, 2009; Kaya et al., 2015b). Then, we analyzed their RLS under control and CR conditions. Elimination of mtDNA caused a decrease in median RLS in all three of them under 2% glucose condition (Figure 13A, Supplementary Material S1). These results indicate mitochondria-specific metabolic adaptation in wild isolates since the elimination of mtDNA has shown mixed effect on median RLS in laboratory-adapted strains (Kaeberlein and Guarente, 2005; Woo and Poyton, 2009). Next, we found that eliminating respiratory deficiency blocked the CR effect on lifespan extension (Figure 13A). Together, these data pinpoint that CR-mediated lifespan extension requires functional mitochondria in wild yeast isolates.

Finally, we wanted to test whether increased GCN4 expression alone can promote lifespan extension in selected CR responding strains (Figure 13B). Our transcriptome data showed that CR induces GCN4 expression and alters GCN4-mediated pathways in CR responding strains (Figure 12). We further found that GCN4 overexpression alone could induce lifespan extension (Figure 13B, Supplementary Material S1). Our data represent unique adaptive changes involving GCN4 function to mitochondria-mediated metabolic and molecular adaptations to the environment as the main factor for lifespan regulation under CR conditions.

Diploid wild isolates of S. cerevisiae were obtained from the Sanger Institute (Carter et al., 2009), and the other budding yeast species were gifted by Mark Johnston (Cliften et al., 2001). We also reported detailed information and lifespan phenotypes of these strains analyzed under high glucose (2% glucose) and glycerol (3% glycerol) conditions (Kaya et al., 2015a; Kaya et al., 2021). For the expression of GCN4, we used a modified p426ADH plasmid by inserting a hygromycin (HYG) cassette along with its promoter and terminator at the XbaI restriction site (Kaya et al., 2021). We omitted the 5′UTR sequence of the GCN4 as described previously (Mittal et al., 2017). Yeast transformation was performed using the standard lithium acetate method. Finally, to isolate rho ⁰ strains, cells were cultured in YPD medium, supplemented with 10 μg/mL and ethidium bromide (EtBr), and incubated at room temperature with agitation for approximately 24 h. This procedure was repeated three times, and after the third growth cycle, cells were diluted (1:100) in water and plated on YPD to obtain single colonies. After that, several individual colonies were selected to test their growth ability on YPG (respiratory carbon source) plates. Colonies that were unable to grow on YPG were selected as rho ⁰ .

RLS assay method and lifespan phenotypes of these strains on YEP (yeast extract, peptone) medium supplemented with 2% glucose or 3% glycerol was described in detail in our previous publications (Kaya et al., 2015a; Kaya et al., 2021; Oz et al., 2022). We modified the YEP plates by supplementing them with 0.05% glucose (YPD-CR) to determine the RLS phenotypes on the CR medium. For each natural isolate, at least 30 individual mother cells were analyzed. Each assay also included the BY4743 strain as a technical control. For RLS analysis of wild isolates harboring expression plasmids for GCN4, several individual colonies were picked up from selection medium (HYG) after transformation, and YPD medium supplemented with 200 μg/mL HYG was used for RLS determination of these cells. Survival analysis and Gompertz modeling were performed using the survival (https://cran.r-project.org/web/packages/survival/index.html) and flexsurv (https://cran.r-project.org/web/packages/flexsurv/index.html) packages in R, respectively.

The RNAseq and metabolomics procedure and the data analyses for cells collected on YPD (2% glucose) were described previously (Kaya et al., 2021). We used the same procedure to collect cells and perform RNAseq analyses for Saccharomyces species and for the cells collected don YPD-CR (0.05% glucose) medium. Data analyses were also done by using similar methodologies for consistency.

To assess the impact of experimental treatments on gene expression, we conducted a Principal Component Analysis (PCA) utilizing the “FactoMineR” package in R (Lê et al., 2008). The PLS-DA analysis was performed by using the function ‘opls’ in R package ‘ropls’ with default parameters (https://bioconductor.org/packages/devel/bioc/vignettes/ropls/inst/doc/ropls-vignette.html#43_Partial_least-squares:_PLS_and_PLS-DA). Additionally, we computed the Spearman correlation coefficient between samples based on their gene expression values. The distance matrix among samples was derived from the 1−correlation coefficient and employed the “hclust” function in R for hierarchical clustering using the “ward.D2” method. For the differential expression analysis, we utilized the “EdgeR” package in R (Robinson et al., 2009) to perform comparisons between positively responding versus non- and negatively responding groups. Briefly, we initially applied the “calcNormFactors” function to calculate the normalization factor for the raw counts of gene expressions by using the value of the total expression of samples. Subsequently, we used the “voomWithQualityWeights” function to estimate the weight of each sample and transformed the expression matrix using voom. Given that there are multiple strains within each group, with each strain having three biological replicates, treating biological repetitions as independent observations would introduce noise into the results. To address this, we employed the “duplicateCorrelation” function to calculate the correlation among biological replicates. We then reapplied the “voomWithQualityWeights” function to re-estimate the sample weights and transform the expression matrix, considering the correlation between biological replicates. Finally, we utilized the “lmFit” function to identify the significant differential expressed gene by setting each strain as a block and incorporating the correlation between biological replicates. The adjustment of the p-value was performed using the Benjamini–Hochberg (BH) method. Gene met adjusted p-value <0.05 and | log2 (Fold change) | > 0.5 as the significant genes. To identify endophenotypes (transcripts and metabolites) correlating with CR-mediated lifespan variation across wild isolates, we performed phylogenetic regression using the generalized least squares method, as we described previously (Kaya et al., 2015a; Kaya et al., 2021).

We performed gene enrichment analysis using R packages “clusterProfiler” (Yu et al., 2012). We selected the items contained in the Gene ontology (GO) and the KEGG databases to analyze up- and downregulated genes for each group. The Benjamini–Hochberg (BH) method was used for the p-value adjustment, and the terms with adjusted p-value <0.05 were selected. For the KEGG database-based pathway enrichments, terms with adjusted p values <0.5 were considered. To identify TFs associated with DEGs identified under CR conditions, we utilized the yeastrack portal (Teixeira et al., 2023), which provides a set of queries to predict transcription regulation networks in yeast from various -omics data.