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 Day 3 adults. F1 progeny are either subject to a constant temperature or an eight-hour shift to 25°C. Beginning in midlarval stages all F1 worms experience the same environment. For each F1 worm, time to egg-laying onset is scored. 24 hours after egg-laying onset, each F1 worm is collected for single-individual mRNA-seq. F2 progeny laid on plates in the first 24 hours 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 red, non-significant genes are in black. Significance determined by a Bonferroni-corrected p-value cutoff of 5.67 x 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.

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 β1 is the coefficient of the linear mixed model. β1 represents the total effect of gene expression on early brood. In Model 2, β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, β1 is the total effect of grsp-4 on egg-laying onset, while in Model 2, β2 represents the effect of grsp-4 that is independent of historical environment.

C. The 448 genes with for which expression variation is significantly associated with early brood, ordered by the magnitude of their β1 effect sizes. The absolute values of β1 and β2 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 β1 effect sizes. The absolute values of β1 and β2 are plotted upward and downward, respectively, for each gene.

C-D. Legend indicates the -log10(p-value) of β1 and β2 for each gene.

E. 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 β2 but parental age does not have a significant effect on expression (‘noise’ genes), and 97 genes that do not have significant values of β2 but parental age does have a significant effect on expression (parental age genes).

Multi-gene models are highly predictive of early brood and egg-laying onset at an individual level.

A. Total R2 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 R2 of top genes in a multiple regression. Gene sets explain early brood (red), ELO (blue), and randomized data (dotted lines and gray shading).

A-B. Dotted lines represent average of randomized data iterations and gray shading indicates standard deviation.

C. Total R2 of top 10 genes in a multiple regression identified in a train set and total R2 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 are shown for early brood.

D-E. Machine learning model elastic net regression and leave-one-out cross validation used to identify a predicted trait for each worm given transcriptomic profile and compared to true trait data. Linear fit used to determine R2. 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.

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

A. 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.

A-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 x 10-6) and a significant main effect of tissue (p < 2 x 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 x 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, 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

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

Variation in 11 genes is significantly associated with egg-laying onset.

Each of 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 5.67 x 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.

Single-individual mRNA-seq reveals most and least variable genes after controlling for environmental perturbations and biological replicate.

A. Log10-transformed unexplained variance in mRNA-seq after controlling for environmental factors and biological replicate is plotted against log10-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 x 10-16