We raised set-2/+ heterozygous animals at 25°C (P0 generation) to obtain set-2 fertile F1 homozygotes that we then grew for one additional generation (F2). F2 animals were allowed to self-fertilize and produce progeny for two additional generations, and F4 progeny were scored and separated in two distinct sets: F4 fertile and F4 sterile. Note that although within each generation the fertility of individual animals is highly variable (Supplementary Figure S1A), by the F6-F8 generation all animals become sterile. We carried out germline dissections followed by RNA sequencing (RNA-seq) on three separate lineages started from set-2/+ single mothers (Figure 1A).

We first used the RNA-seq data sets to identify changes in gene expression across generations, considering lineages as independent replicates and comparing F2, F4 fertile, and F4 sterile animals to the parental samples (P0). This confirmed that loss of germline identity in set-2 mutants is associated with transcriptional deregulation (Robert et al., 2014), and showed that the number of misregulated genes progressively increased between P0, F2, F4 fertile, and F4 sterile animals (Figure 1B). In all generations, we observed significantly more upregulated than downregulated genes, as previously reported for healthy germlines from set-2 animals raised at 20°C (Robert et al., 2014). Furthermore, the total number of misregulated genes was approximately three times higher in sterile than in fertile F4 animals, suggesting that sterility is associated with an increase in transcriptional deregulation over time.

To visualize how changes in gene expression progress over subsequent generations in set-2 mutant germlines, we applied Principal Component Analysis (PCA) on the expression matrix of all genes (Figure 1C). As expected from our initial analysis, the distance on the plot relative to the parental samples increased over generations and was highest for samples from F4 sterile animals. In addition, we observed a significant dispersion of the 3 lineages: inter-lineage distance was more pronounced between F4 fertile samples than between P0 or F2 samples, and was highest for F4 sterile samples (Figure 1C, dotted ellipses).

The divergence between samples indicates that for each lineage, the set of misregulated genes is largely distinct, and prompted us to look for a transcriptional signature common to all lineages across generations. For all genes within each lineage we computed the difference in expression (log ratios) between P0 and F2 samples (early changes), F2 and F4 fertile samples (later changes that partially compromise fertility), and F4 fertile and sterile samples (changes associated with fully penetrant sterility). Strikingly, for each of the three lineages PCA on the resulting matrix revealed clustering between P0 and F2 samples and F4 fertile-F4 sterile samples (Figure 1D). By contrast, F2-F4 fertile sets mapped separately in the PCA plot. This provided a key insight: for each of the 3 lineages, a set of genes is misregulated in early generation (P0 to F2), through a process we define and refer to hereafter as “priming,” and further misregulated in the same direction (up or down) in sterile F4 animals.

We next investigated whether a common set of genes undergoes the 2-step deregulation process described above in all 3 lineages: misregulation in P0-F2, and further misregulation in F4 fertile-F4 sterile, in the same direction. We first discarded genes that were poorly expressed in all samples, since their variation may simply be due to background signal, and focused our analysis on highly expressed genes (see “Materials and Methods”). This left a reduced set of 7,238 genes (Supplementary Figure S1C and Supplementary Table S1). We then considered each lineage as an independent class and performed Within-Class Analysis (WCA), a derivative of PCA that takes into account variation due to a given factor (lineage in this case) (Baty et al., 2006). The first component of WCA captured 67% of the Within-Class variation (Figure 1E). Plotting samples along within-class PC1 revealed that transcriptional deregulation progressed in a continuum across generations in all 3 lineages, first during the priming step (P0-F2), and then toward sterility (from F2 to F4 sterile) (Figure 1F). Based on these results, we conclude that genes that contribute to PC1 define a path to sterility. To identify genes that vary along this path, we plotted the distribution of the coefficients that define PC1, which we called “WCA scores” (Figure 1G). In this analysis, negative WCA scores are attributed to genes globally downregulated, and positive scores to genes globally upregulated across generations. Consistent with a bias toward upregulation, WCA scores were asymmetrically distributed, with significantly more genes having a positive score. The most extreme WCA scores (positive or negative) identify genes that contribute the most to PC1, and hence to sterility.

We next examined how transcription deregulation progressed in each lineage: six of the nine genes with the highest negative WCA score were downregulated between P0 and F2 animals, and eight were more strongly downregulated between F2 and sterile F4 than between F2 and fertile F4 (Supplementary Figure S2). Similarly, eight of the nine genes with the most positive WCA score were upregulated between P0 and F2, and all were more strongly upregulated between F2 and sterile F4 than between F2 and fertile F4 (Supplementary Figure S3). Therefore, our WCA analysis identified genes with similar transcriptional trajectories in all three lineages. We arbitrarily set negative and positive thresholds, and hereafter refer to the genes with the largest negative and positive scores as “downregulated contributors” (250 genes) and “upregulated contributors” (1000 genes), respectively (Figure 1G and Supplementary Table S1). These genes make the largest contributions to within-class PC1, and their deregulation is tightly associated with the loss of germline function in all 3 lineages. We collectively refer to these genes as “contributors.” Note that by definition all “contributors” first become deregulated in the priming event. Interestingly, we found a significant overlap between our list of contributors and genes differentially expressed in set-2 mutant germlines at 20°C (Supplementary Figure S4 and Supplementary Table S2)¹. Because at 20°C transgenerational loss of fertility is only observed after 15–20 generations (Xiao et al., 2011; Herbette et al., 2017), these results suggest that the misregulation of these same genes may also contribute to sterility at 20°C, and heat-stress accelerates this process.

Loss of germ-cell identity in C. elegans germlines depleted of PRC2 or P-granule components is associated with decreased expression of germline genes and increased expression of somatic genes (Gaydos et al., 2012; Knutson et al., 2017). Comparison of our lists of contributors with published lists of gene expression profiles in wild-type animals revealed a similar trend in set-2 mutants, suggesting that this transcriptional reprogramming is an early priming event (Figure 2A and Supplementary Table S1; Knutson et al., 2017). Downregulated germline-specific contributors include genes required for meiosis (apc-10, him-5), the wago-4 argonaute-encoding gene, and genes that encode P-granule components (cey-2, cey-3, andpgl-3).

Progressive loss of fertility in set-2 mutants is associated with loss of P-granules (Robert et al., 2014). We found that set-2 inactivation and loss of P-granule components share similar transcriptional signatures (Figure 2B), with a significant overlap for both down- and up-regulated contributors (Supplementary Table S2; p-value = 1.11e-15 and 1.18e-219, respectively). Commonly downregulated genes are enriched in both ubiquitous and germline genes (Supplementary Table S2; p-value = 2.4e-4 and 5.4e-6), while commonly upregulated genes show no particular bias (Supplementary Table S2). Notably, of the 240 soma-specific genes that are upregulated contributors in set-2 mutant germlines, only 133 are also upregulated in P-granule-depleted germlines (Supplementary Table S2; p-value = 0.07). Therefore, SET-2 and P-granules appear to have both common and unique functions in the repression of somatic gene expression in the germline.

Strikingly, we found that 389 of the 941 transcription factors found in C. elegans [wTF3.0; (Fuxman Bass et al., 2016)] were present in the list of 7,238 genes analyzed by WCA, and of these 101 were identified as contributors (1.5-fold enrichment; p = 6e-06) (Figure 2C and Supplementary Table S3). These include the R-SMAD SMA-3, the THAP transcription factor LIN-15B, the RXF family transcription factor DAF-19, 7 CEH homeodomain transcription factors, 22 Nuclear Hormone Receptors (NHRs), and (CEBP-1/C/EBP) (Savage et al., 1996; Swoboda et al., 2000; Antebi, 2006; Ceol et al., 2006; Yan et al., 2009). These genes consistently displayed increased expression between P0 and F2 (in the priming step) and a further increase between F2 and sterile F4 (Figure 2D). Significantly, downstream targets of LIN-15B and CEBP-1 were also identified as upregulated contributors, including 226 of the 452 genes activated in PGCs following ectopic expression of LIN-15B (Figure 2E, p-value = 1.71e-93; Supplementary Table S3; Lee C. S. et al., 2017), and 38 of 112 CEBP-1 target genes (Kim et al., 2016; Figure 2F, p-value = 4.8e-05; Supplementary Table S3). Targets of transcription factors somi-1 and R06C1.6 were also significantly enriched in the list of upregulated contributors [p-value = 0.006 and 0.0002, respectively (Fuxman Bass et al., 2016)], consistent with deregulation of a broad class of transcriptional networks driving loss of germline identity.

Further analysis revealed that upregulated contributors are highly enriched on the X chromosome (Figure 3A, p-value = 6e-193). These X-linked genes have a significantly higher WCA score than autosomal genes (Figure 3B), showing that they strongly contribute to sterility and loss of germline identity. PRC2 components MES-2 and MES-6 cooperate with MES-4 to repress the X chromosome and maintain germline identity (Xu et al., 2001; Bender et al., 2004, 2006; Rechtsteiner et al., 2010; Gaydos et al., 2012). Comparison of genes misregulated in the absence of both mes-2 and mes-4 (mes-2;mes-4 double mutants) and our lists identified 188 commonly upregulated genes (Figure 3C; hypergeometric p-value = 2.34e-97, Supplementary Table S4), mostly located on the X chromosome (137 genes, hypergeometric p-value = 3.2e-21). Conversely, of the 18 commonly downregulated genes, all are found on autosomes (Figure 3D; hypergeometric p-value = 5.84e-14). This is unlikely to reflect regulation of mes-2 or mes-4 expression levels by SET-2, since neither genes was identified by WCA. Similarly, set-2 is not a transcriptional target of PRC2 (Gaydos et al., 2012). Rather, our results suggest that PRC2/MES-4 and SET-2 act in parallel pathways to silence the X chromosome and promote proper expression of a subset of autosomal genes.

Polycomb Repressive Complex 2-dependent repression of the X chromosomes in the germline is counteracted by increased GLP-1/Notch signaling: in worms bearing the gain-of-function allele glp-1(ar202), increased signaling induces expression of specific genes normally repressed by PRC2 (Seelk et al., 2016). We found that glp-1 is over-expressed in all 3 lineages across generations (WCA score = 0.0054), and that GLP-1/Notch targets are upregulated contributors (39 genes, representing a 5.65-fold enrichment; p-value = 1.25e-24, Figure 3E). Of these targets, 30 are located on the X chromosome (p-value = 7.5e-06, Supplementary Table S1), and 18 of these are upregulated in germlines lacking both MES-2 and MES-4 (p-value = 1e-3). For these genes, SET-2 is likely to act redundantly with PRC2/MES-4 through inhibition of GLP-1/Notch (Figure 3G). An additional subset of upregulated contributors on the X chromosome is enriched in targets of the LIN-15B transcription factor (Figure 3F; p-value = 7.71e-9), and depleted for PRC2/MES-4 targets (p-value = 7e-3). Therefore, SET-2 likely represses their expression through inhibition of LIN-15B and independently of PRC2. Altogether, these data suggest that SET-2 acts upstream of transcription factors and regulatory networks that control X-chromosome repression.

So far our analysis identified a set of genes whose deregulation is associated with loss of fertility and germline identity in all 3 lineages analyzed. We reasoned that if increased expression of upregulated contributors is an essential step leading to sterility and loss of germline identity, then decreasing their expression by RNAi over subsequent generations should delay the onset of sterility.

To test the above hypothesis, set-2 mutants at the L4 larval stage were transferred to RNAi feeding plates at 25°C, and fertility scored in adults. In control experiments (no RNAi), the total number of progeny dramatically decreased between F4 and F8, and animals became completely sterile at the F14 generation (Figure 4). The delay in the onset of sterility compared to previous experiments (F4-F7) and (Xiao et al., 2011; Robert et al., 2014) is most likely due to the E. coli food source (HT115 versus OP50 routinely used for culturing), consistent with bacterial diet influencing fertility (Watson et al., 2015; Heestand et al., 2018). RNAi knock-down of 9 out of 29 upregulated contributors we tested, including the transcription factors cebp-1, daf-19, attf-5, and somi-1, significantly and reproducibly delayed the onset of sterility in set-2 mutant animals (Figures 4B,C and Supplementary Table S5A). Loss of fertility and brood size were extremely variable, as previously reported for other mutants showing transgenerational loss of fertility (Yanowitz, 2008; Robert et al., 2014). In late generations, delayed onset of sterility was reflected in a larger number of plates with progeny (fertile plates), as well as an increase in the number of animals per fertile plates. For instance, when animals were grown on RNAi plates targeting cebp-1, most plates (5/6) contained more than 120 animals at the F12 generation. By contrast, only a single control plate (no RNAi) with less than 30 animals was recovered at the F12 generation (Figure 4B). RNAi knock-down of 11 additional genes (see legends Supplementary Table S5A) resulted in early sterility that prevented further investigation of their role (data not shown). Finally, for the remaining nine genes (flh-2, gfi-3, nhr-48, jmjd-3.1, utx-1, kgb-1, puf-9, miz-1, and ncam-1), no significant transgenerational effect on fertility was observed after RNAi treatment. For these genes, the absence of an effect may be due to reduced efficacy of the RNAi treatment, or redundancy.

The above experiments show that our approach successfully identified genes whose increased expression in the germline over generations plays an active role in preserving germline immortality. Because we identified SMA-3 SMAD as an upregulated contributor to loss of germ-cell fate, we asked whether activation of TGF-β signaling in set-2 mutant germlines actively contributes to loss of germ cell fate. Two canonical TGF-β signaling pathways have been described in C. elegans, defined by the DAF-7 and DBL-1 ligands that function through distinct receptors and their pathway-specific SMADs (Supplementary Figure S5A; Savage-Dunn and Padgett, 2017; Gumienny and Savage-Dunn, 2018). Analysis of the WCA list revealed the presence of components of both pathways, including the Type I receptor DAF-1, R-SMADs (SMA-3, DAF-8, and DAF-14), Co-SMADs (SMA-4 and DAF-3), and downstream transcription factors (SMA-9 and DAF-12). By extending the list to include genes with the top 1500 (instead of 1000) WCA scores, we also identified the Type II receptor DAF-4, the extracellular regulator SMA-10, and the transcription factor DAF-5 (Supplementary Table S1). Additional upregulated contributors associated with TGF-β signaling include the KIN-29 serine/threonine kinase (Maduzia et al., 2005) and OBR-3 (Sugawara et al., 2001). Remarkably, for all of these, transcript levels increased in all 3 lineages between P0 and F2 (Figure 2D and Supplementary Figure S5B). Furthermore, for at least 2 of the 3 lineages, expression showed a more significant increase between F2 and F4 sterile than F2 and F4 fertile animals. RNAi knock-down of sma-3, sma-9, daf-5, sma-10, kin-29, and obr-3 in set-2 mutants resulted in a significant and reproducible delay in the onset of sterility (Figure 5A and Supplementary Table S5B), consistent with ectopic activation of TGF-β signaling actively contributing to loss of germ cell fate.

To confirm by an independent assay that components of TGF-β signaling are upregulated in the germline, we performed single molecule RNA-FISH (smFISH) on dissected gonads (Ji and van Oudenaarden, 2012; Lee C. S. et al., 2017). In agreement with our RNA-seq data, daf-5, and sma-3 transcripts were detected in most germ cells of F4, but not P0 animals (Figure 5B and Supplementary Figure S6A).

In wild-type animals, TGF-β signaling regulates multiple pathways, including somatic development, oocytes quality, and germline proliferation [(Ren et al., 1996; Foehr et al., 2006; Luo et al., 2009; Dalfo et al., 2012; Pekar et al., 2017) (Savage-Dunn and Padgett, 2017; Gumienny and Savage-Dunn, 2018)] through downstream target genes. To test whether downstream components of TGF-β signaling are misexpressed in set-2 mutant germlines, we compared our list to published lists of TGF-β target genes at different developmental stages (Luo et al., 2010). To our surprise, upregulated contributors are enriched in genes upregulated in oocytes from sma-2 mutant animals, in which the pathway is inactive (Luo et al., 2010), as well as a largely non-overlapping set of genes upregulated in sma-2 mutant animals at the L4 stage, prior to oocyte development (Supplementary Figures S6B,C). A total of 137 upregulated contributors are genes normally repressed by TGF-β signaling in oocytes of wildtype animals (p-value = 5.00e-08), or in L4 (p-value = 1.6e-32). This set of upregulated contributors is also enriched in soma-specific genes and LIN-15B activated genes (Supplementary Figure S6D), suggesting that LIN-15B and TGF-β may cooperate at a subset of loci to promote gene expression during loss of germ cell fate. Altogether, these results suggest a complex network whereby ectopic activation of TGF-β pathway components in the germline alters expression of target genes by mechanisms other than canonical signaling.

Nematode strain maintenance was described previously (Brenner, 1974). The strain set-2(bn129) III/qC1(qIs26) used in this study is described in Herbette et al. (2017). Strain set-2(bn129)/qC1(qIs26) was maintained at 20°C and set-2(bn129)/qC1(qIs26) animals were switched to 25°C one generation before starting experiments. During the experiments, animals were maintained at 25°C.

Gonads from set-2 homozygous worms raised at 25°C? were collected from 3 independent lineages derived from set-2(bn129)/qC1(qIs26) [P0 set-2(+)/set-2(−)] animals by initially picking 25 set-2(+)/set-2(−) L4 worms and then cloning single worms from subsequent generations to individual plates. Fertility of each worm was assessed before dissection by monitoring egg laying for 24 h after young adult stage. Worms that laid any live progeny were considered as fertile. For the F4 sterile samples, only worms lacking visible embryos in the uterus were collected. For each replicate, seventy five to one-hundred gonad arms were dissected from set-2(bn129)/qC1(qIs26) [P0 set-2(+)/set-2(−)] and set-2(bn129)/set-2(bn129) [F2 and F4 set-2(−)/set-2(−)] animals. Dissected gonads were cut at the gonad bend with 30 1/2-gauge needles in egg buffer (pH 7.3, 27.5 mM HEPES, 130 mM NaCl, 2.2 mM MgCl₂, 2.2 mM CaCl₂, and 0.528 mM KCl) containing 0.5% Tween 20 and 1 mM levamisole and collected into Trizol. Total RNA was extracted and ribosomal RNA was depleted using an NEBNext rRNA Depletion Kit (Human/Mouse/Rat) (catalog number E6310). Libraries were constructed using an NEBNext Ultra RNA Library Prep Kit for Illumina sequencing (catalog number E7530) and sequenced at the Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley, using Illumina HiSeq 2500 and 4000 platforms. For differential expression analysis, raw sequences were first mapped to transcriptome version WS220 using TopHat2 (Kim et al., 2013). Only reads with one unique mapping were allowed, otherwise default arguments were used. Reads mapping to ribosomal RNAs were then removed. HTSeq (Anders et al., 2015) was used to build a count table of expression levels per transcript. DESeq2 (Love et al., 2014) was used to normalize read counts across samples (Supplementary Figure S1B), leading to matrix M that contained the normalized expression levels of 20261 genes. DESeq2 was also used to determine genes for which expression at a downstream generation differed from expression at P0, with no consideration for lineage-dependent relations (Figure 1C). P-values were adjusted for multiple testing using the Benjamini–Hochberg method as implemented in DESeq2.

We did all further multivariate analysis using R version 3.4.4². To perform PCA (Figure 1C), we first transformed M by applying x = Log₂(x + 1) to all of its values, obtaining matrix LM, which we transposed and processed with the prcomp() function. To perform PCA on expression differences (Figure 1D), we constructed a 20261 × 9 matrix (DLM), where rows were genes and columns were expression differences of interest (e.g., F2-P0 lineage 1) which we computed by subtracting one column of LM (e.g., F2 lineage 1) from another (e.g., P0 lineage 1). We then transposed DLM and processed it with the prcomp() function. Before applying WCA, we first discarded genes that were poorly expressed in all samples, because their variation may simply be due to varying background signal. To determine a meaningful threshold of minimal expression, we plotted the distribution of all genes according to (i) their maximal expression level in all samples and (ii) their standard deviation of expression across all samples, which could be high due to either meaningless background signal variation or meaningful biological differences between samples. The resulting 2-dimensional density plot revealed 2 subpopulations of genes (Supplementary Figure S1C): a major population that was poorly expressed in all samples, and a secondary population with high expression in at least one sample. For the poorly expressed genes, standard deviation increased monotonously with their maximal expression level, as expected for background variation. A different pattern was observed for highly expressed genes (maximal expression greater than 5). First, their standard deviation of expression was not correlated with maximal expression level. This is expected if change in expression is dependent on SET-2 loss and not on background variation. Second, standard deviation was low for the majority of these genes and high for a small subset. This observation is expected if SET-2 acts only on a subset of genes. We therefore chose an arbitrary cutoff of minimal expression that discarded the first subpopulation (red line in Supplementary Figure S1C), leaving a reduced set of 7,238 genes having at least one value greater than 5 (Supplementary Table S1). We then processed the resulting matrix with the dudi.pca() function of the ade4 package [63] (version 1.7-6)³ using parameters scan = FALSE, scale = FALSE, and nf = 4. We then processed the resulting object with the ade4:wca() function using lineages as classes and parameters scan = FALSE, nf = 2. The first column of attribute c1 of the resulting object corresponded to the WCA scores reported in text and figures.

Bacterial clones containing RNAi feeding vectors were collected in the C. elegans RNAi collection (made by J. Ahringer, Source Bioscience). The molecular sequence of insert present in each RNAi clones was checked by sequencing (primer ggtcgacggtatcgataagc) after PCR amplification (single primer in duplicated T7 promoter taatacgactcactataggg) performed directly on colonies. For feeding, bacterial clones were amplified 18 h at 37°C in LB complemented with Ampicilline (50 μg/ml). Transcription (from duplicated T7 promoters) was induced by adding IPTG (1 mM final) and growing liquid cultures for 2 additional hours at 37°C. 200 μl of induced cultures were plated on NGM plates complemented with IPTG (1 mM). At each generation, 6 L4 animals were transferred on RNAi plates and grown at 25°C for 3–4 days. Progeny was briefly counted every two generations to evaluate animas fertility. Animals were considered as sterile when less than 20 progeny were present on a plate.

Single molecule RNA-FISH was performed as described in Lee C. S. et al. (2017), but germline dissections and all subsequent steps were performed on Poly-Lysine coated coverslips instead of test tubes to minimize loss of material. For dissection, worms were placed in a 15 μl drop of dissection buffer, while other steps were carried out in 30 μl of the corresponding buffer.