RNA-seq analysis of germline stem cell removal and loss of SKN-1 in C. elegans

In C. elegans, ablation of germline stem cells (GSCs) extends lifespan, but also increases fat accumulation and alters lipid metabolism, raising the intriguing question of how these effects might be related. Here we show that a lack of GSCs results in a broad transcriptional reprogramming, in which the conserved detoxification regulator SKN-1/Nrf increases stress resistance, proteasome activity, and longevity. SKN-1 also activates diverse lipid metabolism genes and reduces fat storage, thereby alleviating the increased fat accumulation caused by GSC absence. Surprisingly, SKN-1 is activated by signals from this fat, which appears to derive from unconsumed yolk that was produced for reproduction. We conclude that SKN-1 plays a direct role in maintaining lipid homeostasis, in which it is activated by lipids. This SKN-1 function may explain the importance of mammalian Nrf proteins in fatty liver disease, and suggests that particular endogenous or dietary lipids might promote health through SKN-1/Nrf. Overall design: Samples were prepared from ~5,000 synchronized, L1 arrested day-one adult animals cultured at 25°C. Worms were synchronized by sodium hypochlorite (bleach) treatment, as previously described (Porta-de-la-Riva et al., 2012). Bleach solution (9 mL ddH2O; 1 mL 1 N NaOH; 4 mL Clorox bleach) was freshly prepared before each experiment. Worms were bleached for 5 minutes, washed 5x in M9, and arrested at the L1 stage at 25°C in M9 containing 10 µg/mL cholesterol. Feeding RNAi was started at the L1 stage. This approach only partially reduces skn-1 function, but allows analysis of larger samples than would be feasible with skn-1 mutants, which are sterile (Bowerman et al., 1992). Because these animals were not treated with FUdR, the WT adults contained an intact germline and eggs. As is explained in the Results section, we therefore confined our analysis to genes that were overrepresented in glp-1(ts) animals, which lack eggs and most of the germline, and established a high-confidence cutoff for genes that were upregulated by GSC absence as opposed to simply being expressed specifically in somatic tissues. RNA was extracted using the same protocol for qRT-PCR samples. Purified RNA samples were DNase treated and assigned a RIN quality score using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). Only matched samples with high RIN scores were sent for sequencing. Single read 50 bp RNA sequencing with poly(A) enrichment was performed at the Dana-Farber Cancer Institute Center for Computational Biology using a HiSeq 2000 (Illumina, San Diego, CA). FASTQ output files were aligned to the WBcel235 (Feb 2014) C. elegans reference genome using STAR (Dobin et al., 2013). These files have been deposited at the Gene Expression Omnibus (GEO) with the accession number GSE63075. Samples averaged 75% mapping of sequence reads to the reference genome. Differential expression analysis was performed using a custom R and Bioconductor RNA-seq pipeline (http://bioinf.wehi.edu.au/RNAseqCaseStudy/) (Gentleman et al., 2004; Anders et al., 2013; R Core Team, 2014). Quantification of mapped reads in the aligned SAM output files was performed using featureCounts, part of the Subread package (Liao et al., 2013, 2014). We filtered out transcripts that didn't have at least one count per million reads in at least two samples. Quantile normalization and estimation of the mean-variance relationship of the log-counts was performed by voom (Law et al., 2014). Linear model fitting, empirical Bayes analysis and differential expression analysis was then conducted using limma (Smyth, 2005). To identify genes that are upregulated in a SKN-1-dependent manner by GSC loss, we sought genes for which glp-1(ts) expression was higher than WT, and for which glp-1(ts);skn-1(-) expression was reduced relative to glp-1(ts). To test for this pattern, if a gene's expression change was higher in the comparison of glp-1(ts) vs. WT and lower in the comparison of glp-1(ts);skn-1(-) vs. glp-1(ts), then we calculated the minimum (in absolute value) of the t-statistics from these two comparisons, and assessed the significance of this statistic by comparing to a null distribution derived by applying this procedure to randomly generated t-statistics. We corrected for multiple testing in this and the differential expression analysis using the false discovery rate (FDR) (Benjamini and Hochberg, 1995). Heatmaps were generated using heatmap.2 in the gplots package (Warnes et al., 2014). Functional annotations and phenotypes were obtained from Wormbase build WS246. SKN-1 transcription factor binding site analysis of hits was conducted with biomaRt, GenomicFeatures, JASPAR, MotifDb, motifStack, MotIV, and Rsamtools (Sandelin et al., 2004; Durinck et al., 2005; Durinck et al., 2009; Lawrence et al., 2013; Ou et al., 2013; Mercier and Gottardo, 2014; Shannon, 2014). JASPAR analysis was performed with the SKN-1 matrix MA0547.1 using 2 kb upstream sequences obtained from Ensembl WBcel235 (Staab et al., 2013). modENCODE SKN-1::GFP ChIP-seq analysis of hits was performed using biomaRt, ChIPpeakAnno, IRanges, and multtest (Durinck et al., 2005; Durinck et al., 2009; Gerstein et al., 2010; Zhu et al., 2010; Niu et al., 2011; Lawrence et al., 2013). SKN-1::GFP ChIP-seq peaks were generated by Michael Snyder's lab. We used the peak data generated from the first 3 larval stages: L1 (modENCODE_2622; GSE25810), L2 (modENCODE_3369), and L3 (modENCODE_3838; GSE48710). Human ortholog matching was performed using Wormbase, Ensembl, and OrthoList (Shaye and Greenwald, 2011). Gene lists were evaluated for functional classification and statistical overrepresentation with Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.7 (Dennis et al., 2003).

12

Steinbaugh MJ, Narasimhan SD, Robida-Stubbs S, Moronetti Mazzeo LE, Dreyfuss JM, Hourihan JM, Raghavan P, Operaña TN, Esmaillie R, Blackwell TK