MicroRNAs (miRNAs) are small, non-coding RNAs that regulate gene expression and play a critical role in development, homeostasis, and disease. Despite their demonstrated roles in age-associated pathologies, little is known about the role of miRNAs in human aging and longevity.
Results
We employed massively parallel sequencing technology to identify miRNAs expressed in B-cells from Ashkenazi Jewish centenarians, i.e., those living to a hundred and a human model of exceptional longevity, and younger controls without a family history of longevity. With data from 26.7 million reads comprising 9.4 × 108 bp from 3 centenarian and 3 control individuals, we discovered a total of 276 known miRNAs and 8 unknown miRNAs ranging several orders of magnitude in expression levels, a typical characteristics of saturated miRNA-sequencing. A total of 22 miRNAs were found to be significantly upregulated, with only 2 miRNAs downregulated, in centenarians as compared to controls. Gene Ontology analysis of the predicted and validated targets of the 24 differentially expressed miRNAs indicated enrichment of functional pathways involved in cell metabolism, cell cycle, cell signaling, and cell differentiation. A cross sectional expression analysis of the differentially expressed miRNAs in B-cells from Ashkenazi Jewish individuals between the 50th and 100th years of age indicated that expression levels of miR-363* declined significantly with age. Centenarians, however, maintained the youthful expression level. This result suggests that miR-363* may be a candidate longevity-associated miRNA.
Conclusion
Our comprehensive miRNA data provide a resource for further studies to identify genetic pathways associated with aging and longevity in humans.
MicroRNAs (miRNAs) are small non-coding RNAs that regulate gene expression at the post-transcriptional level[1]. Mature miRNAs, between 18–25 bp in length, are transcribed as primary-miRNA (pri-miRNA) molecules which contain a characteristic stem loop structure. This stem loop targets pri-miRNA for processing by a number of RNAses, namely Drosha and Dicer, which produce a short RNA duplex[1]. From the duplex one or both strands is incorporated into the RNA inducing silencing complex (RISC), resulting in an active miRNA. Active miRNAs target the 3′ UTR of a mRNA based on sequence homology[2]. The nucleotides in the 2–7 position of the 5′ end of the mature miRNA comprise a “seed region”. Once an mRNA is targeted by a miRNA it can be degraded or its translation can be repressed through conserved mechanisms leading to downregulation of gene expression[1, 3].
First described in C. elegans[3], many miRNAs have been discovered across a wide range of species, including humans. There is increasing evidence that miRNAs are critical in a number of essential biological processes, including cell differentiation[4], immune response[5], cancer[6, 7], and life span[8]. Thus far, 1048 human miRNA sequences have been identified through cloning, sequencing, or computational analysis (mirBase, release 16, 2010)[9–11].
The multitude of important roles played by miRNAs suggests that they are a critical component of gene regulatory networks. However, the quantification of miRNAs has been technically challenging due to their small size, low copy number, interference from other small RNAs, and contamination by degradation products of mRNAs or other RNA species. Until recently only known and computationally predicted miRNAs have been interrogated using hybridization-based array methods, which suffer from cross-hybridization, and inability to discover novel miRNAs due to limited array content. The increased availability and affordability of massively parallel sequencing now offer an opportunity to gain a high-resolution view of miRNA expression, overcoming past experimental limitations[12]. MiRNA-seq has been utilized to discover novel and quantify expression levels of miRNAs in several species, including humans[13, 14].
Expression levels of genes are heritable in humans as quantitative phenotypes, measurable in a variety of cell types, including B-cells[15]. Recent studies, including our own[16, 17], have demonstrated that B cells reflect functional characteristics of the donor and can be a useful tool for studying genotype-driven molecular endpoints such as gene expression, and expression quantitative trait locus (eQTL) analysis[18, 19]. B cells can act as surrogate tissues whenever there is correlation between the expression levels of B cells and phenotypes of interest[20, 21] and a large number of eQTLs originally identified in B cells can also be detected in multiple primary tissues[22, 23]. Thus B-cells have been increasingly utilized for expression quantitative trait loci (eQTLs) studies[19, 24] and as a cell model to assess gene expression responses[25]. While comprehensive mRNA expression data for human B-cells, obtained by RNA-seq, are available for this purpose[26], to date no such analysis was performed to identify miRNAs expressed in B-cells. Here, we present comprehensive miRNA transcriptome profiles of B-cells from Ashkenazi Jewish centenarians and younger control individuals by miRNA-seq, providing a resource that could serve as a basis for establishing gene regulatory interactions between miRNAs and their target mRNAs in human B-cells.
Results
Discovery of miRNAs expressed in B-cells of centenarians and controls
We generated small RNA libraries of B-cells from 3 Ashkenazi Jewish female centenarians (mean age 104) and 3 Ashkenazi Jewish female younger controls (mean age 63). Sequencing of these libraries by an Illumina technology platform yielded a total of 12.9 × 106 reads from centenarians and 13.8 × 106 reads from controls (Additional file1: Table S1). After removal of low quality reads and redundancy, we had a total of 1.1 × 106 and 1.0 × 106 unique reads for the centenarians and the controls, respectively. Of the unique reads, those with a read number greater than 10 in more than 50% (n = 3) of the individuals sequenced (n = 6) were aligned to a database of known miRNAs (http://www.mirbase.org/) and other known small RNA species. MiRNAs with less than 10 reads were not included due to the error rate of Illumina sequencing and stochastic variation in gene expression[12]. For all comparisons, the number of reads for a given miRNA was normalized by division of the total number of miRNA reads in that library, yielding a normalized read count for each miRNA[12]. We identified a total of 276 known miRNAs (Additional file2: Table S2) with a wide dynamic range of read counts ranging from 10 to over 1 million (Figure1A). About 11% of these miRNAs had a copy number greater than 10,000 and about 67% a copy number less than 1,000. The identification of 276 miRNAs out of the 1048 reported human miRNAs (mirBase, release 16, 2010)[9–11] is in the similar range as observed by others for a particular cell type[14, 27], suggesting a saturation of miRNA-sequencing. The normalized read count for each miRNA indicated that let-7f was the most abundant miRNA detected in B-cells (Figure1B). The top 10 miRNAs by read count corresponded to ~80% of all sequence reads in the B-cells (Figure1B). A complete list of miRNAs sequenced in each library, normalized read counts, and fold differences between centenarian and control cells is provided in Additional file2: Table S2.
Figure 1
Known miRNAs discovered in B-cells from centenarians and controls. (A) Distribution of the miRNAs expressed in B-cells according to their sequence counts in centenarians and controls. On the vertical axis is a bin for miRNAs with sequence read counts within the given range. In the horizontal axis is the numbers of miRNAs found for each of the magnitudes of expression. (B) Top 10 known miRNAs expressed in B-cells from centenarians and controls, which constitute nearly 80% of all miRNAs expressed in the B-cells.
To discover previously unknown novel miRNAs among our short sequence RNAs that did not align to the known miRNA database, we took advantage of the fact that miRNAs are sequentially processed to their mature form from stem-loop containing pre-miRNAs. Utilizing the mirDeep program[28] we identified small RNAs that correspond to processing of characteristic stem loop containing pre-miRNAs. MirDeep uses a number of criteria to select novel miRNAs[29], including a minimum free energy of less than −25 kcal/mol, 25% sequence homology with known miRNA precursors, sufficient GC content, and no alignment to too many different genomic locations. We discovered a total of 8 potentially novel miRNAs with a characteristic stem loop structure that meet these criteria (Table1). None of the novel miRNAs were found to have expression over 100,000 reads in any of our samples. All putative miRNA precursors, with the exception of sRNA8, had a lowest free energy and best folding structure resembling known miRNA (Table1). All of the putative miRNAs and their genomic locations have been submitted to mirBase.
Table 1
Novel miRNAs discovered in B-cells, their chromosome location, mature sequence, and putative folding structure
Novel miRNA were predicted using the mirDeep package. The pre-miRNA sequence of sRNA1 and sRNA5 were present in two genomic locations where miRNA have not been previously predicted. The fold column was generated from mFold[30].
IsomiRs and variability of miRNA processing
Mature miRNAs are not always produced with the same sequence even when originating from the same precursor miRNA. Such mature variants have been termed “isomiRs”[13]. A tissue-specific expression pattern of isomiRs during development has been demonstrated in Drosophila, suggesting a specific role of isomiRs in gene regulatory networks[31]. A failure to identify all isomiRs can alter sequence counts and lead to incorrect conclusions of miRNA abundance[13]. Most of the variation in isomiRs can be described by the variability in Dicer 1 and Drosha processing of the miRNAs[32]. We compared the composition of isomiRs between B-cells originating from centenarians vs. controls and identified a total of 694 variants from 261 miRNAs. Similar isomiR patterns were identified in centenarian and control cells, suggesting that global miRNA biogenesis is not significantly different between centenarians and controls. A subset of identified isomiRs and their compositions within the groups is listed in Table2. We found variants within our most abundant miRNAs (e.g., hsa-let-7a and hsa-miR-21) and within our differentially expressed miRNAs (e.g., hsa-miR148a and hsa-miR-193b). There is evidence that the most abundant isomiR in a sample is not always the same as the reference sequence reported in mirBase[33]. Indeed in our B-cells a variant of hsa-miR-21 not reported in mirBase was the dominant isomiR. A complete list of sequenced isomiRs from each individual miRNA discovered in B-cells is available in Additional file3: Table S3.
Table 2
Representative IsomiRs detected in B-cells
miRNA
Sequence
% Control
% Centenarian
hsa-miR-21
TAGCTTATCAGACTGATGTTGAC
53.8%
52.0%
TAGCTTATCAGACTGATGTTGA
43.3%
44.9%
TAGCTTATCAGACTGATGTTGAA
2.7%
2.8%
TAGCTTATCAGACTGATGTTGAT
0.2%
0.3%
hsa-miR-193b
AACTGGCCCTCAAAGTCCCGCT
84.8%
88.6%
AACTGGCCCTCAAAGTCCCGCTT
8.7%
5.9%
AACTGGCCCTCAAAGTCCCGCTA
6.2%
5.2%
AACTGGCCCTCAAAGTCCCGCTG
0.3%
0.3%
hsa-let-7a
TGAGGTAGTAGGTTGTATAGTT
86.8%
86.9%
TGAGGTAGTAGGTTGTATAGTTA
9.3%
9.2%
TGAGGTAGTAGGTTGTATAGTTT
3.6%
3.5%
TGAGGTAGTAGGTTGTATAGTTG
0.4%
0.3%
hsa-miR-148a
TCAGTGCACTACAGAACTTTGT
94.3%
94.3%
TCAGTGCACTACAGAACTTTGTC
3.9%
3.6%
TCAGTGCACTACAGAACTTTGTT
1.3%
1.4%
TCAGTGCACTACAGAACTTTGTA
0.5%
0.7%
Some of the top isomiRs are shown with the percent contribution of each isomiR to the total read count within each group. Reference sequence is denoted in bold and isomer variation from the reference is noted with bold and underscore.
MiRNAs differentially expressed between centenarians and controls
To identify differentially expressed miRNAs, we performed a Fisher’s exact test and adjusted for multiple testing with Bonferroni correction[33–35]. We identified a total of 24 differentially expressed known miRNAs with a fold change > 2.0 and Bonferroni corrected p-value < 0.05 (Table3). Of the 24 differentially expressed miRNAs, 22 were upregulated and only two downregulated in centenarians as compared to controls (Figure2). Normalized read counts for the differentially expressed miRNAs ranged from the 20s to over 100,000. A complete list of miRNAs sequenced in each library, normalized read counts, and fold changes between centenarians and controls is provided in Additional file2: Table S2.
Table 3
Differentially expressed miRNAs in B-cells from centenarians vs. controls
miRNA
Up regulated in
Control counts (standard error)
Centenarian counts (standard error)
P values
Bonferonni adjusted p values
Fold change
hsa-miR-122
Centenarian
19 (6.26)
102 (6.39)
1.06E-16
2.44E-14
5.83
hsa-miR-363*
Centenarian
184 (23.30)
837 (81.32)
3.6E-112
8.2E-110
4.94
hsa-miR-345
Centenarian
59 (10.65)
227 (32.69)
9.84E-28
2.27E-25
4.18
hsa-miR-20b
Centenarian
306 (25.72)
1108 (240.38)
7.7E-122
1.8E-119
3.93
hsa-miR-454
Centenarian
39 (6.55)
126 (34.59)
1.02E-13
2.36E-11
3.51
hsa-miR-1974
Centenarian
4515 (660.27)
14179 (668.31)
0
0
3.41
hsa-miR-223*
Centenarian
580 (39.17)
1800 (204.29)
2.2E-167
5.1E-165
3.37
hsa-miR-99b
Centenarian
97 (15.32)
276 (83.78)
2.15E-24
4.96E-22
3.09
hsa-miR-181a*
Centenarian
627 (24.30)
1641 (307.63)
4.5E-123
1E-120
2.84
hsa-miR-363
Centenarian
4732 (260.52)
11971 (2286.46)
0
0
2.75
hsa-miR-21*
Centenarian
2529 (286.89)
6313 (1153.62)
0
0
2.71
hsa-miR-92b*
Centenarian
99 (11.67)
245 (42.00)
2.45E-18
5.66E-16
2.69
hsa-miR-20b*
Centenarian
319 (77.91)
708 (84.88)
5.78E-42
1.33E-39
2.41
hsa-miR-148a
Centenarian
11599 (1225.92)
25583 (655.62)
0
0
2.40
hsa-miR-1975
Centenarian
1654 (86.54)
3529 (335.65)
2.1E-188
4.9E-186
2.32
hsa-miR-502-3p
Centenarian
191 (13.48)
387 (14.66)
4.03E-20
9.31E-18
2.20
hsa-miR-181c
Centenarian
153 (6.46)
310 (63.42)
1.8E-16
4.15E-14
2.20
hsa-miR-1259
Centenarian
106 (25.78)
210 (29.26)
4.01E-11
9.27E-09
2.15
hsa-miR-148a*
Centenarian
339 (46.74)
656 (14.26)
3.38E-30
7.82E-28
2.10
hsa-miR-192
Centenarian
7126 (1026.34)
13754 (2446.01)
0
0
2.10
hsa-miR-361-5p
Centenarian
307 (20.35)
587 (89.41)
1.65E-26
3.81E-24
2.08
hsa-miR-9
Centenarian
341 (26.37)
629 (37.72)
3.8E-26
8.78E-24
2.00
hsa-miR-151-3p
Control
4515 (1288.82)
1352 (223.77)
0
0
3.07
hsa-miR-151-5p
Control
257 (72.35)
80 (13.37)
7.22E-20
1.67E-17
2.96
All miRNA determined to be differentially expressed with a Bonferroni adjusted p-value < 0.05 and fold change > 2.0.
Figure 2
Relative abundance and fold changes of miRNAs differentially expressed between centenarians and controls. Each sphere along the X-axis represents a single miRNA that was found differentially expressed between centenarians and controls. The diameter of a sphere represents the relative abundance of the miRNA read count. The Y-axis denotes the fold changes between the two groups. Positive values denote up-regulation, and negative values denote down-regulation, of the miRNAs in centenarians as compared to controls.
Cross platform comparison of differential miRNA expression
Until recently, differentially expressed miRNAs were discovered using miRNA microarray technology, which is subject to certain limitations such as difficulties in detecting low abundance miRNAs and a limited dynamic range, essentially constraining the detection of subtle fold changes. These limitations call for an additional technical validation step, for example, quantitative real time PCR (qRT-PCR), which has been considered the gold-standard of differential analysis of gene expression. We conducted qRT-PCR analysis to validate the differential expression of 7 miRNAs out of the 24 miRNAs found to be significantly different in B-cells from centenarians as compared to controls by Illumina sequencing analysis. Five miRNAs (hsa-miR-363*, hsa-miR-1974, hsa-miR-223*, hsa-miR-148a, hsa-miR-148a*) produced expression patterns consistent with the Illumina sequencing data (Figure3). However, the qRT-PCR data clearly underestimated the fold changes in gene expression when compared with the Illumina sequencing analysis. In addition, qRT-PCR failed to detect the expression of hsa-miR-122 and hsa-miR-502-3p. These results indicate that differential analysis of miRNAs by deep sequencing should be the gold standard for quantification of miRNA expression, especially since there is evidence for both 5′ and 3′ editing of mature miRNAs, which reduces the efficiency of miRNA detection by hybridization-based method such as qRT-PCR[36].
Figure 3
Cross platform comparison of miRNA differential expression. Quantitative RT-PCR with TaqMan probes was used to validate differential expression of miRNAs detected by miRNA-seq. Standard error of the mean for qPCR values are shown.
Targets of differentially expressed miRNAs
Among the 24 differentially expressed miRNAs, several have been previously characterized with known validated targets (Additional file4: Table S4), including hsa-miR-148 and hsa-miR-122[5, 37, 38]. Target sites of the 24 differentially expressed miRNAs were predicted using available computational approaches, in particular TargetScan (http://www.targetscan.org)[39, 40], miRanda (http://www.microrna.org/microrna/home.do)[41], and PicTar (http://pictar.mdc-berlin.de/)[42]. These software packages work by finding an absolute base pair homology of the miRNA seed region (bases 2–8 on the 5′ end of the miRNA) to the 3′ UTR of all mRNAs. In addition, they take into consideration evolutionary conservation of 3′ UTR bases as well as RNA accessibility to the RISC complex. Due to the large number of predicted targets we only considered target genes with 2 predicted sites, an accessible 3′UTR and consensus amongst all prediction algorithms[27, 43]. Based on this approach we identified predicted target genes for 10 out of the 24 differentially expressed miRNAs (Additional file4: Table S4). To explore the potential pathways regulated by these miRNAs, we conducted Gene Ontology analysis utilizing the “GoStat” analysis tool (http://gostat.wehi.edu.au/). We took the total list of predicted and validated targets of the differentially expressed miRNAs and determined if any Gene Ontology categories were over-represented within our list[44]. We found enrichment of functional pathways implicated in the aging process, including cell metabolism, cell cycle, cell signaling, and cell differentiation (Table4) among other significant GO categories (Additional file5: Table S5). In addition, we compared our putative target list to the GenAge database which contains a collection of genes implicated in aging[45] using a standard gene set enrichment approach[46]. We found that our target list has a significant over-representation of these genes (Table4).
Table 4
Representative GO categories enriched in targets of differentially expressed miRNAs
GO TERM
GO ID
p-Value (Bonferonni Adjusted)
cell differentiation
GO:0030154
1.52 × 10-9
cellular metabolic process
GO:0044237
3.82 × 10-6
cell-cell signaling
GO:0007267
1.05 × 10-5
cell proliferation
GO:0008283
6.66 × 10-5
cell cycle
GO:0007049
8.02 × 10-4
GenAge Database
4.5 × 10-2
Cross-sectional analysis of miRNA expression in different age groups
Since 22 miRNAs out of 24 differentially expressed miRNAs were discovered to be upregulated in centenarians (mean age: 103) as compared to controls (mean age: 63), we set out to determine the mode of differential expression over a range of age groups. If upregulation is simply age-related, expression will increase with age (Figure4A). In contrast, if the upregulation is longevity-related, patterns of youthful expression will be preserved in centenarians while its expression levels in control individuals decline with age (Figure4B). Among the differentially expressed miRNAs, we chose hsa-miR-363* for a further validation study because its upregulation in centenarians was detectable by qRT-PCR analysis with > 2 fold change (Figure3). We performed a cross-sectional study to measure miRNA expression levels in B-cells established from 4 age groups of Ashkenazi Jewish individuals; 50–60 (n = 5), 70–80 (n = 5), 80–90 (n = 5), >100 years (n = 12). We found that the expression levels of hsa-miR-363* significantly declined with age in control individuals (p < 0.05), while centenarians maintained comparable “high” expression levels of the youngest age group (Figure4C). The expression differences between centenarians and younger elderly controls were significant (p = 0.05 for 70–80 vs. > 100; p = 0.003 for 80–90 vs. > 100) (Figure4C), suggesting that hsa-miR-363* is a candidate longevity-associated miRNA.
Figure 4
Cross-sectional expression of select miRNA in different age groups. The dotted line shows the age of the discovery group (mean age 63) for differential analysis with centenarians (mean age 104). A) If differential expression is simply age-related, expression will increase with age. B) If the differential expression is longevity-related, patterns of youthful expression will be preserved in centenarians while its expression levels in control individuals decline with age. C) A cross-sectional expression analysis of hsa-miR-363* by quantitative RT-PCR showed the youthful preservation mode: its expression levels significantly decline with age in control individuals (p < 0.05), while centenarians maintain the expression levels of the youngest age group (50–60). The significance of expression differences was assessed by student’s t-test. Each age group had 5 individuals while centenarians had 12 individuals.
Discussion
We employed a massively parallel sequencing technology to identify miRNAs expressed in B-cells established from Ashkenazi Jewish Centenarians and younger elderly controls without a family history of longevity. In this study, which represents the first comprehensive studies to analyze miRNAs in human B-cells, we obtained 26.7 million reads comprising 9.4 × 108 bp from 3 centenarian and 3 control individuals. We found a total of 284 miRNAs expressed in human B-cells, eight of which were previously unknown, putative novel miRNAs.
Profiling miRNA transcriptomes has gained importance with increasing evidence for the role of miRNA expression in defining cellular phenotypes. Prior to the advent of increasingly cost-effective, small RNA sequencing, microarrays were the prevailing methods for miRNA expression analysis. However, apart from noise due to cross-hybridization to probes with similar short-sequences, array technologies are limited in several key ways. First, microarrays are limited to the probe sets, eliminating the possibility to discover novel miRNAs. Second, arrays are only able to identify the relative abundance of miRNAs but not their absolute numbers[12]. Sequencing allows accurate detection of expression levels over a wide dynamic range, including low copy number miRNAs and subtle fold changes between test and control groups. Finally, array-based methods cannot differentiate isomiRs (mature miRNAs originated from the same precursor but differing in one or more bases) resulting from the high variability in miRNA processing[13, 31]. Using miRNA-seq of human B-cells, we were able to identify novel miRNAs and isomiRs, and detected expression of low-copy miRNAs and subtle fold changes between centenarians and controls, which could not be measured using TaqMan qPCR.
Genome-wide transcriptional profiles of miRNAs in C. elegans[47] and human blood cells[48] showed that the majority of miRNAs decreased in expression with age. Likewise, cellular senescence was accompanied by significant down-regulation in miRNA expression[49]. Our results suggest that centenarians may be more resistant to such age-related down-regulation of miRNAs because the overwhelming majority of significant miRNAs are increased in centenarians as compared to controls in B-cells: we found 22 miRNAs significantly upregulated, with only 2 downregulated, in centenarians as compared to younger controls. Strikingly, expression levels of 5 miRNAs shown to be increased in centenarians (hsa-miR-148a, hsa-miR-345, hsa-miR-361-5p, hsa-miR-192, hsa-miR-454) have been demonstrated to be down-regulated during cellular senescence (Table5). Among these, hsa-miR-148a was also found to be down-regulated with age in peripheral blood mononuclear cells (PBMCs)[48]. Together, these results suggest that miRNAs associated with longevity can be detected by cross-sectional expression analysis as utilized in this work.
Table 5
Relevance of miRNAs differentially expressed in centenarians to senescence and aging
By conducting a cross-sectional analysis of miRNA expression in different age groups, we found that expression of hsa-miR-363* significantly declined with age in control individuals, whereas centenarians maintain comparably high expression levels, similar to the observed in the middle age group. The results suggest that hsa-miR-363* is a candidate longevity-associated miRNA. Previously miRNA* sequences were thought to be degraded. However a growing body of work challenges the dogma that miRNA* is simply a non-functional byproduct of miRNA biogenesis, suggesting instead that miRNA* plays a significant role in cellular function and human disease[51]. Furthermore, they are implicated in the aging of C. elegans[52]. The predicted targets of hsa-miR-363* include PTEN, BCL2, AKT1, and IGFBP5 among genes listed in the GenAge database[45]. To our knowledge this is the first cross-sectional analysis of miRNA in a human longevity cohort, identifying a potential longevity-associated miRNA.
Targets for our differentially expressed miRNAs were obtained through the literature and by bioinformatics means (Additional file4: Table S4). The predicted and validated targets of the 24 differentially expressed miRNAs are enriched in biological pathways implicated in the aging process, including cell differentiation, cell cycle, cell signaling, and cell metabolism (Table4), as well as the GenAge database[45]. Of particular interest were the genes involved in cognitive aging (ADAMTS5, APPL1, IGF2, LRP1b), aging-related signaling/transcription (E2F2, FGF7, MLL, MYB , KLF3, PTEN, RICTOR, SMAD2, SOCS5, TCF4), cell cycle/senescence/apoptosis (BCL2L11, CBL , CDC2L6, E2F2, EEF2K), and overlap with the GenAge database of known aging related genes (SHC1, E2F1, STAT3, IGF1, IRS2, PTEN, RPA1, BCL2, SP1, HBP1, MXD1). We searched for potential target genes that were predicted to bind multiple miRNAs found to be differentially expressed in our study because it has been shown that miRNAs can act in an additive manner and the co-targeting by multiple miRNAs would greatly reduce mRNA expression even if the miRNAs were in low copy number[53]. We found that RICTOR (rapamycin-insenstive companion of mTOR) was predicted to bind at least two miRNAs upregulated in centenarians as compared to controls. The 3′ UTR of RICTOR and the binding sites for miR-148a and miR-155 (fold change > 1.5) are shown in Figure5. RICTOR is part of the TOR family of genes which are integral to growth and proliferation, and down-regulation of this pathway is shown to extend lifespan[54].
Figure 5
3′UTR of Rictor showing multiple miRNA target locations.
Longevity is known to have a genetic component in humans. While the heritability of average life expectancy has been estimated to be only ~25%[55, 56], studies of centenarians indicate much stronger heritability at old age. For example, siblings of centenarians have a 4 times greater probability of surviving to age 90 than siblings of people with average life span[57]. Achieving a lifespan of 100 years is 17 and 8 times more likely for male or female siblings of centenarians, respectively, compared to their birth cohort[58]. These findings firmly established the utility of human centenarians as a model system to unravel the genetics of longevity. For this study, we used our unique cohort of centenarians and elderly controls without family history of exceptional longevity, all of genetically homogeneous Ashkenazi Jewish descent. This cohort has been used to successfully discover longevity associated genotypes and phenotypes in the past[16, 59].
Recently, miRNAs have emerged as critical regulators of gene expression and a link between multiple miRNAs and longevity[8, 60] and aging[47] has been demonstrated in C. elegans, implicating their role in regulation of lifespan and in the aging process. Since a significant number of miRNAs are evolutionarily conserved[61, 62], it is conceivable that miRNAs play a role in human longevity as well. Our finding that most of the differentially expressed miRNAs were upregulated in centenarians could point towards increased resilience of centenarians against an age-related decline in gene regulatory control. By conducting a cross-sectional expression analysis, we found a candidate longevity-associated miRNA, hsa-miR-363* (Figure4). Hsa-miR-363* shows reduced expression in control individuals with advancing age while maintaining relatively high expression levels in centenarians. Maintenance of youthful expression patterns may be beneficial and longevity-associated miRNAs may confer robustness to gene expression networks, protecting them against age-related deterioration.
In the field of human genetics, most of the surveys of gene expression have been conducted in B-cells because they are readily available and can be used multiple times under controlled experimental conditions[15, 63]. B cells act as surrogate tissues whenever there is a correlation between expression levels in these cells and phenotypes of interest[20, 21]. However, caution needs to be taken in interpreting the results, especially with negative data, as truly tissue-specific genes will not be detected in B-cells. The lack of a correlation can never be used to infer that the miRNA/gene is not involved in human longevity, and only positive results should be interpreted as in most large-scale discovery-based science.
Conclusions
Massively parallel sequencing technology allowed us to accurately detect miRNAs expressed in B-cells. Considering the increasing use of B-cells for genetic and functional studies[16, 64] our data provides a resource for designing gene expression studies and to study gene regulatory networks mediated by miRNAs. Furthermore, our results from B-cells established from a human longevity cohort may generate an opportunity to explore the possible role of miRNAs in human aging and longevity and to identify genes and pathways that are targets for age-related alteration.
Methods
Population and sample collection
All individuals are enrolled in the Longevity Genes Project, and were recruited as described previously[65]. Informed written consent was obtained in accordance with the policy of the Committee on Clinical Investigations of the Albert Einstein College of Medicine. All blood samples were processed at the General Clinical Research Center at the Albert Einstein College of Medicine in order to produce EBV transformed immortalized B-cells as a source of RNA. Total RNA was extracted from immortalized B-cell lines established from centenarians and controls using TRIZOL reagent as recommended by the manufacturer (Invitrogen).
Ethics
Experimental research reported in this manuscript has been performed with the approval of the Committee on Clinical Investigations of the Albert Einstein College of Medicine. Research carried out in this manuscript is in compliance with the Helsinki Declaration (http://www.wma.net/e/policy/b3.htm).
MiRNA-sequencing
MiRNA-seq was performed as recommended by the manufacturer (Illumina small RNA prep kit v 1.5). 10 μg of total RNA from each sample was resolved on a 15% TBE-Urea polyacrylamide gel followed by the excision of gels corresponding to the 17–35 nucleotides. Small RNAs were isolated from the gel in 300 μl of 0.3 M NaCl for 4 hours at room temperature. The small RNAs were ligated with a biotinylated RNA-DNA 3′-adaptor, gel-purified, and ligated with a 5′-adaptor. Products with both adaptors were gel-purified, reverse-transcribed, and PCR amplified for 14 cycles. Sequencing was performed on an Illumina GA1 analyzer.
Data analysis
The sequencing data was provided from the GA1 sequencer in a standard fastq format[66]. The fastq files were trimmed of adapter sequences and of low quality reads (reads which had more than 3 base-calls below sufficient quality value), through a c++ program. These sequences were then collapsed to remove redundancy using the Galaxy Genome Browser tool fastx[67]. At this point sequences from each of the samples was aligned to the known human miRNA/small RNA database or put into the mirDeep pipeline for the discovery of novel miRNAs[28]. To normalize the samples we determined how frequently the miRNA was annotated per the number of reads reported from the sequencer[12]. The basis of this idea is that for a given number of small sequences isolated from a cell there should be on average the same number of total miRNAs from the sample. Read counts from different libraries were normalized to the total reads in each sample[12].
Statistical analysis
After normalization any reads observed in fewer than 3 samples or with a copy number less than 10 were not considered for analysis; this correction removed extremely low abundance miRNAs. To identify differentially expressed miRNAs, the data was analyzed through Fisher’s exact test using a Bonferonni correction[33–35] for multiple hypothesis testing. Those miRNAs meeting a corrected cutoff with a p-value below 0.05 and with a fold change greater than 2.0 were considered differentially expressed.
Quantitative RT- PCR analysis
Total RNA was isolated from B-cells using RNA isolation kit (Qiagen, Valencia, CA) and then converted to complementary DNA using TaqMan Reverse Transcription kit (Applied Biosystems, Foster City, CA) with microRNA specific RT primer (Applied Biosystems). A TaqMan® microRNA assay was performed using AB StepOneTM real-time PCR system to quantify relative miRNAs expression in these samples. The 20 μl total volume final reaction mixture consisted of 1 μl of TaqMan microRNA specific primer, 10 μl of 2x Universal Master Mix with no AmpErase® UNG (Applied Biosystems) and 1.3 μl of complementary DNA. PCR was performed using the following conditions: 50°C for 2 min, 95°C for 10 min, 40 cycles of 95°C for 15 sec, and 60°C for 1 min. All reactions were run in duplicate to reduce confounding variance. U6 snRNA (Applied Biosystems) was used as an internal control. Means from different conditions were compared using the Student’s t-test. A significance threshold of P < 0.05 was used.
Accession numbers
The sequence data from this study have been submitted to the submitted to the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) under accession no. GSE32493 All novel microRNA sequences have been submitted to miRBase (http://www.mirbase.org).
Declarations
Acknowledgements
This work has been supported by grants to Dr. Y. Suh from the Glenn Award for Research in Biological Mechanisms of Aging, from the US National Institute of Health (RO1 AG024391, PO1 AG027734, and PO1 AG17242), and from a pilot grant from the Einstein Cancer Center.
Authors’ Affiliations
(1)
Department of Systems and Computational Biology, Albert Einstein College of Medicine
(2)
Department of Genetics, Albert Einstein College of Medicine
(3)
Department of Medicine, Albert Einstein College of Medicine
(4)
Institute for Aging Research, Diabetes Research and Training Center, Albert Einstein College of Medicine
(5)
Department of Human Population Genetics, Peking University
(6)
Department of Genetics, Erasmus Medical Center
(7)
Buck Institute for Research on Aging
(8)
Institute of Aging Research, Guangdong Medical College
Filipowicz W, Bhattacharyya SN, Sonenberg N: Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight?. Nat Rev Genet. 2008, 9 (2): 102-114.View ArticlePubMedGoogle Scholar
Du T, Zamore PD: microPrimer: the biogenesis and function of microRNA. Development. 2005, 132 (21): 4645-4652.View ArticlePubMedGoogle Scholar
Tay YM, Tam WL, Ang YS, Gaughwin PM, Yang H, Wang W, Liu R, George J, Ng HH, Perera RJ, et al: MicroRNA-134 modulates the differentiation of mouse embryonic stem cells, where it causes post-transcriptional attenuation of Nanog and LRH1. Stem Cells. 2008, 26 (1): 17-29.View ArticlePubMedGoogle Scholar
Tili E, Michaille JJ, Cimino A, Costinean S, Dumitru CD, Adair B, Fabbri M, Alder H, Liu CG, Calin GA, et al: Modulation of miR-155 and miR-125b levels following lipopolysaccharide/TNF-alpha stimulation and their possible roles in regulating the response to endotoxin shock. J Immunol. 2007, 179 (8): 5082-5089.View ArticlePubMedGoogle Scholar
Huang Q, Gumireddy K, Schrier M, le Sage C, Nagel R, Nair S, Egan DA, Li A, Huang G, Klein-Szanto AJ, et al: The microRNAs miR-373 and miR-520c promote tumour invasion and metastasis. Nat Cell Biol. 2008, 10 (2): 202-210.View ArticlePubMedGoogle Scholar
Silber J, Lim DA, Petritsch C, Persson AI, Maunakea AK, Yu M, Vandenberg SR, Ginzinger DG, James CD, Costello JF, et al: miR-124 and miR-137 inhibit proliferation of glioblastoma multiforme cells and induce differentiation of brain tumor stem cells. BMC Med. 2008, 6: 14-PubMed CentralView ArticlePubMedGoogle Scholar
Boehm M, Slack F: A developmental timing microRNA and its target regulate life span in C. elegans. Science. 2005, 310 (5756): 1954-1957.View ArticlePubMedGoogle Scholar
Morin RD, O'Connor MD, Griffith M, Kuchenbauer F, Delaney A, Prabhu AL, Zhao Y, McDonald H, Zeng T, Hirst M, et al: Application of massively parallel sequencing to microRNA profiling and discovery in human embryonic stem cells. Genome Res. 2008, 18 (4): 610-621.PubMed CentralView ArticlePubMedGoogle Scholar
Bar M, Wyman SK, Fritz BR, Qi J, Garg KS, Parkin RK, Kroh EM, Bendoraite A, Mitchell PS, Nelson AM, et al: MicroRNA discovery and profiling in human embryonic stem cells by deep sequencing of small RNA libraries. Stem Cells. 2008, 26 (10): 2496-2505.PubMed CentralView ArticlePubMedGoogle Scholar
Cheung VG, Spielman RS: Genetics of human gene expression: mapping DNA variants that influence gene expression. Nat Rev Genet. 2009, 10 (9): 595-604.PubMed CentralView ArticlePubMedGoogle Scholar
Suh Y, Atzmon G, Cho MO, Hwang D, Liu B, Leahy DJ, Barzilai N, Cohen P: Functionally significant insulin-like growth factor I receptor mutations in centenarians. Proc Natl Acad Sci U S A. 2008, 105 (9): 3438-3442.PubMed CentralView ArticlePubMedGoogle Scholar
Tazearslan C, Huang J, Barzilai N, Suh Y: Impaired IGF1R signaling in cells expressing longevity-associated human IGF1R alleles. Aging Cell. 2011, 10 (3): 551-554.PubMed CentralView ArticlePubMedGoogle Scholar
Cookson W, Liang L, Abecasis G, Moffatt M, Lathrop M: Mapping complex disease traits with global gene expression. Nat Rev Genet. 2009, 10 (3): 184-194.PubMed CentralView ArticlePubMedGoogle Scholar
Montgomery SB, Dermitzakis ET: From expression QTLs to personalized transcriptomics. Nat Rev Genet. 2011, 12 (4): 277-282.View ArticlePubMedGoogle Scholar
Greenawalt DM, Dobrin R, Chudin E, Hatoum IJ, Suver C, Beaulaurier J, Zhang B, Castro V, Zhu J, Sieberts SK, et al: A survey of the genetics of stomach, liver, and adipose gene expression from a morbidly obese cohort. Genome Res. 2011, 21 (7): 1008-1016.PubMed CentralView ArticlePubMedGoogle Scholar
Mitchell GF, Verwoert GC, Tarasov KV, Isaacs A, Smith AV, Yasmin , Rietzschel ER, Tanaka T, Liu Y, Parsa A: Common Genetic Variation in the 3-BCL11B Gene Desert Is Associated with Carotid-Femoral Pulse Wave Velocity and Excess Cardiovascular Disease Risk: The AortaGen Consortium. Circ Cardiovasc Genet. 2011, 5 (1):Google Scholar
Bullaughey K, Chavarria CI, Coop G, Gilad Y: Expression quantitative trait loci detected in cell lines are often present in primary tissues. Hum Mol Genet. 2009, 18 (22): 4296-4303.PubMed CentralView ArticlePubMedGoogle Scholar
Hovatta I, Zapala MA, Broide RS, Schadt EE, Libiger O, Schork NJ, Lockhart DJ, Barlow C: DNA variation and brain region-specific expression profiles exhibit different relationships between inbred mouse strains: implications for eQTL mapping studies. Genome Biol. 2007, 8 (2): R25-PubMed CentralView ArticlePubMedGoogle Scholar
Lalonde E, Ha KC, Wang Z, Bemmo A, Kleinman CL, Kwan T, Pastinen T, Majewski J: RNA sequencing reveals the role of splicing polymorphisms in regulating human gene expression. Genome Res. 2011, 21 (4): 545-554.PubMed CentralView ArticlePubMedGoogle Scholar
Smirnov DA, Brady L, Halasa K, Morley M, Solomon S, Cheung VG: Genetic variation in radiation-induced cell death. Genome Res. 2012, 22 (2): 332-339.PubMed CentralView ArticlePubMedGoogle Scholar
Kuchenbauer F, Morin RD, Argiropoulos B, Petriv OI, Griffith M, Heuser M, Yung E, Piper J, Delaney A, Prabhu AL, et al: In-depth characterization of the microRNA transcriptome in a leukemia progression model. Genome Res. 2008, 18 (11): 1787-1797.PubMed CentralView ArticlePubMedGoogle Scholar
Friedlander MR, Chen W, Adamidi C, Maaskola J, Einspanier R, Knespel S, Rajewsky N: Discovering microRNAs from deep sequencing data using miRDeep. Nat Biotechnol. 2008, 26 (4): 407-415.View ArticlePubMedGoogle Scholar
Ambros V, Bartel B, Bartel DP, Burge CB, Carrington JC, Chen X, Dreyfuss G, Eddy SR, Griffiths-Jones S, Marshall M, et al: A uniform system for microRNA annotation. RNA. 2003, 9 (3): 277-279.PubMed CentralView ArticlePubMedGoogle Scholar
Fehniger TA, Wylie T, Germino E, Leong JW, Magrini VJ, Koul S, Keppel CR, Schneider SE, Koboldt DC, Sullivan RP, et al: Next-generation sequencing identifies the natural killer cell microRNA transcriptome. Genome Res. 2010, 20 (11): 1590-1604.PubMed CentralView ArticlePubMedGoogle Scholar
Cummins JM, He YP, Leary RJ, Pagliarini R, Diaz LA, Sjoblom T, Barad O, Bentwich Z, Szafranska AE, Labourier E, et al: The colorectal microRNAome. Proc Natl Acad Sci U S A. 2006, 103 (10): 3687-3692.PubMed CentralView ArticlePubMedGoogle Scholar
Romualdi C, Bortoluzzi S, D'Alessi F, Danieli GA: IDEG6: a web tool for detection of differentially expressed genes in multiple tag sampling experiments. Physiol Genomics. 2003, 12 (2): 159-162.View ArticlePubMedGoogle Scholar
Bail S, Swerdel M, Liu HD, Jiao XF, Goff LA, Hart RP, Kiledjian M: Differential regulation of microRNA stability. RNA. 2010, 16 (5): 1032-1039.PubMed CentralView ArticlePubMedGoogle Scholar
Braconi C, Huang N, Patel T: MicroRNA-dependent regulation of DNA methyltransferase-1 and tumor suppressor gene expression by interleukin-6 in human malignant cholangiocytes. Hepatology. 2010, 51 (3): 881-890.PubMed CentralPubMedGoogle Scholar
Tsai WC, Hsu PW, Lai TC, Chau GY, Lin CW, Chen CM, Lin CD, Liao YL, Wang JL, Chau YP, et al: MicroRNA-122, a tumor suppressor microRNA that regulates intrahepatic metastasis of hepatocellular carcinoma. Hepatology. 2009, 49 (5): 1571-1582.View ArticlePubMedGoogle Scholar
Lewis BP, Burge CB, Bartel DP: Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005, 120 (1): 15-20.View ArticlePubMedGoogle Scholar
Subramanian A, Tamayo P, Mootha V, Mukherjee S, Ebert B, Gillette M, Paulovich A, Pomeroy S, Golub T, Lander E, et al: Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005, 102 (43): 15545-15550.PubMed CentralView ArticlePubMedGoogle Scholar
Ibanez-Ventoso C, Yang M, Guo S, Robins H, Padgett RW, Driscoll M: Modulated microRNA expression during adult lifespan in Caenorhabditis elegans. Aging Cell. 2006, 5 (3): 235-246.View ArticlePubMedGoogle Scholar
Hooten NN, Abdelmohsen K, Gorospe M, Ejiogu N, Zonderman AB, Evans MK: microRNA Expression Patterns Reveal Differential Expression of Target Genes with Age. PLoS One. 2010, 5: 5-Google Scholar
Brosh R, Shalgi R, Liran A, Landan G, Korotayev K, Nguyen GH, Enerly E, Johnsen H, Buganim Y, Solomon H, et al: p53-Repressed miRNAs are involved with E2F in a feed-forward loop promoting proliferation. Mol Syst Biol. 2008, 4: 229-PubMed CentralView ArticlePubMedGoogle Scholar
Dhahbi JM, Atamna H, Boffelli D, Magis W, Spindler SR, Martin DI: Deep sequencing reveals novel microRNAs and regulation of microRNA expression during cell senescence. PLoS One. 2011, 6 (5): e20509-PubMed CentralView ArticlePubMedGoogle Scholar
Bhayani MK, Calin GA, Lai SY: Functional relevance of miRNA sequences in human disease. Mutat Res. 2012, 731 (1–2): 14-19.View ArticlePubMedGoogle Scholar
Kato M, Chen X, Inukai S, Zhao H, Slack F: Age-associated changes in expression of small, noncoding RNAs, including microRNAs, in C. elegans. RNA. 2011, 17 (10): 1804-1820.PubMed CentralView ArticlePubMedGoogle Scholar
Wu S, Huang S, Ding J, Zhao Y, Liang L, Liu T, Zhan R, He X: Multiple microRNAs modulate p21Cip1/Waf1 expression by directly targeting its 3 ' untranslated region. Oncogene. 2010, 29 (15): 2302-2308.View ArticlePubMedGoogle Scholar
Zoncu R, Efeyan A, Sabatini DM: mTOR: from growth signal integration to cancer, diabetes and ageing. Nature Rev Mol Cell Biol. 2011, 12 (1): 21-35.View ArticleGoogle Scholar
Mcgue M, Vaupel JW, Holm N, Harvald B: Longevity Is Moderately Heritable in a Sample of Danish Twins Born 1870–1880. J Gerontol. 1993, 48 (6): B237-B244.View ArticlePubMedGoogle Scholar
Herskind AM, McGue M, Holm NV, Sorensen TIA, Harvald B, Vaupel JW: The heritability of human longevity: A population-based study of 2872 Danish twin pairs born 1870–1900. Human Genet. 1996, 97 (3): 319-323.View ArticleGoogle Scholar
Perls TT, Bubrick E, Wager CG, Vijg J, Kruglyak L: Siblings of centenarians live longer. Lancet. 1998, 351 (9115): 1560-View ArticlePubMedGoogle Scholar
Perls TT, Wilmoth J, Levenson R, Drinkwater M, Cohen M, Bogan H, Joyce E, Brewster S, Kunkel L, Puca A: Life-long sustained mortality advantage of siblings of centenarians. Proc Natl Acad Sci U S A. 2002, 99 (12): 8442-8447.PubMed CentralView ArticlePubMedGoogle Scholar
Martin GM, Bergman A, Barzilai N: Genetic determinants of human health span and life span: progress and new opportunities. PLoS Genet. 2007, 3 (7): e125-PubMed CentralView ArticlePubMedGoogle Scholar
de Lencastre A, Pincus Z, Zhou K, Kato M, Lee SS, Slack FJ: MicroRNAs both promote and antagonize longevity in C. elegans. Curr Biol. 2010, 20 (24): 2159-2168.PubMed CentralView ArticlePubMedGoogle Scholar
Ibanez-Ventoso C, Vora M, Driscoll M: Sequence Relationships among C-elegans, D-melanogaster and Human microRNAs Highlight the Extensive Conservation of microRNAs in Biology. PLoS One. 2008, 3: 7-View ArticleGoogle Scholar
Monks SA, Leonardson A, Zhu H, Cundiff P, Pietrusiak P, Edwards S, Phillips JW, Sachs A, Schadt EE: Genetic inheritance of gene expression in human cell lines. Am J Hum Genet. 2004, 75 (6): 1094-1105.PubMed CentralView ArticlePubMedGoogle Scholar
Smirnov DA, Morley M, Shin E, Spielman RS, Cheung VG: Genetic analysis of radiation-induced changes in human gene expression. Nature. 2009, 459 (7246): 587-591.PubMed CentralView ArticlePubMedGoogle Scholar
Cock PJ, Fields CJ, Goto N, Heuer ML, Rice PM: The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 2009, 38 (6): 1767-1771.PubMed CentralView ArticlePubMedGoogle Scholar
Blankenberg D, Von Kuster G, Coraor N, Ananda G, Lazarus R, Mangan M, Nekrutenko A, Taylor J: Galaxy: a web-based genome analysis tool for experimentalists. Curr Protoc Mol Biol. 2010, Chapter 19: 11-21. Unit 19 10Google Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
