Deep sequencing identifies circulating mouse miRNAs that are functionally implicated in manifestations of aging and responsive to calorie restriction

Analysis of circulating miRNA sequencing reads

To investigate the potential effects of aging and calorie restriction (CR) on the circulating levels of miRNAs, we used small RNA-Seq (Illumina reads of 50 nt) to analyze the serum levels of miRNAs from 3 young (7-month) and 3 old (27-month) control mice, and 3 old (27-month) mice subjected to CR from 1 month of age. The raw sequencing reads of circulating miRNAs from the different experimental groups were pre-processed and analyzed with miRDeep2 [47]. After removing the 3' adapter sequences, discarding reads shorter than 18 nucleotides, and aligning the processed reads to the mouse genome, miRDeep2 detects both known and novel miRNAs. The miRDeep2 algorithm is based on the miRNA biogenesis model; it aligns reads to potential hairpin structures in a manner consistent with Dicer processing, and assigns scores that represent the probability that hairpins are true miRNA precursors. In addition to detecting known and novel miRNAs, miRDeep2 estimates their abundance.

Discovery of novel miRNAs circulating in the mouse serum

miRDeep2 predicted 79 potential novel miRNAs, and detected 553 known miRNAs, at the relatively stringent score cut-off of 4 and signal-to-noise ratio of 19.6 (Table 1). After filtering of the predicted novel miRNAs by removal of loci matching other RNA genes, and keeping only novel miRNAs with significant randfold p-value, the list was reduced to 23 candidate novel miRNAs (Supplementary Table 1). Illustrative examples of novel miRNAs are depicted in Fig. 1. One of the depicted novel miRNAs is located in a conserved genomic region and is derived from an intron of the Gnb2 gene. The other example maps to a genomic region with no features annotated in the Ensembl and RefSeq Gene tracks of the UCSC genome browser. The targets (tissues and mRNAs) and functions of these novel miRNAs remain to be discovered.

Both age and CR alter the abundance of circulating known miRNAs

To determine potential effects of age and/or CR on the circulating levels of known miRNAs, the expression values generated by miRDeep2 were analyzed with the Bioconductor package edgeR [48]. The samples were first examined with plotMDS, an edgeR function that produces a multi-dimensional scaling plot in which distances reflect the biological coefficient of variation between samples. One dimension of the plotMDS adequately separated the young group from both old groups (old control and old CR), while the other dimension adequately separated the old control group from the old CR group (Fig. 2). This analysis confirms the homogeneity of the replicates, and reveals distinct effects of age and CR on the abundance of circulating miRNAs.

We performed pairwise comparisons between the young and old control groups to measure differences in the circulating miRNAs associated with old age, and pairwise comparisons between the old control and old CR groups, to unravel any potential effect of CR on the age-associated changes in circulating levels of miRNAs. Differences were considered significant if the miRNAs achieved a minimum of 10 counts per million (cpm) reads in at least one of the 3 experimental groups (young, old control, and old CR), the fold change between any two groups was ≥ 1.5, and the p-value of this difference was < 0.05. Application of these criteria revealed that aging increased the circulating levels of 45 known miRNAs, but decreased the circulating levels of only 3 known miRNAs. CR either completely or partially mitigated these age-associated changes. A subset of miRNAs with the largest age-associate fold increases is presented in Fig. 3, while a complete list of all circulating miRNAs affected by both aging and CR are reported in Table 2. Aging was associated with the increase or decrease in the levels of 28 and 44 circulating miRNAs, respectively, without CR having any significant effects on these age-associated changes (Supplementary Table 2). On the other hand, CR altered the circulating levels of a group of miRNAs that were not affected by age; CR decreased the levels of 18, and increased levels of 2 circulating miRNAs (Supplementary Table 3).

The functional significance of the observed age- and/or CR-associated changes remains to be established: they would be biologically meaningful only if circulating miRNAs are taken up by peripheral tissues, and retain functional mRNA targeting capabilities that regulate gene expression in recipient tissues. Pertinent to this point, circulating miRNAs can bind high-density lipoprotein (HDL), and are taken up by recipient cells where they directly target mRNAs [12]; there are significant differences in the circulating HDL-miRNA profile between normal subjects and patients with familial hypercholesterolemia, and delivery of atherosclerotic HDL-miRNAs to hepatocytes alters the expression of genes related to lipid metabolism, inflammation, and atherosclerosis [12]. In another study, a set of brain miRNAs (miR-137, -181c, -9, -29a/b) was reported to circulate at decreased levels in the serum of Alzheimer's Disease (AD) patients and AD animal models, and these miRNAs are potentially involved in AD through the regulation of ceramides [49]. This evidence for function of circulating miRNAs raises the possibility that they may have utility as diseases markers and even therapeutic targets.

Pathways relevant to aging are associated with circulating miRNAs whose serum levels are increased with age and decreased by CR

The evidence that circulating miRNAs may be delivered to and directly regulate gene expression in peripheral tissues prompted us to carry out a functional characterization of circulating miRNAs, to assess their potential impact on gene expression in recipient tissues. We focused on the group that shows an age-associated increase in serum levels that is antagonized by CR. First, we identified the cellular mRNAs that may be targeted by miRNAs in the group. Cellular mRNAs were identified with miRDB [50], an algorithm that uses machine learning to predict miRNA targets. Only those mRNAs with a prediction score > 80 were considered, since they are most likely to be real targets [50]. Second, we used DAVID and PANTHER [51, 52] to perform a functional annotation of the potential mRNA targets. The pathways identified by this strategy (Table 3; the pathways in Table 3 are the only pathways with significant scores and p values) reveal that the circulating miRNAs of interest may be involved in the regulation of two GO biological processes (‘positive regulation of macromolecule biosynthetic process’ and ‘negative regulation of apoptosis’), and one PANTHER pathway (‘Wnt signaling’).

The GO term ‘Positive regulation of macromolecule biosynthetic process’ represents the anabolic pathways that use ATP to synthesize the four classes of macromolecules needed by the cell: polysaccharides, lipids, nucleic acids, and proteins. Both aging and CR are potent modulators of metabolism. Aging decreases macromolecular turnover which may underlie the age-related accumulation of oxidative damage, while CR is thought to extend lifespan by reducing metabolic rate and lowering the production of toxic by-products of metabolism [29]. Aging decreases the expression and activity of enzymes required to mobilize proteins for the production of metabolic energy, while CR increases the catabolism of protein and lipid derived from proteolysis and autophagy to generate substrates for energy generation. Moreover, gene expression profiling studies suggested that CR may retard aging by causing a metabolic shift toward increased protein turnover and decreased macromolecular damage [24-26]. In addition, the most widely accepted theory explaining the lifespan-extending effects of CR is that it shifts energy usage away from growth and reproduction and toward maintenance and stress resistance during times of nutritional stress [53]. Thus the miRNAs whose circulating levels are increased with age, and decreased by CR, are associated with functions that are closely associated with known changes in cellular metabolism that occur with aging. The causal relationship between these miRNAs and manifestations of aging remains to be determined.

The second GO term obtained by the functional analysis of age- and CR-regulated circulating miRNAs is ‘negative regulation of apoptosis’. Aging increases apoptosis in postmitotic tissues, including brain, skeletal and cardiac muscle, and germ cells, impeding the homeostasis of somatic organs and stem cell self-renewal [54-56]. In contrast, the suppression of apoptosis during tumorigenesis may underlie the age-associated increased prevalence of cancers [57]. CR increases apoptosis and decreases cellular proliferation in mitotic tissues, where it selectively eliminates preneoplastic and neoplastic cells, which are more sensitive to apoptosis than normal cells (reviewed in [29]). Three months of CR decrease the number and volume of chemically induced preneoplastic foci by 85%, suggesting that CR may exert its anticarcinogenic activity by preferentially inducing apoptosis in tumors and preneoplastic foci [58, 59]. Given these complex effects of aging and CR on apoptosis, our findings of age- and CR-associated changes in the circulating levels of miRNAs may reveal another layer of apoptosis regulation in potential recipient tissues during aging and in response to CR.

Wnt signaling regulates cell proliferation and differentiation, apoptosis, and stem cell renewal, and plays a complex role in aging (reviewed in [60]). The outcome of its effects on aging depends on whether ß-catenin associates with the transcription factor FOXO to induce senescence, or with TCF/LEF to stimulate stem cell renewal. Since FOXO transcription factors regulate the rate of aging and may mediate the antineoplastic effects of CR [61, 62], it is tempting to speculate that the CR-associated changes in the circulating levels of miRNAs may contribute to fine tuning of Wnt signaling in peripheral tissues, to preclude the age-induced senescence of stem and proliferating cells and delay the onset of age-related disorders.

Old age diseases are associated with circulating miRNAs whose serum levels are increased with age and decreased by CR

miRNAs regulate a wide range of pathological processes, including age-associated impairments such as neoplasia, inflammation and neurodegeneration [63-66]. There is by now a sizable body of data linking dysregulation of miRNA expression to human disease. Databases [67-69] containing detailed information on microRNA-disease relationships provide a resource for study of the potentially pathogenic role of microRNA dysregulation. We searched the miR2Disease and HMDD databases [67, 68] for associations between specific diseases and the circulating miRNAs whose serum levels we have found to be increased with age and decreased by CR. As shown in Supplementary Table 4, by far the largest number of diseases associated with these circulating miRNAs are cancers; other associated diseases include neurodegenerative, cardiovascular, and inflammatory disorders. All of these pathologies are linked to old age. This association implies that circulating miRNAs whose serum levels are increased with age and decreased by CR may participate in the pathogenesis of age-induced diseases, and that their modulation by CR may underlie the anti-aging effects of CR.

Mice and diets

One-month-old male mice of the long-lived B6C3F1 strain were purchased from Harlan (Indianapolis, IN). One week after arrival, mice were individually housed and randomly assigned to one of two groups, control or calorie restricted (CR). Control mice were fed 93 kcal/wk of a defined control diet (AIN-93M, diet no. F05312, BIO-SERV). CR mice were fed 52.2 kcal/wk of a defined CR diet (AIN-93M, diet no. F05314, BIO-SERV). The CR mice consumed <40% fewer calories than the control group, but the CR diet was enriched so that CR mice consumed approximately the same amount of protein, vitamins, and minerals per gram of body weight as control mice. All mice had free access to water. Mice were maintained at 20-24°C and 50-60% humidity with lights on from 0600 to 1800 h. Sentinel mice were kept in the same room as the experimental mice, and serum samples were screened every 6 months for titers against 11 common pathogens. No positive titers were found during these studies. At 27-months of age, mice were euthanized, and blood was collected through cardiac puncture and processed immediately. A group of control mice were euthanized at 7 months of age and used as a young control group. Each group consisted of 3 mice. The Institutional Animal Care and Use Committee of the University of California, Riverside, approved animal protocols.

Serum collection, RNA isolation, and small RNA library construction

Immediately after collection, blood was transferred to BD Microtainer tubes (Becton, Dickinson and Company), incubated for 30 min at room temperature to allow blood clotting, and then centrifuged at 5,000 g for 10 min. The serum supernatant was transferred to new tubes, centrifuged at 16,000 g for 15 min to remove any residual cells and cell debris, and stored at −80 °C before use. Isolation of total RNA including small RNA was performed with miRNeasy kit (Qiagen) according to the manufacturer's protocol with the exceptions of mixing 2 mL of Qiazol reagent with 0.4 mL serum, loading the entire aqueous phase onto a single column from the MinElute Cleanup Kit (Qiagen), and eluting the RNA in 20 μL of RNase-free water.

One fourth (5 μL) of the RNA isolated from each serum sample was used to construct sequencing libraries with the Illumina TruSeq Small RNA Sample Prep Kit, following the manufacturer's protocol. Briefly, 3′ and 5′ adapters were sequentially ligated to small RNA molecules and the obtained ligation products were subjected to a reverse transcription reaction to create single stranded cDNA. To selectively enrich those fragments that have adapter molecules on both ends, the cDNA was amplified with 15 PCR cycles using a common primer and a primer containing an index tag; this allows multiplexing and sequencing of different samples in a single lane of a flowcell. The amplified cDNA constructs were gel purified, and validated by checking the size, purity, and concentration of the amplicons on the Agilent Bioanalyzer High Sensitivity DNA chip. The libraries were pooled in equimolar amounts, and sequenced on an Illumina HiSeq 2000 instrument to generate 50 base reads. Image deconvolution and quality values calculation were performed using the modules of the Illumina pipeline.

MiRDeep2 analysis of sequencing reads

Raw sequencing reads were analyzed with miRDeep2 [47]. Briefly, miRDeep2 pre-processed raw sequencing reads by removing the 3' adapter sequence and discarding reads shorter than 18 nucleotides, before aligning reads to the mouse genome (NCBI37/mm10). For the purpose of analyzing the sequenced miRNAs, the known miRNA input was from miRBase v.19, and Rattus norvegicus was designated as the related species. MiRDeep2 estimates expression levels of known miRNAs, and identifies novel miRNAs. Small RNAs identified by short read sequencing are derived from longer RNAs, and may or may not be true miRNAs. The miRDeep2 algorithm is a probabilistic algorithm based on the miRNA biogenesis model, and is designed to detect miRNAs from deep sequencing reads. It aligns reads to potential hairpin structures in a manner consistent with Dicer processing of hairpins to produce mature miRNAs, and assigns scores that measure the probability that hairpins are true miRNA precursors. A small RNA derived from the 5' end of a predicted precursor is considered to be a miRNA if it is highly abundant relative to small RNAs derived from loop and star regions of the precursor, and less likely to be a true miRNA if it is present in similar proportions to the loop and star. The miRDeep2 algorithm uses this principle to produce a log-odds score that a small RNA is a true miRNA; it outputs a scored list of known and novel miRNAs as well as their expression levels.

Statistical analysis of differential miRNA expression

In addition to identifying mature miRNAs in deep sequenced small RNA samples, miRDeep2 also generates expression values for the detected miRNAs. To test for differential abundance of circulating miRNAs between the experimental groups (young control, old control, and old CR), expression data for known miRNAs produced by miRDeep2 were used as input for the Bioconductor package edgeR [48] to quantify the differential abundance of circulating miRNAs. The algorithm of edgeR fits a negative binomial model to the count data, estimates dispersion, and measures differences using the generalized linear model likelihood ratio test which is recommended for experiments with multiple factors, such as the simultaneous analysis of the effects of age and diet in our study. The fitted count data was analyzed by performing pairwise comparisons between the different experimental groups. Young and old control groups were compared to measure the differential abundance in circulating miRNAs associated with old age. Old control and old CR groups were compared to determine whether CR has an effect on any age-associated changes. The differentially abundant miRNAs identified by edgeR were further filtered to remove miRNAs with less than 10 counts per million (cpm) in at least one of the 3 experimental groups.

Prediction and functional annotation of peripheral genes potentially targeted by differentially abundant circulating miRNAs

mRNA transcripts that are potentially targeted by circulating miRNAs were extracted from miRDB database, which uses machine learning to predict miRNA targets; we included only mRNAs with prediction score > 80, since they are most likely to be real [50]. The predicted mRNAs were functionally annotated with DAVID and PANTHER [51, 52]. PANTHER uses pathways compiled by experts to find pathways that are overrepresented in a list of genes by comparing it to a reference list, e.g. all genes in the mouse genome as used here. A Bonferroni correction for multiple testing was applied to the analysis. We also used DAVID to interrogate representation in Gene Ontology (GO). The DAVID algorithm measures the similarities among GO terms based on the extent of their associated genes and assembles the similar and redundant GO terms into annotation clusters. Each GO term in a cluster is assigned a Fisher Exact p-value representing the degree of enrichment of the GO term in the input gene list. Each cluster is assigned an enrichment score to rank its biological significance. Thus, a biologically significant cluster (high enrichment score) is generated only when most of its GO term members have significant enrichment values (low Fisher Exact p-values). The resulting clusters were further curated to keep only GO terms with p-values < 0.05.