Reconstructing a metazoan genetic pathway with transcriptome-wide epistasis measurements

Abstract

RNA-sequencing (RNA-seq) is commonly used to identify genetic modules that respond to perturbations. In single cells, transcriptomes have been used as phenotypes, but this concept has not been applied to whole-organism RNA-seq. Also, quantifying and interpreting epistatic effects using expression profiles remains a challenge. We developed a single coefficient to quantify transcriptome-wide epistasis that reflects the underlying interactions and which can be interpreted intuitively. To demonstrate our approach, we sequenced four single and two double mutants of Caenorhabditis elegans From these mutants, we reconstructed the known hypoxia pathway. In addition, we uncovered a class of 56 genes with HIF-1-dependent expression that have opposite changes in expression in mutants of two genes that cooperate to negatively regulate HIF-1 abundance; however, the double mutant of these genes exhibits suppression epistasis. This class violates the classical model of HIF-1 regulation but can be explained by postulating a role of hydroxylated HIF-1 in transcriptional control.

Keywords: epistasis; gene expression; genetics; hypoxia; transcriptomics.

Fig. 1.

Genetic and biochemical representation of the hypoxia pathway in C. elegans. Red arrows are arrows that lead to inhibition of HIF-1, and blue arrows are arrows that increase HIF-1 activity or are the result of HIF-1 activity. EGL-9 is known to exert VHL-1–dependent and –independent repression on HIF-1 as shown in the genetic diagram. The VHL-1–independent repression of HIF-1 by EGL-9 is denoted by a dashed line and is not dependent on the hydroxylating activity of EGL-9. RHY-1 inhibits CYSL-1, which in turn inhibits EGL-9, but this interaction was abbreviated in the genetic diagram for clarity.

Fig. 2.

Analysis workflow. After sequencing, reads are quantified using Kallisto. Bars show estimated counts for each isoform. Differential expression is calculated using Sleuth, which outputs one $\beta$ coefficient per isoform per genotype. $\beta$ coefficients are analogous to the natural logarithm of the fold-change relative to a wild-type control. Downstream analyses are performed with $\beta$ coefficients that are statistically significantly different from 0. $q$ values less than 0.1 are considered statistically different from 0.

Fig. 3.

PCA of various C. elegans mutants. Genotypes that have a constitutive hypoxia response [i.e., egl-9(lf)] cluster far from genotypes that do not have a hypoxic response [i.e., hif-1(lf)] along the first principal component. The second principal component separates genotypes that do not participate in the hypoxic response pathway.

Fig. 4.

Interacting genes have correlated transcriptional signatures. The rank order of transcripts contained in the shared transcriptional phenotype is plotted for each pairwise combination of genotypes. Correlations between in-pathway genotypes are strong, whereas comparisons with a fog-2(lf) genotype are dominated by noise. Comparisons between some genotypes show populations of transcripts that are anticorrelated, possibly as a result of feedback loops. Plots are color-coded by row. Comparisons with genotypes with a constitutive hypoxia response are in blue, comparisons with genotypes negative for hif-1(lf) are in black, and comparisons involving fog-2(lf) are in red. The x and y axes show the rank of each transcript within each genotype.

Fig. 5.

Quantification of epistasis transcriptome-wide. (A) Schematic diagram of an epistasis plot. The x axis on an epistasis plot is the expected coefficient for a double mutant under an log-additive model (null model). The y axis plots deviations from this model. Double mutants that deviate in a systematic manner from the null model exhibit transcriptome-wide epistasis ($s$). To measure $s$, we find the line of best fit and determine its slope. Genes that act log-additively on a phenotype (Ph) will have $s=0$ (null hypothesis, orange line), whereas genes that act along an unbranched pathway will have $s=-1/2$ (blue line). Strong repression is reflected by $s=-1$ (red line), whereas $s>0$ correspond to synthetic interactions (purple line). (B) Epistasis plot showing that the egl-9(lf); vhl-1(lf) transcriptome deviates significantly from a null additive. Points are colored qualitatively according to density (purple, low; yellow, high) and size is inversely proportional to the SE of the y axis. The green line is the line of best fit from an orthogonal distance regression. (C) Comparison of simulated epistatic coefficients against the observed coefficient. Green curve shows the bootstrapped observed transcriptome-wide epistasis coefficient for egl-9 and vhl-1. Dashed green line shows the mean value of the data. Simulations use only the single mutant data to idealize what expression of the double mutant should look like. $a>b$ means that the phenotype of $a$ is observed in a double mutant $a^{-}{b}^{-}$.

Fig. 6.

Transcriptomes can be used to order genes in a pathway under certain assumptions. Arrows in the diagrams above are intended to show the direction of flow and do not indicate valence. (A) A linear pathway in which rhy-1 is the only gene controlling egl-9, which in turn controls hif-1, does not contain information to infer the order between genes. (B) If rhy-1 and egl-9 have transcriptomic effects that are separable from hif-1, then the rhy-1 transcriptome should contain contributions from egl-9, hif-1 and egl-9– and hif-1–independent pathways. This pathway contains enough information to infer order. (C) If a pathway is branched both upstream and downstream, transcriptomes will show even faster decorrelation. Nodes that are separated by many edges may begin to behave almost independently of each other with marginal transcriptomic overlap or correlation. (D) The hypoxia pathway can be ordered. We hypothesize the rapid decay in correlation is due to a mixture of upstream and downstream branching that happens along this pathway. Bars show the SE of the weighted coefficient from the Monte Carlo Markov Chain computations.

Fig. 7.

Fifty-six hif-1–dependent genes show nonclassical antagonistic effects of vhl-1 and egl-9. (A) A total of 56 genes in C. elegans exhibit nonclassical epistasis in the hypoxia pathway, characterized by opposite effects on gene expression, relative to the wild type, of the vhl-1(lf) compared with egl-9(lf) [or rhy-1(lf)] mutants. Shown are a random selection of 15 out of 56 genes for illustrative purposes. (B) Genes that behave noncanonically have a consistent pattern. vhl-1(lf) mutants have an opposite effect to egl-9(lf), but egl-9 remains epistatic to vhl-1 and loss-of-function mutations in hif-1 suppress the egl-9(lf) phenotype. Asterisks show $\beta$ values significantly different from 0 relative to wild type ($q<{10}^{-1}$).

Fig. 8.

A hypothetical model showing a mechanism where HIF-1-OH antagonizes HIF-1 in normoxia. (A) Diagram showing that RHY-1 activates EGL-9. EGL-9 hydroxylates HIF-1 in an oxygen-dependent manner. HIF-1 is rapidly hydroxylated, and the product, HIF-1-OH, is rapidly degraded in a VHL-1–dependent fashion. EGL-9 can also inhibit HIF-1 in an oxygen-independent fashion. In our model, HIF-1 and HIF-1-OH have opposing effects on transcription. The width of the arrows represents rates in normoxic conditions. (B) Table showing the effects of loss-of-function mutations on HIF-1 and HIF-1-OH activity, showing how this can potentially explain the ftn-1 expression levels in each case. S.S., steady state.