Progressive changes in non-coding RNA profile in leucocytes with age

It has been observed that immune cell deterioration occurs in the elderly, as well as a chronic low-grade inflammation called inflammaging. These cellular changes must be driven by numerous changes in gene expression and in fact, both protein-coding and non-coding RNA expression alterations have been observed in peripheral blood mononuclear cells from elder people. In the present work we have studied the expression of small non-coding RNA (microRNA and small nucleolar RNA -snoRNA-) from healthy individuals from 24 to 79 years old. We have observed that the expression of 69 non-coding RNAs (56 microRNAs and 13 snoRNAs) changes progressively with chronological age. According to our results, the age range from 47 to 54 is critical given that it is the period when the expression trend (increasing or decreasing) of age-related small non-coding RNAs is more pronounced. Furthermore, age-related miRNAs regulate genes that are involved in immune, cell cycle and cancer-related processes, which had already been associated to human aging. Therefore, human aging could be studied as a result of progressive molecular changes, and different age ranges should be analysed to cover the whole aging process.


INTRODUCTION
In the last years, aging has become one of the most studied biological processes, due to its increasing incidence in our population and also to its correlation to other highly prevalent diseases, such as cancer or neurodegenerative diseases [1]. It has been estimated that 30 % of the total population in Europe will be older than 65 years in 2050, becoming a huge challenge from the medical and social point of view.
From a biological viewpoint, aging has been defined as a progressive postmaturational decline in physiolo-gical capacity, accompanied by an increased susceptibility to disease and an increased mortality risk [2].
It is a complex and multifactorial process in which genetic, lifestyle and environmental factors play a role. The genetic configuration of each individual defines how the organism faces our lifestyle and the interaction with the environment in the aging process [1]. Aging could be defined as the phenotypic expression of cellular senescence in our tissues. Cellular senescence is characterized by cell cycle arrest and it is associated to major changes in gene expression patterns. Moreover, AGING aging reaches most tissues and organs, which is measured as accumulation of senescent cells [3].
A high amount of cells of the immune system acquire senescent features, in the process termed immunosenescence [4]. These age-related changes result in increased susceptibility to viral and bacterial infection, poor vaccination response and higher incidence of agerelated diseases, as stated before. Apart from immune cell deterioration, a chronic low-grade inflammation has also been described in the elderly, the so called inflammaging. Inflammaging results from immune dysfunction and accumulation of senescent cells that promote proinflammatory signals, such as elevated secretion of proinflammatory cytokines and activation of NK-κB transcription factor [5,6].
These cellular and serological changes must be driven by numerous changes in gene expression. Previous works have shown that genomic DNA methylation declines with age [7] and also that there is a change in gene expression at transcriptomic level [8] in peripheral blood cells. Additionally, cellular changes could be led by products of the non-coding part of the genome, noncoding RNAs (ncRNAs).
ncRNAs include a wide variety of RNAs that can be classified according to their size in long ncRNAs (>200 bp) or small ncRNAs (<200 bp). Apart from the wellknown ribosomal RNA or transfer RNA, other less characterized ncRNAs are being described to be involved in many cellular processes through the regulation of gene expression. This regulation can be achieved taking part in different processes, such as transcription, mRNA processing or epigenetic control [9].
Among small ncRNA family (sncRNA), microRNAs (miRNA) are the best described class of ncRNAs. They also participate in gene expression regulation, through the binding of target mRNA, triggering either their destabilization or translation inhibition. During the last decade, the implication of miRNAs in virtually all biological processes has been reported, including immune cell development, function [10][11][12] and immunosenescence [13].
Moreover, several works have studied miRNA expression in peripheral blood cells from donors of different ages, including preterm infants, young, adults, octogenarians and centenarians, in order to find healthy aging miRNA expression patterns [14][15][16][17][18]. Three of these reports, include centenarians in their studies, adding a valuable information about human healthy aging [15][16][17]. Conversely, Noren Hooten and colleagues compare old subjects (mean age = 64.5 years) with young subjects (mean age = 30 years) [14] while Lai et al compare preterm infants with adults (mean age= 37 years old) [18]. In all these works, the strategy to find age-related miRNAs has been to compare the expression pattern from the oldest subjects to that from youngest subjects.
Small non-coding RNAs include other subtypes apart from miRNAs and, among them, the small nucleolar RNAs (snoRNAs) have been the most widely studied. They perform sequence-specific 2′ -O-methylation and pseudouridylation of ribosomal RNA which takes place in the nucleolus after forming the small nucleolar ribonucleoprotein (snoRNPs) complex [19]. A role for snoRNAs in cancer has been described and they have also been proposed as good candidates for involvement in other diseases [20]. Recently, it has been reported the presence of snoRNA in human plasma samples and identified a subset of 9 snoRNA showing a decreased expression in individuals being older than 66.5 years [21]. This work supports the idea that several types of non-coding RNA may play a role in human aging.
In the present work, we wanted to focus on how blood miRNA expression changes progressively across the life course, from adulthood to the elderly, apart from the classical analysis of comparing elder subjects' expression profile to that of young subjects'. Using a miRNA microarray, we analyzed the expression of 1769 human miRNA and other sncRNAs and we identified 69 sncRNAs that progressively changed with age. Furthermore, this analysis suggests that 50 years could be a turning point in human aging process.

RESULTS
In the present work we wanted to explore how the global expression of sncRNA changes in human aging. To do that, we analysed the expression of sncRNA by microarrays in leucocytes from healthy donors of different ages, ranging from 24 years to 79 years old ( Figure 2A). We took two approaches to analyse data. On the one hand, we performed a correlation analysis aiming at identifying progressively changing sncRNA and on the other hand, we applied a more classical differential expression analysis. Finally, we carried out validation experiments both in the same discovery cohort and in an independent cohort of 44 donors (age range: 24-87). This study design, which is depicted in Figure 1, allowed us to study the expression of sncRNA in all the stages of the adulthood and part of the elderly.

Correlation and differential expression analysis
The first approach was to assess how the global sncRNA expression changes with chronological age and whether the sncRNA expression was correlated with AGING human aging. The correlation analysis between gene expression and age resulted in 69 sncRNAs (56 miRNA and 13 snoRNA) (Table 1) (Figure 3) that progressively change with age (R>0.75) and have a fold-change greater than 1.5 between the oldest and the youngest sample. Forty-eight sncRNAs show an increasing expression pattern (R>0.75) while 21 have a decreasing expression pattern (R<-0.75) ( Figure 2C). The representation of the mean z-score of expression for the decreasing and increasing sncRNAs ( Figure 2B), shows that at 50 years old the two lines cross, highlighting that this could be a critical age.
In a second approach, we applied a classical differential expression analysis, comparing the expression of sncRNA between elder and young donors. To set up the age at which we should consider the segregation of two groups, we looked at the correlation analysis. According to those results, 50 years seems to be critical, so we set this age as the cut-off between the two groups. Therefore, we had 21 young samples (age ≤ 50) and 17 elder samples (age >50).
Differential expression analysis identified 39 overexpressed and 18 underexpressed sncRNAs in the elder group compared to the young group (57 sncRNAs in total) (limma test, FDR corrected p<0.05)( Table 2). Among those, there are 40 miRNA and 17 snoRNA.
Then, we wanted to compare whether there was any overlap between the sncRNA identified with correlation analysis and differential expression analysis and we found that 25% of the miRNA and 6% of the snoRNA were common between both analyses ( Figure 3).

In silico analysis
After identifying age-related miRNAs, we wanted to see whether those miRNAs had already been associated to aging. Among the 56 miRNAs that progressively change with age (correlation analysis), 18 had already been associated with aging process and interestingly, 11 had the same expression pattern as in the present study (increased or decreased with age). When we move to the 40 differentially expressed miRNAs between extremes of the cohort (older vs young, differential expression analysis) and compare them with the database, we also found 18 overlapping miRNAs, 8 being in the same direction of deregulation ( Figure 4 and Table 3).

AGING
As stated before, the main function of miRNAs is to regulate gene expression post-transcriptionally and they act binding to their target mRNAs and causing mRNA degradation or translation inhibition. Therefore, to know which functions were affected by the 56 miRNAs that progressively increase or decrease with age, we searched them in miRTarBase database and found 9370 interactions for 54 miRNA with 4082 genes. Among all those interactions, 458 were supported by strong experimental evidences and included 222 unique genes regulated by miRNAs that decrease with age and 26 unique genes regulated by miRNAs that increase with age ( Figure 4). These target genes are involved in immune, cell-cycle and cancer related processes, UV-d image response and other biological processes ( Table 4).
As for age-related snoRNA, even if they are not included in Digital Aging Database, there are three (U43, U52 and U91) which had already been related to aging process by other authors ( Table 3).

Validation of candidate miRNA
Once we identified a subset of miRNA and snoRNA related to human aging in leucocytes, we validated some candidates in an independent cohort. The selection of the miRNA for validation was made on the basis of fold-change and previous association with aging. Therefore, taking the list of 69 sncRNA that showed to be correlated with aging, we selected those having a fold change higher than 2 (up or downregulated) between the last and the first sliding window, which yielded a list of 15 sncRNA. Finally, we applied an additional filter, selecting those sncRNA that had already been associated to aging according to Digital Aging Atlas. This selection resulted in a list of 5 candidate miRNA: miR-15a, miR-30b, let-7i, let-7g and miR-1281 (Table 5).  AGING First, we assessed by RT-qPCR the expression of candidate miRNA in the same discovery cohort of microarray analysis as a technical validation, which showed that the expression of let-7i and let-7g was significantly lower in the group which was over 50 years (RQ=0.566; p-value= 0.0286 and RQ= 0.578; p-value=0.0281, respectively) ( Table 5). miR-15a and miR-30b also had a slightly lower expression in this group, but the difference was not statistically significant (RQ=0.632; p-value=0.121 and RQ=0.777; p-value=0.0591, respectively) ( Table 5). On the other hand, miR-1281 did not show any difference in its expression between the two groups (RQ=1.127; p-value=0.672) ( Table 5). This RT-qPCR analysis confirmed the results observed in microarray experiment, except for miR-1281 ( Figure 5). Even if miR-15a and miR-30b did not show a statistically significant change, we decided to include them in further validation experiments. Therefore, we studied the expression of the four downregulated miRNA in an independent cohort and found a statistically significant underexpression in subjects being older than 50 years ( Figure 5A and Table 5). AGING Furthermore, we segregated the cohort in five age ranges to see whether the expression decrease was a gradual process with age, as identified by microarray experiment. We created 5 groups as follows: 24-40 (n=5), 41-60 (n=10), 61-70 (n=7), 71-80 (n=15) and 81-90 (n=7). Then, we represented the expression values of each miRNA (in DCq) for each group and it can be observed that their expression decreases gradually in this ranges, specially in miR-15a, let-7i and let-7g. However, it is true that the group composed of individuals between 80 and 90 years old, have a similar expression to that of 71-80 ( Figure 5B).

DISCUSSION
In the present work we have identified a sncRNA signature that shows a progressive expression change with age in peripheral blood leucocytes. This signature involves 69 sncRNAs (56 miRNAs and 13 snoRNAs), 48 of them showing an increasing expression while 21 exhibit a decreasing expression pattern. Our data present that the expression of this sncRNAs subset starts to change around 40 years old and at 47 years, this alteration suffers and acceleration, the main change occurring from 47 to 54 years old. These results highlight that the expression change occurring with age is a progressive process rather than being a sudden event happening in the elderly. Yet, it is true that this progressive change does not follow a uniform pattern, given that we have observed an acceleration between 47 to 54 years.
This expression pattern described for sncRNAs, is similar to that found with mRNA expression analysis previously by our group [8]. Both studies have found a subset of RNAs whose expression changes progressively with chronological age. Yet, the most prominent change described for mRNA starts at a later age than that seen in small RNAs (49 years old in transcriptome and 47 in miRNome) and it is slightly more pronounced. Furthermore, both works determine the critical age in the sixth decade of life, based on the age at which increasing and decreasing mRNA and sncRNA cross ( Figure 2B).
The correlation analysis carried out in the present work has been able to identify a subset of 69 probes, the expression of which gradually changes with chronological age. Furthermore, a more classical differential expression analysis has stood out an overlap of 35.2% and 25% of miRNA and snoRNA, respectively, with the correlation analysis, highlighting the validity of the correlation analysis to identify agerelated miRNA. It is true, however, that there is an important fraction of sncRNA that have been identified only with one of the analyses. This observation suggests that each of the methods is able to address different questions: the correlation analysis would identify genes showing a gradual expression change correlated with chronological age, and would fail to detect those that suffer a big sudden change in the elderly, whereas differential expression analysis would identify genes presenting the opposite situation.
Among the 71 miRNAs and 24 snoRNAs identified in the present study (with both types of analysis), 46 have previously been reported to change with age by different authors and in different tissues or cell types. Interestingly, 16 sncRNA identified in the present study are in the same direction of change as reported previously (Table 3). Surprisingly, some of those miRNAs had been identified in an in vitro model of foreskin and our research highlights that those miRNAs may have an impact also in the aging of immune system cells. The rest of the miRNAs had also been described in blood cells, reinforcing the role of those miRNAs in identified by correlation and differential expression analysis and the overlapping probes between two analyses.
AGING senescence and reflecting the validity of our approach. Of note, let-7a, let-7g, let-7i, miR-126, miR-130a and miR-19b have been reported to be decreased and miR-423-3p increased in the elderly by more than one author apart from the present study, supporting their role in healthy aging (Table 3). This observation, points at those miRNAs as candidate molecules to be studied in the future in order to stablish their biological function and understand their implication in human healthy aging.   AGING Moreover, we have confirmed the decreasing expression of miR-15a, miR-30b, let-7i and let-7g with age in and independent cohort by RT-qPCR. As stated before, those miRNA had already been associated to human aging by several authors. Hackl and colleagues, for instance, had also found the downregulation of these miRNA in several cell types and tissues including T cells, which is in accordance with our results [22]. Serna et al, had also related those miRNA to human aging but, they found and overexpression of miR-30b, miR-15a and let-7i in centenarians compared to octogenarians [17]. These data show the opposite trend than what we have observed in our cohort, where those miRNA appear to decrease with age. Nonetheless, they have studied centenarians, octogenarians and young people and they conclude that as for miRNA expression is concerned, centenarians are more similar to young people than to octogenarians. Therefore, the different trend they observed in miR-30b, miR-15a and let-7i expression comparing to our results could be conferred to the fact that the expression of those miRNA in centenarians is very similar to young people.
In order to assess the possible implication of age-related miRNA in human aging, miRNA-target interaction analysis was carried out, which revealed that age associated miRNAs regulate genes involved in immune, cell-cycle and cancer related processes, UV-induced damage response and others. Interestingly, many of the enriched pathways in age-related decreasing miRNA target genes are the same as those found among increasing miRNA targets. Therefore, those biological processes might be regulated either by repressing the expression (through increasing miRNA) of some genes or reducing miRNA-exerted repression (through decreasing miRNA) of some other genes in the same pathways.
Inflammatory pathways had already been associated to aging process [1,23] and our results indicate that some of those pathways are regulated by miRNAs. For instance, NFKB1 appears to be a target of miR-146a-5p, which decreases with chronological age, being a candidate that would explain the increase in this transcription factor activity that have already been described in immunosenescence [6]. The increase of IL-6 and TNFα plasma levels have been described in elder people as well [24][25][26][27][28], and our results point out that IL-6 and TNFα induced signalling pathways via NFkB might be regulated by a subset of miRNAs.
Moreover, the target genes of miRNAs that exhibit a gradual change with chronological age (found in miRTarBase and before filtering by strong evidence of interaction) include several genes whose expression also show a progressive pattern with age, as found in our previous work [8]. Namely, 40 genes (21 downregulated and 19 upregulated) appear to be agerelated miRNA targets and show a concordant progressive expression change (increasing miRNA and decreasing target gene, and vice versa) (Figure 4). This observation points out that age-related miRNAs identified in the present work might be responsible for the expression pattern described for a subset of proteincoding genes.
In conclusion, our study identifies a subset of 69 small non-coding RNAs (54 miRNAs and 15 snoRNAs) of which expression in peripheral blood leucocytes changes progressively with chronological age. Furthermore, our data suggest that the age range from 47 to 54 is critical from a molecular point of view, given that sncRNA expression changes start to accelerate. Moreover, these miRNAs regulate several genes that are involved mainly in cancer and immune system-related pathways, which had already been associated to human aging. Therefore aging could be studied as a result of progressive molecular changes, and different age ranges should be analysed to cover the whole aging process.

Sample donors and sample collection
In order to perform the study, healthy donors were recruited in the Department of Neurology at Donostia University Hospital and the Basque Biobank. The discovery cohort included 38 donors with an age range from 24 to 79 and the validation cohort included 44 donors ranging from 24 to 87 years old. Donors who were over 60 years were examined by a neurologist to ensure their healthy status. Samples from all donors were collected after receiving written informed consent. The study was approved by the hospital's ethics committee.
Peripheral blood (10 ml) was collected by venipuncture after receiving informed consent. Blood samples were collected in tubes containing EDTA as anticoagulant (Vacutainer, Becton Dickinson) and samples were processed and stored at the Basque Biobank.

RNA isolation
Total RNA was isolated from peripheral blood leucocytes with the LeukoLOCK kit (Ambion, ThermoFisher Scientific, 168 Third Avenue, Waltham, MA 02451, USA) using the alternative protocol to capture small RNAs. RNA concentration was measured using a NanoDrop ND-1000 spectrophotometer and RNA integrity was assessed using Agilent 2100 Bioanalyzer with the RNA 6000 Nano Assay Protocol (Agilent Technologies, 5301 Stevens Creek Blvd. Santa Clara, CA 95051, USA). We only included samples with a RNA integrity number higher than 6.

Data normalization and filtering
We first performed raw data analysis, including a detection step (a probe set is detected above background with an associated p-value) resulting in a true/false call and a quantile normalization step using the miRNA QC Tool software v1.0.33.0 (Affymetrix, 3420 Central Expressway, Santa Clara, CA 95051, USA). In a subsequent filtering step, all non-human probesets were removed, resulting in a final list of 1769 probesets for downstream analysis. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [29] and are accessible through GEO Series accession number GSE89042 (www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89 042).

Correlation analysis
The workflow followed in this study for data analysis is shown in figure 4. In order to detect the genes whose expression progressively increase or decrease with age, the Pearson's correlation was calculated between the age and each of the sncRNA expression. To do that, we first performed a sliding window analysis. Briefly, we ordered the samples by age from the youngest to the oldest, made a group of first 10 samples (first window) and calculated the average expression for each probeset in this window. Then, we slid the window one sample to the right and calculated the average expression for each probeset in this second window. We carried out this computation to a total of 28 windows, obtaining at the end a matrix of 28 values for each probeset, being each value the average expression of 10 samples. Z-scores were computed for this matrix and these values were used for calculating the correlation with age. Finally, to identify the most significant sncRNAs, we selected those having a R > 0.75 and with a > 50% expression difference between the first and the last sliding window.

Differential expression analysis
With the aim of identifying any sncRNA that was differentially expressed between elder subjects and young subjects, these two groups were compared using the "limma" package from Bioconductor [30]. According to the correlation analysis, the cut-off to establish the two groups (young and elder) was set at age 50. Differentially expressed sncRNAs were defined as those passing the false discovery rate (FDR) corrected p-value, set at 0.05 (Figure 4). Venn diagrams were created using Venny 2.1 (http://bioinfogp.cnb.csic.es/tools/venny/index.html).

Search for microRNA targets and biological function
In order to get the targets of miRNAs, miRTarBase database was used, which gathers experimentally validated miRNA-target interactions (MTI) [31]. First, using a tool available at miRandola web (http://atlas.dmi.unict.it/mirandola/index.php) [32] we converted the miRNA names from miRBase v11 (in which the microarray used is based) to miRbase v20 to match the nomenclature of miRTarBase. Afterwards, we found the miRNA-target interactions registered in miRTarBase for the miRNAs of interest and removed the interactions that have weak evidence according to the database. Afterwards, we computed overlaps of miRNA target genes with Hallmark gene sets from Molecular Signature Database v5.0 on Gene Set Enrichment analysis website [33], considering only the resulting top twenty gene sets from significant list (FDR < 0.05).
Finally, we wanted to evaluate whether the miRNAs identified in the present study had already been related to aging process in other works. To do that, we compared our age-related miRNAs with those compiled in Digital Aging Atlas [34] and the article by Serna et al [17], which is not included in this database. The snoRNA are not included in this database either.

Validation of candidate miRNA
Validation of candidate miRNA were performed by reverse transcription quantitative PCR (RT-qPCR) using miScript PCR system (Qiagen, QIAGEN Strasse 1 40724 Hilden, Germany). Total RNA (200 ng) was reverse transcribed using miScript II Reverse Transcription kit using the HiSpec buffer and following the manufacturer's protocol in a Veriti thermal cycler (Applied Biosystems, ThermoFisher Scientific, 168 Third Avenue, Waltham, MA 02451, USA). Afterwards, qPCR was carried out with miScript SYBR Green PCR kit, adding 3 ng of cDNA as input, as specified in the protocol. We used the following miScript primer Assays to amplify miRNAs of interest: Hs_let-7g_2, Hs_let-7i_1, Hs_1281_1, Hs_miR-15a_1, Hs_miR-30b_1 and Hs_miR-191_1 as the reference gene. The PCR was carried out in a CFX384 thermal cycler (BioRad, 1000 Alfred Nobel Drive, Hercules, CA 94547, USA).
Raw data were processed in Bio-Rad CFX Manager software and the subsequent analysis to calculate relative expression was carried out in Excel software using 2 -DDCT method [35]. Statistical analysis and plots were performed in RStudio software version 1.0.44 running R version 3.3.2. To compare the expression of candidate miRNA between donors who were over 50 years and those who were under 50 years old we applied a Wilcoxon runk sum test.