PHA-4/FOXA-regulated microRNA feed forward loops during Caenorhabditis elegans dietary restriction

Next generation sequencing reveals eat-2(ad1116) miRNA profile

In order to gain an insight into the complex interplay of miRNAs during DR, we performed NGS analysis to compare miRNA profiles of wild-type (WT) N2 Bristol and eat-2(ad1116). The eat-2(ad1116) strain lives considerably longer than WT; however, none of the strains have significant mortality at day 1 (equivalent to young-adults) or day 8 of adulthood (aging worms) (Figure S1). The eat-2 mutants have defective pharyngeal pumping and consume less food since hatching. We hypothesized that the system in eat-2(ad1116) is already geared up for long life and would have significantly modified post-transcriptional cellular environment on day 1 of adulthood i.e., at the young-adult stage. Additionally, we wanted to study the changes that occur later in their life, as compared to WT, in order to understand the late-life cellular response to DR. Synchronized worms were harvested on day 1 or 8 of adulthood and small RNA sequencing libraries prepared (see experimental procedures). Each sample was run in a single lane of a flow cell in the Genome Analyzer IIx (GAIIx) platform (Illumina Inc., USA) that led to the generation of around 20 million reads per lane with an average trimmed read size of ~22 nucleotides (Figure S2A). All reads mapped to the C. elegans genome while 81-95% readily mapped to the miRBase release 19 that has a collection of 368 worm miRNAs (Figure S2B, 1). We detected 185 miRNA in day 1/young-adult worms and 187 in day 8/aging WT samples, while 224 and 161 miRNA were detected in eat-2(ad1116) samples in day 1 and day 8 samples, respectively, all miRNA having greater than 10 read counts (Table S1). We were also able to detect between 0.5-1.2 million reads that mapped to known 21U-RNAs (Table S1). After eliminating rRNA, tRNA and ncRNA reads, the remaining un-annotated reads were processed for novel miRNA discovery, as discussed below.

We compared the miRNA profiles of WT worms collected at days 1 or 8 with a previously published study on changes in miRNA profile during aging in WT nematode [31] (Table S2). Out of a total of 21 miRNA reported to be upregulated in day 8 compared to day 1 in their study, 16 matched to our data. We found an additional 20 miRNA upregulated in day 8 samples. When we compared the downregulated miRNAs from the two studies, 7 out of 18 matched our data while we obtained an additional 24 miRNAs, probably because of increased depth of sequencing. It may also be due to the less number of miRNAs (174) present in the miRbase release 14, to which the earlier data was mapped. Together, the sequencing depth and coverage attained in our study includes most of the known miRNAs and has the potential to discover new and low expressing ones.

Unique changes in miRNA profile in eat-2(ad1116) worms

We compared the normalized miRNA expression (TPM, see experimental procedures) of WT with eat-2(ad1116) either in day 1/young-adults or in day 8/aging worms and found that most of the miRNAs were distinctly upregulated in young-adult eat-2(ad1116) (Figure 1A, Figure S3 for analysis flow chart, Table 1, Table S3). For this, we considered all miRNAs that had read counts over 10 in both the strains on a particular day, with minimum fold change of ±1.5 having a p<0.05. Out of a total of 184 miRNA common between young-adult WT and eat-2(ad1116) on day 1, 105 were significantly upregulated in eat-2(ad1116). Interestingly, none were significantly downregulated. Seventy nine miRNA had no significant change in expression pattern. We found that 40 miRNA were exclusively expressed in eat-2(ad1116) (Table S4); out of them cel-miR-41-5p, cel-miR-40-5p, cel-miR-39-5p, cel-miR-4813-5p had read counts more than 100. For further bioinformatic analysis using these upregulated miRNAs in day 1/young-adult worms, we did not consider the strandedness of the miRNAs and only analysed 82 out of the 105 as uniquely upregulated miRNAs.

In order to validate the sequencing results, we chose 9 miRNAs having varying expression levels. Using quantitative real time PCR, we could validate the expression pattern of eight of these miRNAs; they were significantly upregulated in eat-2(ad1116) compared to WT on day 1 (Figure 1C).

Next we compared the miRNA profiles of aging WT and eat-2(ad1116) worms collected at day 8 of adulthood (Figure 1B, Table S3). We found 157 miRNAs common between WT and eat-2(ad1116); 15 were significantly upregulated in eat-2(ad1116) while 43 were downregulated compared to WT. Twenty miRNAs had no significant change in expression pattern. Four miRNAs were exclusively expressed in eat-2(ad1116) while 30 were present only in WT; out of them 4 had read counts more than 100 (Table S5).

From the above observations, it appeared that many miRNAs that are upregulated in young-adult eat-2(ad1116) on day 1 returned to the levels of WT or lower with progressing age. When we compared the expression patterns of miRNAs that were expressed in all the four samples with a minimum read count of 10, we found that 61 of the 152 (~40%) were upregulated in eat-2(ad1116) on day 1 but returned to near WT levels on day 8 (Cluster A in Figure 2A; Table S6). These included miRNAs cel-miR-70, cel-miR-51, cel-miR-56, cel-miR-230, cel-miR-238, cel-miR-83, etc. There were 10 others that were upregulated in day 1 of eat-2(ad1116), but remained upregulated even in day 8 (Cluster B). These included miRNAs cel-miR-237, cel-miR-238, cel-miR-229, cel-miR-791 and cel-miR-72. Levels of miRNAs cel-miR-241, cel-miR-53, cel-miR-57 and cel-miR-228 were decreased in day 8 in WT compare to day 1 levels, but failed to so in eat-2(ad1116) (Cluster E). Several miRNAs had expression patterns that were reverse of that in WT worms. For e.g., 11 miRNA, including cel-miR-63, cel-miR-53, cel-miR-36-41, cel-miR-254 and cel-miR-61 increased in expression from day 1 to 8 in WT worms; however, in eat-2(ad1116) it had an inverse expression pattern (Cluster C). As a corollary, only cel-miR-58 had the opposite expression pattern-decreased in WT but increased in eat-2(ad1116). Apart from miRNAs expressing in all samples, another 79 miRNA had more than 10 read count in any one or more of the 4 samples and were categorized separately (Table S7). Most of these miRNAs had low abundance as indicated by their read counts that were less than 1000. Only cel-miR-86 had 1258 read counts (TPM 58.72) on day 1 but had no detectable count on day 8. In this list also, we found 29 miRNAs were exclusively expressed in eat-2(ad1116) on day 1. These observations point to the profound changes in miRNA expression during eat-2(ad1116)-mediated dietary restriction, possibly to support metabolic changes associated with the long life span. MicroRNAs were, in general, found to be upregulated in eat-2(ad1116) young adult worms on day 1 of adulthood.

The miRNAome of eat-2(ad1116) worms is different from the long-lived insulin-like signalling mutant

Next we asked whether the profile obtained with eat-2(ad1116) is a general pattern for all long-lived mutants. In worms, the insulin-like signalling pathway is a major regulator of longevity and development [32]. Mutations in the insulin-like receptor daf-2 leads to dramatic increase in life span [33]. However, the mechanisms by which DR and insulin-like signalling increase life span are independent of each other [6]. We profiled the miRNAome of daf-2(e1370) at day 1 of adulthood and found that only 24 miRNAs were upregulated significantly while 33 were downregulated (Table S1, S18). All 24 upregulated miRNAs were also commonly detected in eat-2(ad1116) young adult worms, although the levels were higher in the latter in most cases (Figure 2B). In stark contrast to daf-2(e1370), none of the miRNAs were downregulated in eat-2(ad1116) as mentioned above. Together these data showed that the profile of miRNA in eat-2(ad1116)-mediated DR worms is unique and does not represent a common trend of long-lived mutants. However, some of the commonly upregulated miRNAs may still correspond to a signature for long-lived mutants, but needs to be validated in other strains that have enhanced longevity.

Novel miRNA that change expression in eat-2(ad1116)

Since DR induces a dramatic reprogramming of development and metabolism, we hypothesized that novel miRNAs may contribute towards this effect. Additionally, with increased depth of sequencing, we expected to detect novel miRNAs in our WT samples also. We focused our attention on the 877496 (in young adults on day 1) and 682280 (in day 8/aging worms) non-redundant un-annotated reads that aligned to the C. elegans genome in eat-2(ad1116) (Table S1). In order to discover novel miRNAs, we used the miRNA discovery package miRDeep2 (see Supporting information for reference). We found 48 novel miRNAs in young-adults and 9 in aging eat-2(ad1116) worms, with a read count greater than 10. Additionally following analysis of unassigned reads, in WT we identified 21 in young-adults and 35 in aging worms (8 and 9). Out of the 13 commonly predicted miRNAs between day 1/young-adult WT and eat-2(ad1116) (without read count or p value cut-off) worms, 9 were upregulated and 3 down-regulated as calculated on the basis of read count (Table S10). On the other hand, 2 out of 3 common miRNAs on day 8 were downregulated in eat-2(ad1116). Thus, similar to the known miRNA, we found distinct upregulation of novel miRNA population in young-adult eat-2(ad1116) worms.

In order to select predicted novel miRNA for experimental validation, we used one additional criterion. Apart from read count filter of ≥10 and presence of a star sequence, we used MiRDeep Score >1 as a cut-off. Using these criteria, we identified 7 predicted miRNA in WT and 15 in eat-2(ad1116) young-adult worms on day 1, and 5 in aging WT as well as 1 in aging eat-2(ad1116). From these, we selected 3 predicted novel miRNA for validation (Figure 3A, Table S11). We grew WT, alg-1(gk214) or alg-2(ok304) and checked the expression of the predicted miRNAs at young-adult stage. The ALG-1 and ALG-2 proteins are required for miRNA maturation [34]. We found that all of the candidate miRNAs expressed in an alg-2-dependent manner while 2 expressed in alg-1-dependent manner as determined by quantitative real-time PCR analysis (Figure 3B), confirming that they are indeed bona fide miRNAs. Out of these three novel miRNA, two are down-regulated in eat-2(ad1116) young-adult worms while one is upregulated when compared to WT (Figure 3C), matching their expression as determined by sequencing. Together, it appears that the profound transcriptional and post-transcriptional changes required following DR in eat-2(ad1116) involves the role of additional novel miRNAs.

PHA-4 regulates majority of the upregulated miRNA in eat-2(ad1116)

The FOXA transcription factor PHA-4 is an absolute requirement for DR-induced longevity in eat-2 mutants [11]. We asked whether PHA-4 could be directly regulating the expression of the miRNAs upregulated in eat-2(ad1116). For this, we used ChIPBase (See Supporting information for reference), a database that catalogues transcription factor binding maps on miRNA (as well as other small RNAs) and protein-coding gene promoters using ChIP-seq data from published work [29, 30]. We queried the database with all the miRNAs that were upregulated in young-adult eat-2(ad1116) worms collected on day 1. We found that 65 of the 82 promoters of upregulated miRNA have one or more PHA-4 binding peaks (Table 1), potentially pointing at direct transcriptional regulation of these miRNA by PHA-4. In order to validate this data, we chose 9 of these miRNAs that are induced in eat-2(ad1116) as compared to WT and checked their expression in eat-2(ad1116) mutants, both in presence or absence of PHA-4. We found that all these miRNAs, except cel-mir-47, are regulated in a PHA-4-dependent manner (Figure 4A). Thus, PHA-4 not only controls mRNA levels of coding genes, like sod-2 and sod-4 during DR [11], it possibly regulates the transcription of a large part of the miRNA population.

DR-induced miRNAs in eat-2(ad1116) may target a large part of the transcriptome

We discovered dramatic changes in miRNA expression levels with DR in eat-2(ad1116). We also identified PHA-4 as a transcription factor responsible for regulating most (65 out of 82) of these miRNAs.

MiRNAs target the 3′UTR of mRNA and lead to translational arrest or degradation depending on the extent of base pairing in the ‘seed’ sequence. To throw light on the biological processes being regulated by the PHA-4-dependent miRNAs in eat-2(ad1116), we predicted their target genes. Since the prediction of miRNA target is not always accurate, we used two different programs, miRanda as well as Targetscan (see Supporting information for references) and considered only the ones that came up with both the searches (Table S12). The two programs together predicted a total of 5145 possible unique targets of the 65 upregulated miRNAs (Table S13). We found that while miR-71 had the largest number of predicted targets (782), miR-243 has 34 targets. We also found that multiple miRNAs target the same mRNA (Table S13). For example, syg-1 is targeted by 18 miRNA, ptc-1 by 16, whereas ima-3 and ina-1 are targeted by 13 and 15 of the upregulated miRNA, respectively. SYG-1 is required cell autonomously in the HSNL neuron to determine synaptic specificity with vulval muscles and VC neurons [35]. PTC-1 is a Drosophila patched homolog, required for germline cytokinesis [36]. INA-1 is required during development for migration of neurons, coelomocyte precursor and distal tip cells of somatic gonad [37] while IMA-3 is the worm homolog of importin alpha involved oogenesis [38]. These genes may therefore be potentially important targets for miRNAs during DR.

Mature miRNAs, in most cases, bind to the 3′UTR of target mRNAs leading to their translational arrest [15]. Under this premise, we would expect to see down-regulation of protein levels corresponding to these predicted targets. In a previous study [39], the complete proteome of WT and eat-2(ad1116) were compared using quantitative proteomics. Out of the 167 proteins that were significantly downregulated in eat-2(ad1116), we found that 74 of them are predicted targets of the miRNAs (hypergeometric analysis, p=2.67e-8; Table S14). This suggests that PHA-4-regulated miRNAs may translationally inhibit many of their target mRNAs.

Many genes are upregulated in eat-2(ad1116) in a PHA-4-dependent manner

Next, in an effort to study the transcriptional response to DR, we compared the transcriptome of young-adult WT and eat-2(ad1116) worms collected on day 1 of adulthood, using NGS. We generated ~ 41 million reads for WT and ~42.2 million reads for eat-2(ad1116) using the GAIIx platform, as described in the experimental procedure. About 96.4% of WT reads and 97.2% of eat-2(ad1116) reads mapped to the C. elegans genome. We found that similar to the miRNA results, many more mRNAs were upregulated (3607 mRNAs) in eat-2(ad1116) compared to those that were downregulated (231 mRNAs)(Table S15).

Since PHA-4 plays such a central role in DR-mediated longevity, we asked whether these upregulated genes are under the direct control of PHA-4. For this, we looked for PHA-4 binding in the 5 kb promoter proximal region of the genes upregulated in eat-2(ad1116) using ChIPbase. We found that out of the 3607 upregulated transcripts, 2037 have at least one PHA-4 binding site (p=6.18e-171, by Hypergeometric test)(Table S15). To validate whether the presence of binding sites correlate with PHA-4 dependence in transcription, we randomly chose 10 genes. Using qRT-PCR, we found that 8 transcripts are dependent on PHA-4 (Figure 4B). These data suggests that PHA-4 may bind and upregulate a large number of genes in eat-2(ad1116); this may be one of the reasons why knocking down PHA-4 completely suppresses DR-induced life span.

Feed forward loops involving PHA-4 along with its target mRNA and miRNA

We have found that the FOXA transcription factor PHA-4 may be a critical requirement in upregulation of mRNAs as well as miRNA in eat-2(ad1116). Next, we asked whether there was an overlap between genes that are transcriptionally upregulated by PHA-4 and the targets of PHA-4-regulated miRNA. By comparing 2037 PHA-4-regulated genes to 5145 genes that are predicted targets of PHA-4-regulated miRNA, we found a significant overlap (1073 genes with a p = 3.16e-176 by Hypergeometric test; Table S16). Thus, PHA-4 forms a large number of incoherent feed forward loops [40] during DR in eat-2(ad1116) (Figure 5A, B). While on one hand, PHA-4 transcriptionally upregulates the expression of a many transcripts, the same set of transcripts are potentially targeted by miRNAs that it transcribes (Figure 5B). In fact, some of the target genes are controlled by upto 20 miRNAs indicating that they might potentially play important roles during DR.

Next we asked whether the insulin-like signalling pathway also uses such FFLs to control gene expression. The insulin-like signalling pathway components phosphorylate and inactivate the FOXO transcription factor DAF-16 that is a major output of this cascade [32]. Under low insulin signalling conditions, such as in daf-2(e1370), DAF-16 translocates into the nucleus and transactivates a large assortment of genes [32]. Using ChIPbase, we first determined that 22 of the 24 miRNA upregulated in daf-2(e1370) are potential direct targets of DAF-16. Next, using a published daf-2(e1370) transcriptomics data [41], we found that 687 transcripts are upregulated in daf-2(e1370); out of these 336 have DAF-16 binding sites within 5 kb of promoter, according to ChIPbase. Thereafter, we overlapped the predicted targets of the 22 DAF-16-regulated miRNA (2235 transcripts) with that of the DAF-16-dependent mRNA. We found that only 47 genes overlapped (13.98%; p=0.01 by hypergeometric test) in contrast to that of eat-2(ad1116) where PHA-4 regulated 52.67% (p = 3.16e-176) of the upregulated transcripts by FFLs. This data suggested that PHA-4-mediated FFLs are a unique feature of DR-induced longevity in eat-2(ad1116) and may not be used extensively by DAF-16/FOXO downstream of the insulin-like signalling pathway.

PHA-4 target genes regulated through FFLs control important biological functions

Next we determined the functions of all the genes that are potentially transcriptionally as well as post-transcriptionally regulated by PHA-4. Since these genes are under such tight regulation, we can expect them to be important for DR-mediated longevity. We categorized these 1073 genes based on their Gene Ontology (GO) terms as well as determined which molecular pathways these genes function in. We used “Gene Ontology enRIchment anaLysis and visualization” tool (GOrilla) and visualized the results using Revigo. We found that the GO terms associated with metabolism, protein folding and response to various chemical stimuli are enriched (Figure 6A). These genes are involved in ubiquitin mediated protein degradation, lysosomal autophagy, inositol phosphate metabolism, phosphatidylinositol signalling, ErbB signalling as determined by KEGG pathway analysis using DAVID (see Supporting information for references) (Figure 6B). These 1073 genes are also highly connected with each other through protein-protein interaction as determined by STRING database (Figure 5C). Further, these genes are enriched for transcription factors (~10% of the 1073 genes code for transcription factors as against 4.4% found in the genome; p=1.68e-12 by Hypergeometric test). Interestingly, when we compared this list of genes with the GeneAge database (http://genomics.senescence.info/genes/), we found that 86 genes are already known to be important in longevity determination having p=1.66e-12 (Table 2). Together, the genes that are regulated by PHA-4 directly or indirectly through miRNAs most likely are important players in DR-mediated longevity in eat-2(ad1116).

Strain maintenance

All strains were maintained at 20°C using standard C. elegans techniques. All RNAi experiments were initiated using synchronized L1 stage worms. Strains used in the study are: N2 Bristol as WT, eat-2(ad1116), rrf-3(pk1246), eat-2(ad1116);rrf-3(pk1246), alg-1(gk214) and alg-2(ok304).

For collecting worms for NGS analysis, synchronized L1 larvae were placed on Nematode Growth Media (NGM) agar plates and allowed to grow till L4 stage. The plates were overlaid with FuDR to stop the hatching of the eggs. Worms were collected from the plates by washing with 1 X M9 the following day (day 1/young-adult worms) or on day 8 (aging worms). The worm pellet was washed three times with 1 X M9 before Trizol reagent (Invitrogen, USA) was added and RNA isolated as described below.

RNA isolation

RNA isolation was performed using Trizol. Briefly, worms grown on OP50 bacteria, vector or RNAi of interest were washed off the plates with M9 buffer. Thereafter, 0.4 ml of Trizol reagent was added and the worms lysed by vigorous vortexing. RNA was purified by phenol:chloroform:isoamylalcohol extraction and ethanol precipitation. RNA integrity were confirmed by analysis on an Agilent 2100 bioanalyzer, using the RNA 6000 nano kit (Agilent) to confirm that the RIN numbers are above 8.0. Alternatively, the quality of the ribosomal RNA 28 S and 18 S as determined on an agarose gel was used as a measure of integrity and the absorbance at 260/280 nm was used to determine quantity.

RNA library preparation, Next Generation sequencing and data analysis. The small RNA libraries were constructed using Small RNA sample preparation kit v1.5 according to the manufacturer's instructions (Illumina Inc., USA). Briefly, the total RNA (2μg) were ligated to 3′ RNA adapter using RNA ligase truncated and 5′ adapter using T4 RNA ligase 2 (New England Biolabs, USA). The ligation products were reverse transcribed using Superscript II Reverse Transcriptase (Life Technologies, USA) and amplified with 12 cycles of PCR. The PCR products constituting the small RNA cDNA libraries were resolved on 6% Novex TBE PAGE Gel (Invitrogen, USA) and ~150 bp fragments excised. The library was eluted from the acrylamide gel and analyzed on Agilent 2100 Bioanalyzer using DNA high sensitivity kit (Agilent Technologies, USA). NGS of cDNA libraries were performed using Illumina GAII_X for 36 cycles. A total of 2.8 GB of raw sequence data, comprising of WT and eat-2(ad1116) strains at both day 1 and day 8, was imported into the CLC Genomics Workbench 6.5.1 (CLC Bio, Denmark). The reads were trimmed off the adapter sequences and reads containing low quality bases eliminated. The trimmed raw sequences were mapped to the miRBase release 19, allowing for a maximum of two gaps or mismatches. Unpaired group comparisons, based on Transcript Per Million (TPM), were used as expression values to compare different samples. A fold change ±1.5 with a minimum read count of ≥10 were used to filter the differentially expressed miRNA. The p value cutoff was set at p≤0.05 based on Kal's Z test statistical. The sequencing data is available at GEO repository with Series record number GSE60155 and GSE61112.

RNA-Sequencing (RNA-seq) libraries of WT and eat-2(ad1116) Day 1 samples were prepared as recommended by the Illumina TruSeq™ RNA Sample Preparation kit using Low-Throughput (LT) Protocol (Illumina, Inc.,USA). NGS of libraries was performed using Illumina GA II_X for 78 cycles including 6 additional cycles for index read. Sequence reads were aligned using CLC Genomics Workbench 6.5.1 with default setting against C. elegans genome assembly (WS231). Unpaired group comparisons, based on RPKM (Reads Per Kilobase per Million mapped reads), were chosen as expression values for comparing the samples. A fold change ±2.0 and P value ≤0.05 (Kal's Z test) were used to filter the differentially expressed genes.