Systems-level design principles of metabolic rewiring in an animal

Abstract

The regulation of metabolism is vital to any organism and can be achieved by transcriptionally activating or repressing metabolic genes^1-3. Although many examples of transcriptional metabolic rewiring have been reported⁴, a systems-level study of how metabolism is rewired in response to metabolic perturbations is lacking in any animal. Here we apply Worm Perturb-Seq (WPS)-a high-throughput method combining whole-animal RNA-interference and RNA-sequencing⁵-to around 900 metabolic genes in the nematode Caenorhabditis elegans. We derive a metabolic gene regulatory network (mGRN) in which 385 perturbations are connected to 9,414 genes by more than 110,000 interactions. The mGRN has a highly modular structure in which 22 perturbation clusters connect to 44 gene expression programs. The mGRN reveals different modes of transcriptional rewiring from simple reaction and pathway compensation to rerouting and more complex network coordination. Using metabolic network modelling, we identify a design principle of transcriptional rewiring that we name the compensation-repression (CR) model. The CR model explains most transcriptional responses in metabolic genes and reveals a high level of compensation and repression in five core metabolic functions related to energy and biomass. We provide preliminary evidence that the CR model may also explain transcriptional metabolic rewiring in human cells.

Extended Data Fig. 1:. WPS dataset.

a, Genes covered in the WPS experiments. b, Coverage of metabolic pathways and categories. Small pie-chart indicates all genes in iCEL1314 excluding six pseudogenes. c, Scatter plot showing the relation between the number of DEGs in a WPS experiment versus the level of knockdown of the RNAi target gene. Most WPS-targeted genes exhibited a solid reduction in mRNA levels and genes with low efficiency tended to be lowly expressed, which suggest that they may not be confidently quantified and/or may not be functional under the conditions tested. Color indicates level of expression in wild type (WT) animals at the adult stage. d,e, WPS responsiveness for nonredundant and redundant iCEL1314 genes. Nonredundant genes are defined as those that are not buffered by any isoenzymes as determined by Gene-Protein-Reaction associations in iCEL1314 (Supplementary Methods). Bar plot shows the proportion of genes that are responsive when perturbed. Box plot shows the number of DEGs in each perturbation, where the box bounds the IQR divided by the median, and whiskers extend to a maximum of 1.5 × IQR beyond the box. f, WPS responsiveness for iCEL1314 genes by WormPaths Level 1 metabolic pathway or category. g, WPS responsiveness for iCEL1314 genes by WormPaths Level 4 metabolic pathway or category. In both (f) and (g), the gene targeted by RNAi was excluded from counting DEG for each pathway. h, Animal body size in WPS perturbations measured one hour before RNA extraction. Each dot represents the average body size of all animals in three replicates for each perturbation. i, Showing the body size distribution of every animal of three example perturbations (red dots in (h)). n = 164 animals (vector RNAi), n = 85 animals (fum-1 RNAi), n = 59 animals (elo-5 RNAi), n = 110 animals (let-754 RNAi). Animals from the same biological replicate are labeled with dot colors. Box bounds the IQR divided by the median, and whiskers extend to a maximum of 1.5 × IQR beyond the box. A notch is added to the box approximating a 95% confidence interval for the median. a.u., arbitrary unit.

Extended Data Fig. 2:. mGRN characteristics.

a, Distribution of number of interactions versus the Log2(Fold Change, FC) in the mGRN. b, Plot depicting fraction of total interactions versus fraction of individual genes showing that 20% of genes that are differentially expressed contribute 60% of the total interactions in the mGRN (red dashed lines). c, Log-log plot showing the power-law distribution of k-in in the mGRN. d, Relationship between level of gene expression in wild type (WT) animals and whether the gene is decreased (down) or increased (up) in expression when called as DEG. e, Bar graph of FC consistency, which is defined by the fraction of DEG interactions that is increased or decreased, whichever is greater, among all interactions of a gene in the mGRN. The red dashed line indicates 70% of genes that showed FC consistency over 90%. f, Volcano plot of k-in versus median log2(FC) of each gene when identified as a DEG across perturbations. Hub genes affected by many perturbations are labeled and their association with WormCat categories is colored as indicated. g, WormCat level 1 enrichment of DEGs for each perturbation, with level 2 indicated in inset (adjusted P value (P_adj) < 0.01). Only perturbations with more than 10 DEGs (excluding the perturbed gene) were analyzed.

Extended Data Fig. 3:. WPS perturbation clusters.

Detailed view of Fig. 2c with gene names indicated.

Extended Data Fig. 4:. WPS perturbation-perturbation similarity and gene-gene correlation.

a, Association of cosine similarities of perturbations with indicated features of gene pairs. b, Cosine similarity for pairs of perturbations corresponding to genes encoding protein complex members, in flux coupled reactions, in the same metabolic pathways, and for all other pairs. Only pairs with high cosine similarity (> 0.2) are shown. When a gene pair appears in more than one of the three functional categories, we prioritized the labeling in the order of ‘complex’, ‘flux coupling’, and ‘WormPaths pathway’. c, Tissue expression annotation of the gene-gene correlation hierarchical tree using z-transformed tissue TPM data from L2 animals. d, Number of total mGRN interactions associated with the two gene-gene correlation clusters. e, Annotation of gene-gene correlation clusters and GEPs. Most GEPs were defined by a single, functionally annotated gene-gene correlation cluster, while a few included multiple clusters that contained functionally related genes and stemmed from the same branch of the tree. Technical details can be found in Supplementary Methods. f, Bubble plot showing the regulation of GEPs by nuclear hormone receptors (NHRs), SREBP (sbp-1), MAP kinase (mpk-1) and TOR (let-363) (Supplementary Table 3). The size of the bubble shows the enrichment P values of genes in a GEP among the DEGs in the WPS profile of a regulator. P values were calculated using a one-tailed hypergeometric test and adjusted for multiple comparisons using the Benjamini-Hochberg (BH) method. Color shows the average of corrected Wald statistic of the GEP genes in the corresponding WPS profile. Only GEPs that showed significant association (P_adj < 0.01) are visualized.

Extended Data Fig. 5:. Effect sizes of perturbation clusters and GEPs.

a, Relationship between number of associated GEPs (defined by median P_adj < 0.05) and number of associated DEGs for the 22 perturbation clusters. The number of associated DEGs for a cluster is visualized by the mean ± s.d. of DEG numbers for perturbations in the cluster. b, Relationship between the number of perturbations and perturbation clusters that are associated with the 44 GEPs. The associated perturbations and clusters are defined by P_adj and median P_adj less than 0.05, respectively (details in Supplementary Methods).

Extended Data Fig. 6:. mGRN decomposition.

a, Diagram illustrating how the plots in Fig. 3a and 3b were derived. b, Bubble plot showing the tissue enrichment of the 44 GEPs using single cell RNA-seq data from adult animals. The size of the bubble indicates the proportion of genes from each GEP enriched in the indicated tissue and the color indicates the significance of the enrichment (one-tailed hypergeometric test, multiple testing adjusted by BH method). Tissue-enriched genes were defined by thresholding the relative tissue expression (z-score > 2, Supplementary Methods). c, Same as (b) but for the 22 perturbation clusters. d, Network diagram showing potential tissue-tissue interactions revealed by the mGRN. Each edge represents an interaction between the tissue in which a perturbed gene is enriched (source node) and the tissue in which each of its DEG are enriched (target node). The width of the edge is proportional to the total number of interactions. Red indicates activation and blue indicates repression. Tissue enriched genes were defined by thresholding the z-score of tissue expression levels (z-score > 2, Supplementary Methods).

Extended Data Fig. 7:. Examples of GEP scatter plots.

a, Example of detailed view of one GEP where each dot is a perturbation, and where perturbations are separated on the x-axis by the number of associated DEGs and on the y-axis by the GEP Wald statistic. The GEP Wald statistic is a test statistic following standard normal distribution and describing the GEP changes in the perturbation (Extended Data Fig. 6a, Supplementary Methods). Larger dots indicate significant change (P_adj < 0.05, Supplementary Methods) of the GEP whose targeted gene name and cluster are labeled by corresponding color. b-e, Four additional examples of GEP scatter plots indicating the perturbations associated with indicated GEPs. All 44 GEP scatter plots are available on the WormFlux-WPS website (https://wormflux.umassmed.edu/WPS).

Extended Data Fig. 8:. Transcriptional rewiring of metabolism.

a, b, Enrichment P value (a) and fraction (b) of all metabolic and iCEL1314 genes in the DEGs of each of the 385 responsive WPS experiments. P values were derived from hypergeometric test and directly visualized without multiple testing adjustment. c, Metabolic rewiring of Met/SAM cycle genes by Met/SAM cycle perturbations. The top panel shows the pathway cartoon with the genes perturbed by WPS indicated with a black dot, and genes upregulated in red. The middle bar plot indicates perturbations on x-axis and on the y-axis the corrected Wald statistic of the red genes that are numbered according to the top pathway graph. The scatter plot at the bottom shows all perturbations and their associated number of DEGs versus the average Wald statistic of the red genes. d, WormPaths level 4 annotation of within-reaction and within-pathway compensation. e, Activation of genes in de novo purine synthesis pathway by perturbations in purine salvage and interconversion pathways. The scatter plot shows all perturbations and their associated number of DEGs versus the average Wald statistic of the red genes in the pathway cartoon. f, Transcriptional rerouting from branched chain amino acid breakdown to phenylalanine/tyrosine breakdown (left) where knockdown of genes indicated by a black dot activates genes indicated in red (middle). The scatter plot on the right shows the average Wald statistic of red genes in the network cartoon for all perturbations. g, Fractional contribution of tyrosine to TCA cycle metabolites (L_{TCA metabolite}/L_Tyr, L is the isotope labeling fraction) for indicated knockdowns, measured by [U-¹³C]Phe tracing. Each bar represents the mean (± s.d.) and each dot indicates a biological replicate (n= 3) collected at isotopic steady state (4 to 6 hours post tracing). P values were derived by comparing RNAi and vector control samples for each metabolite, calculated using two-tailed Student’s t-test and labeled on top of each bar.

Extended Data Fig. 9:. Examples of metabolic rewiring related to network coordination.

a, The scatter plot on the right shows all perturbations and their associated number of DEGs versus the average Wald statistic of the red genes in the network cartoon on the left. b, The subnetwork in GEP Metabolism: lipid synthesis (intestinal) is activated by perturbations in lipid metabolism (green), PPP (blue), energy metabolism (orange) and peroxisomal beta-oxidation (magenta). The subnetwork genes are shown as the rows of the heatmap. The scatter plot shows the average corrected Wald statistic of the subnetwork genes (rows of the heatmap) versus number of DEGs for each perturbation. c, Fractional contribution of glucose-6-phosphate (g6p) to serine in [U-¹³C]glucose tracing experiments of two perturbations in genes from the GEP ‘Metabolism: lipid synthesis’. Each bar represents the mean (± s.d.) and each dot indicates a biological replicate (n= 4). P values were calculated using two-tailed Student’s t-test. d, e, Similar graphs as (b) for two additional GEPs where (d) shows that the subnetwork in GEP Metabolism: energy production #1 is activated by perturbations in energy metabolism, e.g., TCA and ETC genes that belong to the subnetwork (red and orange). It is also repressed by perturbations in lipid synthesis (green) and oxidative PPP but activated by that in non-oxidative PPP (blue); and (e) shows that the subnetwork in GEP Metabolism: ketone, BCAA and SCFA oxidation is partially activated by perturbations in the GEP genes, however, broadly activated by perturbations in other genes that locate in the subnetwork but not in the GEP. f, Scatter plot showing the GEP expression changes upon perturbations of the corresponding GEP genes in each of the 14 metabolic GEPs. Each dot is a perturbation with the y-axis showing the corrected Wald statistic of the GEP its targeted gene resides in, and the x-axis indicating the number of DEGs obtained with this perturbation. Red dots indicate perturbations whose residing GEP was significantly (P_adj < 0.05, Supplementary Methods) activated when perturbed. The grey line and area indicate the background variation, which is defined as median plus and minus one Median Absolute Deviation (MAD) of the Wald statistics of GEPs in which the targeted gene does not reside.

Extended Data Fig. 10:. Analysis of C. elegans core metabolic functions using FBA

a, Manually identified interactions among three core metabolic functions (bottom left) based on perturbation clusters and GEPs. The size of the bubble indicates the significance of the change (P_adj) in a GEP at the perturbation cluster level and the color indicates corresponding cluster-level corrected Wald statistic for a GEP (statistical details can be found in Supplementary Methods). b, Schematic outlining the FBA simulation that assigns core metabolic functions to iCEL1314 genes. The FBA simulations produced a core function score that captures the effect of the removal of each gene (and associated reaction) on the maximal flux capacity and flux efficiency for each of the five core metabolic functions (Supplementary Table 9). Thresholding this score associates metabolic genes with their core metabolic functions. To facilitate interpretation we also refined some associations to reduce unspecific associations resulted from dominating effects. For instance, many energy production genes such as ETC genes scored high on all five core metabolic functions because energy production underlies all these functions. We assigned genes only to energy production if their perturbations influenced ATP production more than other processes. Similarly, we assigned genes primarily influencing protein synthesis to this function only, avoiding their assignment to ECM synthesis. Details are explained in Supplementary Methods. c, DEG rewiring patterns for all perturbations whose targeted genes are associated with any core metabolic function. The rows are ordered based on the hierarchical clustering of perturbations using the proportions of each core function in the upregulated DEGs (i.e., the height of each colored bars in the ‘Up DEG’ column). Cosine distance and centroid linkage were used for this clustering. The core metabolic function associations for each perturbation are shown on the right with the same color code. The colors of DEGs are determined based on the same majority-vote strategy as in Fig. 4c (Supplementary Methods).

Extended Data Fig. 11:. Analysis of C. elegans core metabolic functions and CR-model using FBA.

a, Visualizing DEGs for each RNAi perturbation annotated to the five core metabolic functions by FBA, only using the genes and perturbations that associate with unique core metabolic functions. b, Bar plot showing the fraction of iCEL1314 DEGs explained by the CR model and FBA in all 223 analyzed perturbations. c, d, e, Randomization (10,000) of gene-core metabolic function association showing that the observed fraction (red line) of condition-wise average (c), activated (d), and repressed (e) DEGs explained by the CR model is statistically significant. The empirical P values are shown on the plot. The condition-wise average is the average of explained DEG proportions of each individual condition. The core function score threshold used is 0.001. f, Network diagram showing tissue-tissue interactions for DEGs explained by the CR model. Each edge represents an interaction between a tissue-enriched perturbed gene (source node) and a tissue-enriched DEG affected by this perturbation (target node). The width of the edge is proportional to the total number of interactions between two tissues that are explained by the CR model. Red indicates activation and blue indicates repression. Tissue-enriched genes were defined by thresholding the z-score of tissue expression levels in single cell data from adult animals (z-score > 2, Supplementary Methods).

Extended Data Fig. 12:. Analysis of human core metabolic functions and CR-model using single-cell Perturb-seq data.

a, Euler diagram showing the number of genes assigned to each core metabolic function in the Human 1 metabolic network model. b, Distribution of normalized gene expression levels for DEGs of 1,769 perturbations in Human 1 genes covered by Perturb-seq. The normalized gene expression level is a z-score value representing the number of standard deviations from the mean in control cells, with positive (negative) values indicating up(down)-regulated genes. DEGs were defined as Anderson-Darling (AD) test P_adj < 0.01 (Supplementary Methods). c,d, Zoom-in version of Fig. 5a with gene names and numbers labeled. e, Bar plot showing the fraction of Human 1 DEGs explained by the CR model using FBA for all 190 perturbations analyzed. f, g, h, Randomization (10,000) of gene-core metabolic function association showing that the observed fraction (red line) of condition-wise average (c), activated (d), and repressed (e) DEGs explained by the CR model are statistically significant. The empirical P values are shown on the plot. The condition-wise average is the average of explained DEG proportions of each individual condition. The core function score threshold here is 0.001 and DEG P_adj threshold is 0.01. i, Fraction of total DEGs explained and the corresponding empirical P values as a function of the metabolic core function score threshold in the human CR model analysis.

Fig. 1:. A genome-scale, whole-animal mGRN for C. elegans.

a, Modes of transcriptional rewiring of metabolism, including pathway-level rerouting (top right) and compensation (middle right). b, WPS coverage for different types of genes and reactions. Nonredundant genes are defined as those that are not buffered by any predicted isoenzymes in iCEL1314. Perturbable reactions are defined as reactions associated with nonredundant genes. Red dashed line indicates 95%. c, Out-degree (k-out) distribution of perturbed genes. Out-degree is defined as the number of DEGs per knockdown. A cutoff of five DEGs was used as a threshold for responsiveness based on false discovery rate (FDR) benchmarking. d, Single and multiple RNAi knockdown followed by RNA-seq for two EC families. The top panel shows number of DEGs, the bottom panel shows change in mRNA level of the knocked down gene(s). e, The mGRN contains 110,930 perturbation-DEG interactions (FDR < 0.1, fold change (FC) > 1.5) among 385 perturbations and 9,414 genes. f, Fraction of upregulated genes versus the number of WPS perturbations. g, In-degree (k-in) distribution of the 9,414 genes in mGRN. In-degree is defined as the number of perturbations affecting each gene. For (f) and (g), red arrows indicate upregulated DEGs and blue arrows indicate downregulated DEGs.

Fig. 2:. Modular organization of the C. elegans mGRN.

a, Heatmap of 385 responsive perturbations versus 9,414 genes that are differentially expressed in at least one perturbation (FDR < 0.1, FC > 1.5). Values indicate the corrected Wald statistic that is a test statistic for differential expression following standard normal distribution (Supplementary Methods). b, Perturbation-perturbation similarity matrix (measured by cosine similarity as indicated by color scale. See Supplementary Methods for details). c, 2-dimensional visualization of perturbation clusters by densMAP that preserves the local density information in the reduced dimension. We rigorously checked the quality of the clusters to avoid potential problems introduced by the embedding (Supplementary Fig. 2b–f, Supplementary Note 2). Unclustered genes are colored grey. d, Gene-gene correlation matrix of genes identified as differentially expressed in any perturbation. The matrix shows Pearson correlation coefficients between genes calculated using their corrected Wald statistics in the 385 responsive perturbations. The matrix was organized by a hierarchical tree, which is identical to that in (e) but omitted here to save space. e, Tissue expression of all genes in the matrix of (d). Tissue expression levels were derived from a previous single-cell RNA-seq dataset of adult stage C. elegans. f, Bubble plot of the 22 metabolic perturbation clusters versus the 44 GEPs. The size of the bubble indicates the significance of the change (adjusted P value (P_adj)) in a GEP at the perturbation cluster level and the color indicates corresponding cluster-level corrected Wald statistic for a GEP (details can be found in Supplementary Methods). The average tissue expression levels (from same data in (e)) of genes in corresponding GEPs or perturbation clusters are shown on the right and top, respectively. mal, male; int, intestine; neu, neuron; hyp, hypodermis; mus, muscle; ger, germline; gon, gonad; pha, pharynx.

Fig. 3:. Modes of transcriptional rewiring of C. elegans metabolism.

a, Cartoon illustrating compensation within a reaction or pathway. b, Number of genes whose perturbation triggered within reaction, pathway, or both types of transcriptional compensation. c, Breakdown of numbers in (b) for top WormPaths pathways. d, e, Randomization of mGRN by edge swapping (Supplementary Methods, 10,000 times randomizations) shows that the observed frequency (red lines) of reaction-level (d) and pathway-level (e) compensation are statistically significant. The empirical P value is indicated. f, Cartoon illustrating metabolic rerouting. g, Example of metabolic rerouting in purine metabolism indicating the pathways (left) with the genes perturbed by WPS indicated with a black dot, and genes upregulated in red. Bar graphs indicate perturbations (bottom) and corrected Wald statistic of red DEGs numbered according to the left pathway graph. h, Isotope tracing using heavy glycine (heavy atom in red) in vector control and hprt-1 knockdown animals. Each bar represents the mean (± s.d.) and each dot indicates a biological replicate (n= 4). P values were calculated using two-tailed Student’s t-test. i, Cartoon illustrating network coordination rewiring. j, Example of compensation within a subnetwork. Genes with black dots indicate WPS perturbations and those in red are iCEL1314 genes in the GEP (left). The heatmap shows changes in the subnetwork genes, including red genes and acly-1 in the network cartoon, on y-axis versus perturbations on the x-axis. The Wald statistics of knocked down genes were masked to zero. k, Compensation rate of 14 metabolic GEPs, which is defined as the frequency of a gene residing in a GEP it activates (P_adj < 0.05, Supplementary methods) when perturbed. The frequency was calculated with respect to responsive perturbations of genes from each GEP (the denominator of the fraction shown on top of each bar).

Fig. 4:. A compensation/repression (CR) model of transcriptional rewiring of C. elegans metabolism.

a, Assigning genes to five core metabolic functions by FBA. The cartoon on the left shows the biomass fraction (w/w) of different biomass components based on iCEL1314. For algorithmic details, see Extended Data Fig. 10b. b, Example showing the core metabolic functions of DEGs (in columns) elicited by gspd-1 knockdown. The color of the most prevalent core metabolic function was used when a gene has multiple associations. c, Visualization of WPS-derived DEGs for 160 perturbations of genes associated with a single core metabolic function. d, Cartoon illustrating interactions among the five core metabolic functions, statistically supported by randomization of perturbation-DEG associations (left), or core function assignments (right). Each edge indicates the up (red) or down (blue) regulation of target node genes in response to perturbations of source node genes. Nodes are colored following (c). Only edges with empirical P value < 0.1 are shown. DEG contribution indicates the average fraction of DEGs contributed by the edge. e, CR model concept. f, An example for gspd-1 RNAi illustrating the validation of CR model. Note that the numbers are slightly different from a manual calculation using numbers in (b), as (b) shows the refined core metabolic function associations for better data visualization (see Supplementary Note 8). g, Statistical significance of CR model. Histogram shows the proportion of total DEGs explained by the CR model in 10,000 randomizations of the gene-core metabolic function associations and the red vertical line showing the fraction observed in the real data. h, Fraction of total DEGs explained and the corresponding empirical P values as a function of the metabolic core function score threshold. Grey dashed line indicates the cutoff used in (g).

Fig. 5:. Testing CR model in human cells.

a, Visualization of metabolic rewiring DEGs for 148 perturbations of genes associated with a unique core metabolic function by FBA. DEGs were defined by Anderson-Darling (AD) test of the perturbation and control cell populations (P_adj < 0.01, Supplementary Methods). The color code of core functions is the same as for C. elegans in Fig. 4. b, Network cartoon illustrating the interactions among the five core metabolic functions that are statistically supported by randomization of core function assignments. The network is visualized as in Fig. 4d. Only edges with empirical P < 0.1 are shown, with a single exception for the compensation edge of ECM perturbation (yellow node), whose empirical P value is 0.11. c, Proportion of total DEGs explained by the CR model for 190 perturbations of genes with FBA-defined core metabolic functions. The bar graph indicates results from 10,000 randomizations of the gene-core metabolic function associations and the red vertical line showing the fraction observed in the real data. d, Fraction of total DEGs explained and the corresponding empirical P values as a function of the DEG P_adj threshold.