Gene expression variation across genetically identical individuals predicts reproductive traits

Abstract

In recent decades, genome-wide association studies (GWAS) have been the major approach to understand the biological basis of individual differences in traits and diseases. However, GWAS approaches have limited predictive power to explain individual differences, particularly for complex traits and diseases in which environmental factors play a substantial role in their etiology. Indeed, individual differences persist even in genetically identical individuals, although fully separating genetic and environmental causation is difficult in most organisms. To understand the basis of individual differences in the absence of genetic differences, we measured two quantitative reproductive traits in 180 genetically identical young adult Caenorhabditis elegans roundworms in a shared environment and performed single-individual transcriptomics on each worm. We identified hundreds of genes for which expression variation was strongly associated with reproductive traits, some of which depended on individuals' historical environments and some of which was random. Multiple small sets of genes together were highly predictive of reproductive traits, explaining on average over half and over a quarter of variation in the two traits. We manipulated mRNA levels of predictive genes to identify a set of causal genes, demonstrating the utility of this approach for both prediction and understanding underlying biology. Finally, we found that the chromatin environment of predictive genes was enriched for H3K27 trimethylation, suggesting that gene expression variation may be driven in part by chromatin structure. Together, this work shows that individual, non-genetic differences in gene expression are both highly predictive and causal in shaping reproductive traits.

Keywords: C. elegans; environmental variation; epigenetics; gene expression; genetics; genomics; reproduction.

Figure 1.. Differences in early brood across isogenic individuals are associated with expression levels of hundreds of genes.

(A) Experimental design for single-individual mRNA-seq. P0 worms are allowed to lay progeny as day 1 or 3 adults. F1 progeny are either subject to a constant temperature or an 8-hr shift to 25°C. Beginning in mid-larval stages, all F1 worms experience the same environment. For each F1 worm, time to egg-laying onset is scored. 24 hr after egg-laying onset, each F1 worm is collected for single-individual mRNA-seq. F2 progeny laid on plates in the first 24 hr of egg-laying are counted to determine the early brood of each F1 worm. (B) Each of 8824 expressed genes was analyzed in a mixed model to assess the strength of its association with early brood. Significant genes are in pink (positively associated) or purple (negatively associated), non-significant genes are in black. Significance determined by a Bonferroni-corrected p-value cutoff of 0.05 (nominal p-value of 5.67 × 10^–6). Phenotypic and expression data for each worm is shown for the top two genes, col-20 and phb-1, with linear model fits in black and 95% confidence intervals in gray. (C) Heatmap of 448 brood-associated genes for all 180 individuals. Worms are sorted from left to right from lowest to highest early brood. Early brood data and environmental information is shown in the bars at the top. Genes are sorted into two clusters based on whether the association was positive or negative.

Figure 1—figure supplement 1.. Isogenic worms at the same developmental stage exhibit variability in egg-laying onset and early brood that partially depends on previous environmental conditions.

(A) Univariate histograms of early brood and egg-laying onset trait values for all worms. (B) Bivariate scatterplots of early brood and egg-laying onset trait for all worms. The same data is plotted twice but color-coded according to two previous environmental perturbations (age of parent and early-life temperature). Center lines for box plots indicate median; box limits indicate upper and lower quartiles. Significance determined using a linear mixed model with environmental perturbation as a fixed effect and biological replicate as a random effect, ***p < 0.001, *p < 0.05.

Figure 1—figure supplement 2.. Variation in 11 genes is significantly associated with egg-laying onset.

Each of the 8824 expressed genes was analyzed in a mixed model to assess the strength of its association with egg-laying onset. Significant genes are in red, non-significant genes are in black. Significance determined by a Bonferroni-corrected p-value cutoff of 0.05 (nominal p-value of 5.67 × 10^–6). Phenotypic and expression data for each worm is shown for two exemplar genes, anmt-3 and rnp-8. Black lines from linear model fit and gray shading indicates 95% confidence intervals.

Figure 2.. Associations between gene expression and reproductive traits are partially driven by both noise and historical environment.

(A) Two example path analyses showing the relationship between gene expression and early brood. In Model 1, gene expression of a particular gene, puf-5, drives differences in early brood, and β₁ is the coefficient of the linear mixed model. β₁ represents the total effect of gene expression on early brood. In Model 2, β₂ represents the effect of puf-5 on early brood that is independent of the historical environment. (B) Two example path analyses showing the relationship between gene expression and egg-laying onset, analogous to A. In Model 1, β₁ is the total effect of grsp-4 on egg-laying onset, while in Model 2, β₂ represents the effect of grsp-4 that is independent of historical environment. (C) The 448 genes for which expression variation is significantly associated with early brood, ordered by the magnitude of their β₁ effect sizes. The absolute values of β₁ and β₂ are plotted upward and downward, respectively, for each gene. (D) The 11 genes for which expression variation is significantly associated with egg-laying onset across individuals, ordered by the magnitude of their β₁ effect sizes. The absolute values of β₁ and β₂ are plotted upward and downward, respectively, for each gene. For C and D, legend indicates the −log₁₀(p-value) of β₁ and β₂ for each gene. (E) Gene Ontology (GO) Terms for the 448 genes associated with early brood with a q-value <0.01 are shown. GO Terms for the two subsets of the 448 genes with overlapping GO Terms are also shown: 114 genes that have significant values of β₂ but parental age does not have a significant effect on expression (‘noise’ genes), and 97 genes that do not have significant values of β₂ but parental age does have a significant effect on expression (parental age genes).

Figure 3.. Multi-gene models are highly predictive of early brood and egg-laying onset (ELO) at an individual level.

(A) Total R² of gene expression principal components (PCs) in a multiple regression. PCs explain early brood (red), ELO (blue), and randomized data (dotted lines and gray shading). (B) Total R² of top genes in a multiple regression. Gene sets explain early brood (red), ELO (blue), and randomized data (dotted lines and gray shading). For A and B, dotted lines represent average of randomized data iterations and gray shading indicates standard deviation. (C) Total R² of top 10 genes in a multiple regression identified in a train set and total R² of the same 10 genes used in a test set. Train and test sets randomly selected 500 times and top 10 genes identified in each iteration. Box plot center lines indicate median and box limits indicate upper and lower quartiles. An example of one iteration of train and test sets and corresponding model fits is shown for early brood. For D and E, machine learning model elastic net regression and leave-one-out cross-validation were used to identify a predicted trait for each worm given transcriptomic profile and compared to true trait data. Linear fit used to determine R². Lines indicate linear fit, and gray indicates 95% confidence intervals. (D) Early brood model. (E) ELO model. (F) Proportion of 500 iterations from C in which a given gene is selected as one of the top 10 predictive genes for early brood, ordered by those selected from most to least often, and an inset showing the genes that are most frequently selected as predictive.

Figure 4.. Genes predictive of reproductive traits causally affect early brood and are enriched for H3K27me3.

(A) RNA interference (RNAi) knockdown of selected predictive genes compared to empty vector (ev) and effects on early brood. (B) Dose response of puf-5 and puf-7 RNAi together with empty vector and effects on early brood. For A and B, linear mixed model with RNAi as fixed effect and biological replicate as a random effect. Center lines, median; box limits, upper and lower quartiles, *p < 0.05, ***p < 0.001. (C) Effect sizes on brood for genes categorized by tissue (soma or germline) and chromatin domain (active or regulated). Two-way ANOVA showed significant interaction between tissue and chromatin domain (p = 1.58 × 10^–6) and a significant main effect of tissue (p < 2 × 10^–16). Post hoc Tukey tests corrected for multiple testing showed a significant difference between chromatin domains within somatic tissues (adjusted p = 0) but not germline tissues (adjusted p = 0.98). Within somatic genes, active somatic genes effects on average do not differ from 0 (p = 0.3, one-sample t-test with μ = 0), while regulated somatic genes on average have a significant negative association with brood (p < 2.2 × 10^–16, one-sample t-test with μ = 0), ***p < 0.001. (D) Histogram of all iterations identifying predictive gene sets and the proportion of each set in a regulated chromatin domain marked with H3K27me3 compared to randomly selected sets of genes from the same background set. Early brood iterations are shown in pink, egg-laying onset (ELO) in blue, and the same gray control is shown for both histograms. The vertical lines are color-coded to indicate the median proportion of genes in regulated domains in the corresponding condition. Kolmogorov–Smirnov test used to determine statistical significance, ***p < 0.001.

Figure 4—figure supplement 1.. Single-individual mRNA-seq reveals most and least variable genes after controlling for environmental perturbations and biological replicate.

(A) Log₁₀-transformed unexplained variance in mRNA-seq after controlling for environmental factors and biological replicate is plotted against log₁₀-transformed average CPM across all worms and a loess fit line is shown in black with a 95% confidence interval in gray. Each point is a different gene and genes are color-coded according to their distance from the fit line. Genes in purple are most variable given their level of expression, while genes in pink are the least variable. Exemplar genes with comparable expression levels are shown on the right. Center lines for box plots indicate median; box limits indicate upper and lower quartiles. (B) Variance Z-scores for genes with regulated domains (H3K27me3) compared to variance Z-scores for genes with active domains (H3K36me3). Wilcoxon test, ***p < 2 × 10^–16.