Reduced sphingolipid biosynthesis modulates proteostasis networks to enhance longevity

Four hours of myriocin treatment induce robust transcriptional changes

To examine transcriptional changes induced by Myr treatment, we diluted stationary phase prototrophic BY4741 cells (50–60% of cells synchronized for first two cell cycles over a 6 h time frame [7]) into fresh culture medium (Time 0), with and without Myr treatment. Samples for analysis of mRNA abundance by RNA seq were taken at 1 h intervals over a 6 h time course (Figure 1A). The normalized RNA seq data for the time course contained transcripts, expressed as transcripts per million (TPM), mapped to 6198 genes, 5169 of which were uniquely annotated and of sufficient signal intensity for subsequent analysis (Supplementary Table 1, Filter tab). Filtered data were analyzed by two-way ANOVA with a statistical cutoff of p = ≤ 0.01 to give a set of 4964 significant genes (Figure 1B). These were further sorted into drug (D), time (T), both time and drug (T&D) and Interactions (I) (Figure 1B, Euhler diagram), and assigned a p-value. Genes found to be significant in any of these three categories are referred to as differentially expressed (DEG).

DEGs that change up or down over time were identified by using a ‘template assignment tool’ (MATERIALS AND METHODS, example in Figure 2A, center panel). Most DEGs (99.86%) were assigned to a template. For DEGs significant by Time (T), Drug (D), or both (T&D), post hoc analysis categorized each DEG into one of 11 different temporal patterns plus one with no temporal effect (Figure 1C, Supplementary Table 1 - Pattern Graph tab can be used to plot the pattern of any gene). Within each temporal pattern (other than ’00’) a gene was placed in one of two horizontal reflections (Figure 1C, columns ‘a’ or ‘b’). Three types of Drug effects are shown in each pattern (Figure 1C, Y-axis: light blue, black or orange curves) along with the number of genes in each category (see Legend for Figure 1C for more details). These 11 patterns explained the majority (77%) of all DEGs. Results of DEGs in significantly enriched patterns are shown as a heatmap and a graph for Vehicle and Myriocin treatments (Figure 1D, see Legend for details and pattern labels). A graph of any gene transcript across the 1–6 h time-frame can be plotted using the Graph Reporter tab in Supplementary Table 1. The most significant GO terms enriched in these patterns are shown in Table 1 and a complete list of genes in each GO term is presented in Supplementary Table 2 (Webgestalt patterns combined tab). These GO terms primarily represent major processes, cellular compartments or molecular functions needed to drive cell growth and division such as chromatin organization, transcription regulation, mitotic cell cycle and biochemical pathways.

In contrast to the Time x Drug sector, the 1009 genes in the Time sector of the Euler diagram fall primarily into three patterns (Figure 1C and 1D, T01b(↔), T10A (↔) and T10b(↔)). Genes in these patterns are not significantly changed by Myr treatment but are significantly changed with time and serve as an example of time-related changes. All of these patterns will require further effort to determine their significance, if any, in Myr-induced longevity.

The Drug sector of the Euler diagram (Figure 1C and 1D, pattern D(↑)) contains 329 up-regulated genes (Table 1 and Supplementary Table 2). Because of our experimental design, we cannot be sure if some genes in this sector are driving the changes associated with time that we see, but they could be. The most enriched GO terms relate to chromosomes, and transcription (Table 1), indicating a strong response to Myr treatment independent of time in culture.

Next, we analyzed the 1570 genes in the Interaction group for temporal patterns responding to Myr treatment (Figure 2A). Five patterns contained more than 100 DEGs (Figure 2A, center panel, numeric results for this grid are available in Supplementary Table 1, tab called ‘I intersections’). The right-hand graph in Figure 2A indicates the number of up- and down-regulated genes in these 5 patterns. Five templates for these patterns (Figure 2B, left-side graphs) along with heatmaps of representative genes in each pattern for Vehicle and Myriocin treated samples are shown in Figure 2B (center panels). The right-most panels in Figure 2B represent average Z-scored gene expression profiles for the 5 patterns examined. Both the heatmaps and the data plots on the right side of Figure 2A reveal that Myr treatment enhances transcription at the 5 and 6 hour time points. This observation was tested by examining all genes found to be significant by the two-way ANOVA’s in the interaction term. For this test we plotted the log 2-fold changes of Myr treatment/vehicle treatment at each time point and color-coded the results if they were significantly down-regulated (blue) or up-regulated (orange) (p ≤ 0.05, pairwise Fisher’s LSD). As predicted, the predominant effect of Myr was to upregulate genes at the 5 and 6 h time points (Figure 2C). Functional overexpression analysis of significant DEGs at these two time points (called: I: 5 and 6(↑)) revealed a more generalized picture of the biological processes and cellular compartments responding to Myr (Figure 2C, right panel, all genes listed in Supplementary Table 2). As we discuss in detail below, up-regulation of genes starting between the 4th and 5th hour occurs when the majority of cells enter their second cell division cycle [7]. Thus, these data explain much of the effect of time and drug components within the interaction gene set.

We recently reported that Myr treatment has a notable effect on the size of most amino acid pools which remain significantly smaller in drug-treated cells starting before the 1 h time point, and remaining smaller over the six hour time course studied [7]. To determine if transcription has roles in lowering and maintaining smaller amino acid pools, we performed functional enrichment analysis on the 1570 genes in the Interaction group (Figure 2A) because this group captures a large fraction of genes responding differently over time with Myr treatment. We hypothesized that down-regulation of genes would imply a contribution to lower amino acid pools while up-regulation would imply the opposite effect or an attempt to restore pools to the size found in vehicle-treated cells. The KEGG pathway term ‘sce01100:Metabolic pathways’, containing 314 genes, is highly enriched in this gene set (p = 1,6E-27, FDR 1.54E-25) (Supplementary Table 2, INTERACTION group and 314 genes Metabolic pathways tabs). Importantly, there is strong enrichment for genes in other KEGG pathways including tryptophan, methionine and branched chain amino acid metabolism. Additionally, within the Interaction group the term ‘GO:0003333~amino acid transmembrane transport,’ with 16 genes, is enriched. A distinguishing feature of most genes in these three pathways plus the transporters is up-regulation in Myr-treated cells with many genes peaking at 5 and 6 h. These data suggest that Myr-treated cells are attempting to increase amino acid uptake and de novo synthesis as a way to increase amino acid pool levels, however, we consider other interpretations in the Discussion.

Correlation analysis of transcriptomics data identifies clues to amino acid metabolism

As another approach to understand effects of Myr on the transcriptome, we analyzed transcript data using the R program package of Weighted Gene Correlation Network Analysis, WGCNA [28], to identify clusters (modules) of highly correlated genes showing a similar response to Myr over the 1–6 h time course. A gene dendrogram produced by average linkage of hierarchical clustering revealed many modules with the seven most correlated having from 18 to 2490 genes (Supplementary Table 3, column L in Cluster Analysis tab, see also the Module Eigengenes tab). We then examined the relationship of these gene modules to free amino acid pools over the same 1–6 hr time course as an alternative way to determine if transcription played a role in lowering or maintaining smaller amino acid pools in drug-treated cells as we had previously found [7]. Results of this analysis are presented in Figure 3A as a heatmap plot. For amino acids (abbreviation shown on the X-axis), the numbers in a column represent the correlation p-value and the correlation (in parentheses) between a module eigengene (ME, Y-axis) expression value and an amino acid level. Here a ME corresponds to the first principal component of the expression matrix of the module. The ME can be considered the most representative gene in a module and, as used here, it represents the average effect of Myr on genes compared to the untreated drug control (Supplementary Table 3, Module Eigengenes tab). The heatmap shows that the Brown and Green modules with 1126 and 144 genes, respectively, have the greatest negative correlation to Myr treatment. That is, Myr maintains low pool levels but enhances the level of a transcript at one or more time points. The values in the ‘InQ’ column of the heatmap represent the correlation p-value and correlation value (parentheses) between the module eigengene expression and incubation time course (1–6 h). For these association calculations, the control (no drug) samples were set as 0 and Myr-treated samples were set as 1. Thus, the MEturquoise row has a strong correlation with time (Figure 3A, column inQ) as shown graphically in Figure 3F. The values in the 'myriocin' column represent the association between the module eigengene expression and the absence/presence of Myr and they indicate that the Green and Brown MEs are the most negatively correlated MEs.

As another way to represent the effects of Myr treatment on MEs and to detect interactions between MEs, we performed network analysis and visualized the results by using Cytoscape. One such network, representing the relationship between genes in the Brown and Green modules, includes the Turquoise module because genes interact with multiple genes in both the Brown and Green modules (Figure 3B). We further analyzed these three MEs by plotting the ME for each time point with and without Myr treatment (Figure 3C). For example, the ME representing the Brown module reveals Myr induces transcription starting at 4 h (Figure 3C) whereas there is no increase in amino acid pool size in drug-treated cells [7]. This opposite effect of Myr on transcription and amino acid pools corresponds to the negative correlation represented by the green shading in Figure 3A. The highest-ranking GO term in the Brown module is Regulation of Transcription (Figure 3D) where drug treatment enhances transcript levels particularly after 4 h, consistent with the ANOVA analysis (Figure 2B and 2C). A different negative correlation pattern is seen for the ME representing the small Green module where genes are up-regulated by Myr across the 1–6 h time frame compared to vehicle control but gene expression drops in both control and drug-treated samples starting at 2 h (Figure 3E). The main GO term in the Green module is Translation (Biological Process, p = 1.7E-88) along with Ribosome Assembly and related processes responding to drug-induced slowing of protein synthesis and growth rate [7]. Lastly, the large Turquoise module with 2490 genes captures Myr-induced transcriptional events involving processes or pathways or cell components each with more than 200 genes (Figure 3F, 3G). Interestingly, transcripts represented in the Turquoise module increase across the 1–6 h time frame with drug-treated cells always having higher transcript levels, substantiating the ANOVA analysis represented in Figure 1D (patterns T10a(↔), T&D04a(↑)., T&D05a(↑)), Figure 2B (patterns I05a(04a(↑)), I10a(04a(↑)) and I10a(10a(↑))) and Figure 2C. Functional enrichment analysis for all MEs in presented in Supplementary Table 3.

Deubiquitination of Mup1 impairs Myr-enhanced longevity

A drawback of WGCNA and pathway and pattern analyses is a reduced ability to identify significant biological features involving smaller numbers of genes or to identify genes that do not fit a specific pattern across the 1–6 h time-frame. To circumvent these limitations and to discover transcript changes with potential roles in lowering amino acid pools or novel roles in longevity enhancement, we examined genes with the highest possible significance (1E-16) in the Drug, Time and Interactions columns of data from the two-way ANOVA analysis of mRNA levels. This approach identified only 16 genes out of the 4964 genes (Supplementary Table 1, Filter tab, columns M, N and O). The UBI4 gene, which is stress-induced and encodes for ubiquitin, captured our attention because stresses of various types have known roles in aging and longevity and we previously identified increased stress responses in Myr-treated cells [27, 29]. A defining feature of UBI4 transcript abundance in our studies is a rapid increase at the 5–6 h time period (Figure 4A, statistical significance of Area Under Curve (AUC) 95% CI (difference of the means): −189.081 to −137.609, p-value = 0.00006). The path of the UBI4 transcript level corresponds to the 5–6 h time period in which Myr has its most significant effect on transcription (Figure 2B and 2C) and it falls within Interaction pattern I: 5 and 6(↑) (Figure 2C and Supplementary Table 2, Webgestalt patterns combined tab). Proteolysis is a highly enriched GO term in this pattern of gene expression (Figure 2C, table on the right). The potential biological significance of the UBI4 transcript pattern may relate to our previous analyses of Mup1-pHluorin trafficking. We found that the fluorescent signal decreases in the PM more rapidly during the 4–7 h time-frame in Myr-treated cells, in a dose-dependent manner, which requires ubiquitin conjugation for endocytosis to occur [7].

To determine if ubiquitination has roles in Myr-enhanced chronological lifespan (CLS), we examined untreated ubi4Δ cells and found a faster loss of viability compared to untreated BY4741 cells (Figure 4B, 50% survival day 3 vs. day 5, respectively). Additionally, the CLS of ubi4Δ cells treated with Myr does not increase and remains the same as untreated cells. Untreated wild-type (WT, BY4741) cells show 50% survival at day 5–6 which is extended to day 11 by Myr treatment (AUC, 95% CI: −6.828 to −2.840, p-value = 0.0025). From these data we conclude that UBI4 is required for Myr to enhance CLS.

Since ubiquitin has many functions, we studied its role in ubiquitin-mediated endocytosis of Mup1 using the same strains as used in our published analyses of endocytosis [7]. Specifically, CLS was assayed using cells with chromosomal MUP1 tagged with pHluorin-UL36 (catalytic active or catalytic dead, CD). UL36 is a viral deubiquitinase (DUB) that can reverse localized ubiquitin conjugation activity, thus protecting the fusion protein from ubiquitin-mediated degradation or trafficking events [30]. CLS is strongly enhanced by Myr treatment in two of these strains, namely WT cells (BY4741, MUP1-pHluorin) and cells with the catalytically dead UL36-CD domain (MUP1-pHluorin-UL36-CD) (Figure 4C, AUC of Myr-treated vs. untreated cells, 95% CL: −2.710 to −1.018, p-value = 0.0036 for WT and CI: −3.450 to −1.346, p-value = 0.0032 for UL36-CD). These significance values are probably an underestimate because both strains tend to regrow or gasp after day 8 which is why data are not shown beyond this time point. In contrast, CLS is only slightly enhanced by Myr treatment in cells with a catalytically active UL36 domain (Figure 4C, AUC, 95% CI: −1.864 to −0.0284, p-value = 0.05). Statistical significance was also evaluated by using a t-test for the day 8 survival values for Myr treated cells: WT vs. UL36 (p-value = 0.024) and UL36-CD vs. UL36 (p-value = 0.032). Taken together, these data show that Myr-enhanced longevity depends on Mup1 ubiquitination.

Our results reveal that Myr-enhanced longevity is prevented by Mup1 deubiquitination, which is a gain-of-function condition since it promotes PM stability. To determine how MUP1 loss of function affects drug-enhanced longevity, we compared the CLS of mup1Δ cells to wild-type cells and cells with the mup1Δ allele replaced with wild-type MUP1. Myr treatment significantly enhanced CLS in each of the three strains (Figure 4D, for WT BY4741, mup1Δ and MUP1 the AUC for Myr-treated vs. untreated cells, p-values ≤ 0.009, 0.026 and 0.0024, respectively). The viability curves for untreated mup1Δ and MUP1 cells do not differ significantly (AUC, 95% CI: −1.682 to 0.860, p-value = 0.42). We conclude from these data that deleting MUP1 does not prevent Myr-enhanced longevity.

Given the increased UBI4 transcript level after 4–6 h of Myr treatment (Figure 4A), we hypothesized that Myr treatment might affect total ubiquitin and its cellular distribution. To test this, we affinity purified FLAG-ubiquitin (from yeast cells harboring N-terminal FLAG fusions at two endogenous ubiquitin-encoding loci (RPS31 and RPL40B)) from untreated or Myr-treated yeast cells. Immunoblot analysis revealed that Myr treatment for 2 h or 4 h did not significantly alter the amount of FLAG-ubiquitin recovered or the amount of K48-linked ubiquitin polymers recovered (Figure 4E). In contrast, we observed a significant increase in K63-linked ubiquitin polymers after 4 h of Myr treatment (Figure 4E, right lane in bottom immunoblot), suggesting that increased conjugation of ubiquitin in K63-linked polymers is part of the cellular response to depletion of sphingolipids.

To confirm this result, and to resolve other linkage types, we performed stable isotope labeling with amino acids in cell culture (SILAC) on yeast cells and subjected heavy-labelled and light-labelled cells to mock-treatment and Myr-treatment, respectively. At 2 hours and 4 hours of treatment, we collected cells, prepared lysates, and affinity purified FLAG-ubiquitin, followed by mixing, tryptic digestion, and processing of peptides for analysis by mass spectrometry. This SILAC-MS analysis resolved various peptides corresponding to the different linkage types of ubiquitin polymers, including K6-linked, K11-linked, K48-linked, and K63-linked polymers. Interestingly, we detected only modest changes in linkage types after 2 hours of Myr treatment, while 4 hours of Myr treatment resulted in increased formation of several polymer types, most notably for K63-linked ubiquitin polymers (Figure 4F). This finding is consistent with the immunoblot results shown in Figure 4E, confirming that yeast cells increase the formation of K63-linked ubiquitin polymers in response to Myr treatment. Taken together, these results suggest that yeast cells respond to Myr treatment by remodeling ubiquitin pools, partly by increasing production of ubiquitin via transcription of UBI4 and partly by deploying existing ubiquitin pools to promote increased formation of K63-linked ubiquitin conjugates.

Strains, culture conditions, lifespan assays and statistical significance

Strains (Table 3), culture conditions, lifespan assays and their statistical significance were similar to ones described previously [7]. Specifically, yeast cells, made prototrophic by transformation with pHLUM (carries HIS3, LEU2, URA3 and MET, Addgene, Watertown, MA, USA [37]) and stored at −80°C, were streaked onto minimal glucose plates and incubated at 30°C for 3 days. One to three colonies were inoculated into 5 ml of SD complete medium (SDC) [38] and incubated on a gyratory water bath shaker for 18–24 h at 30°C. Typical culture densities were 7–9 A600 nm/ml (1.5 × 10⁷ cells/A600 nm). Cells were diluted into 25 ml fresh SDC warmed to 30°C in a 125 ml long-neck flask, containing vehicle (EtOH) or Myr (details described in the next paragraph), to give a starting cell density of 0.15 A600 nm. Flasks were incubated at 30°C on a gyratory air-bath shaker (200 rpm). After 72 h of incubation cells were diluted (1:50 × 1:100) and from 10 to 40 μL was spread on a YPD plate [38] and incubated 48–72 h at 30°C followed by colony counting. CLS assays were performed at least twice using triplicate cultures. Concentrations of Myr used in CLS assays ranged from 475–525 ng/ml depending on the sensitive of strains compared to wild-type, prototrophic BY4741. Drug-sensitivity was measured by culture density (A600 nm) after 24 h of growth using cultures started at 0.15 A600 nm units of cells/ml. All BY4741 yeast strains used for lifespan assays and for RNA extraction were made prototrophic by transformation with the pHLUM plasmid. The DUB (UL36) fusion yeast strains used in this study were generated by homologous recombination in BY4741 and SEY6210 background strains using the reagents and strategy previously described (Hepowit et al. 2022). The MUP1 knock-in strain (NHY945) was generated in BY4741 by swapping the endogenous MUP1 coding region with NATMX (NHY930), shuffling the chromosomally integrated NATMX with URA3 (NHY938.1), and snipping out URA3 by homologous reintegration of a PCR-amplified MUP1 coupled with counter selection on 5-fluoroorotic acid (5-FOA) synthetic media plate. A two-tailed Student’s t-test and Area Under the Curve (AUC) (Sigma Plot) were used to evaluate statistical significance of CLS assays done with at least 3 biological replicates. Results were verified by one or more repeat experiments.

Budding yeast mRNA-enriched profiles were generated from 42 individual samples of prototrophic BY4741 cells, including 3 replicates at 7 time points and 2 treatments (control and drug-treated). Culture conditions were similar to ones described previously (Hepowit et al. 2021, Aging) except for the following modifications. Prototrophic BY4741 yeast cells were grown in 200 ml of SDC culture medium in a 1 L flask with the medium heated to 30°C before addition of EtOH (final concentration of 0.3% for control samples). For drug-treated samples, myriocin was added after addition of EtOH to give final concentrations of 0.3% EtOH and 700 ng/ml myriocin. Lastly, cells from a saturated overnight culture were diluted into the 200 ml cultures to give an initial A600 nm units/ml of 0.15. Flasks were incubated at 30°C and 200 rpm on a rotary shaker. Control and myriocin-treated cells were harvested at time 0, 1, 2, 3, 4, 5 and 6 h (all time points for each replication were from the same culture flask).

RNA extraction

Culture conditions, medium and prototrophic BY4741 cells were the same as for lifespan and amino acid pool assays [7]. RNA was extracted from 5 A600 nm unit/ml of yeast cells harvested by filtration on a membrane filter at each time point (1–6 h, Figure 1A) [7]. Filtered cells were washed once with 5 ml of ice-cold nanopure water and the filter was quickly transferred to a chilled 1.5 ml microfuge tube containing 0.5 ml cold nanopure water. Tubes were vortexed 10 sec followed by centrifugation for 15 sec. Supernatant fluid (450 μl) was transferred to a new tube and frozen in a dry-ice EtOH bath followed by storage at −80°C. Cold acid-washed glass beads (300 μl, 0.5 mm dia.) were added to a frozen cell pellet followed by addition of 300 μl of RLT buffer (RNAeasy mini kit, Qiagen, Germantown, MD, USA). Tubes were vortexed 5 min at room temperature and placed on ice for 1 min. After 4 cycles 300 μl of ice-cold RLT buffer was added followed by mixing and centrifugation at 13,000 × g for 2 min at room temperature. Supernatant fluid (450 μl) was transferred to a new microfuge tube and then mixed with 1 ml of 70% (v/v) EtOH before transfer to a RNeasy spin column and processed according to the manufacturer’s instructions. Processed samples were frozen in a dry-ice EtOH bath and stored at −80°C.

RNA-seq analysis

RNA sequencing was performed on total RNA samples at the Roy J. Carver Biotechnology Center at the University of Illinois. Two different mixes of ERRC spike-in RNA controls were added to the samples; one mix for the control samples, and one mix for the drug-treated samples. Libraries were constructed with Illumina’s ‘TruSeq Stranded mRNAseq Sample Prep kit’. Each library was quantitated by qPCR and sequenced on one lane for 151 cycles from each end of the fragments on a HiSeq 4000 using a HiSeq 4000 sequencing kit version 1. Fastq files were generated and demultiplexed with the bcl2fastq v2.17.1.14 Conversion Software (Illumina). Each sample’s pair of fastq files were run through trimmomatic 0.36 to first remove any remaining standard Illumina PE v3 adapters, then trim bases from both ends with quality scores below 28, and finally to remove individual reads shorter than 30 bp and their paired read, regardless of length (parameters ILLUMINACLIP:/home/apps/trimmomatic/trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:15:10 TRAILING:28 LEADING:28 MINLEN:30). Paired reads per sample were pseudo-aligned to the Yeast R64 transcriptome and 92 ERRC spike-in control sequences using Salmon 0.8.2 with parameters -l A –numbootstraps = 30 --seqBias –gcBias.

Resulting FASTQ files mapped to the yeast genome (R64), resulted in count files and normalized using the transcripts per million (TPM) algorithm [39] using WebMev [40]. Resulting data were downloaded as flat files and loaded into Excel for further analysis. From a total of 6198 mapped genes, 5169 were uniquely annotated with gene symbols and had sufficient non-zero readings for further analysis. The filtered data were analyzed by two-way ANOVA for the main effects of drug and time, as well as for interaction and assigned a p-value. A gene found to be significant in any of these three categories was referred to as differentially expressed (DEG). Significant by the time term or both the drug and time term, data were further analyzed by post-hoc pairwise Fisher’s protected Least Significant Difference (pLSD), and log 2-fold change comparison to further isolate the effects of drug over time.

To identify DEGs that change up or down over time we built a ‘template assignment tool’ that took each DEG and attempted to assign it to each of 63 patterns. The pattern to which a DEG fit best (by Pearson’s correlation, had to be R > 0.85) became the pattern to which it was assigned. This is reported in Figure 1C (99.86% of all DEGs were assigned to a template). To determine which templates had significantly more DEGs than expected by chance, we ran a binomial test- this is the basis for some templates being bolded and having an asterisk in Figure 1C- they had a binomial p < 0.01. The 12 patterns of expression that had more DEGs than expected by chance assigned to them explained the majority (77%) of all DEGs. Data have been deposited in the GEO (GSE199904) [NCBI tracking system #22817261]. GO terms were analyzed by using DAVID software [41] or online WebGestalt.

Gene regulatory network using WGCNA

The WGCNA (v1.70-3) [28] was used to identify gene modules and build unsigned co-expression networks, which include negative and positive correlations. Briefly, WGCNA constructs a gene co-expression matrix, uses hierarchical clustering in combination with the Pearson correlation coefficient to cluster genes into groups of closely co-expressed genes termed modules, and then uses singular value decomposition (SVD) values as module eigengenes to determine the similarity between gene modules or to calculate association with sample traits (for example, incubation time or amino acid levels). The top 2,000 variable genes were used to identify gene modules and network construction. Soft power 8 was chosen by the WGCNA function pick SoftThreshold. Then, the function TOMsimilarityFromExpr was used to calculate the TOM similarity matrix via setting power = 8, networkType = "signed. The distance matrix was generated by subtracting the values from similarity adjacency matrix by one. The function flashClust (v.1.01) was used to cluster genes based on the distance matrix, and the function cutreeDynamic was utilized to identify gene modules by setting deepSplit = 3. Cytoscape (v.3.8.2) was applied for the network visualizations.

Analysis of ubiquitin and its linkage types

For immunoblotting assays 5 A600 nm unit of cells were precipitated in 10% trichloroacetic acid on ice for 30 min. The protein precipitate was washed twice with ice-chilled acetone, lyophilized by vacuum centrifugation, and resuspended in urea sample buffer (75 mM Tris-HCl [pH 6.8], 6 M urea, 1 mM EDTA, 3% SDS, 20% glycerol, bromophenol blue), heated at 65°C for 5 min, and vortexed for 5 min. Proteins were separated by SDS-PAGE and transferred onto Immobilon-P^SQ membrane (0.2 μm; Millipore). Immunoblotting was performed using the following primary antibodies: anti-ubiquitin (1:10,000; LifeSensors; MAb; clone VU-1), anti-K48 (1:10,000; Cell Signaling; RAb; clone D9D5), and anti-K63 (1:4000; EMD Millipore; RAb; clone apu3). Secondary antibodies were anti-mouse (IRDye 680RD-Goat anti-mouse) or anti-rabbit (IRDye 800CW-Goat anti-rabbit) obtained from LI-COR Biosciences). Blots were scanned using Odyssey CLx and signal fluorescence was visualized using Image Studio Lite (LI-COR Biosciences). Ubiquitin linkage types in JMY1312 yeast cells were determined and quantified by SILAC-based mass spectrometry as previously described except for a change in the cell lysis buffer: 50 mM Tris-HCl [pH7.5], 150 mM NaCl, sodium pyrophosphate, 20 mM β-glycerophosphate, 2 mM sodium orthovanadate, 1 mM phenylmethylsulfonyl fluoride, 0.2% NP-40, 10 mM iodoacetamide, 20 μM MG132, 1 mM 1,10-phenanthroline, 1× EDTA-free protease inhibitor cocktail (Roche), 1× PhosStop (Roche) [42].