MicroRNA, mRNA, and protein expression link development and aging in human and macaque brain

Abstract

Changes in gene expression levels determine differentiation of tissues involved in development and are associated with functional decline in aging. Although development is tightly regulated, the transition between development and aging, as well as regulation of post-developmental changes, are not well understood. Here, we measured messenger RNA (mRNA), microRNA (miRNA), and protein expression in the prefrontal cortex of humans and rhesus macaques over the species' life spans. We find that few gene expression changes are unique to aging. Instead, the vast majority of miRNA and gene expression changes that occur in aging represent reversals or extensions of developmental patterns. Surprisingly, many gene expression changes previously attributed to aging, such as down-regulation of neural genes, initiate in early childhood. Our results indicate that miRNA and transcription factors regulate not only developmental but also post-developmental expression changes, with a number of regulatory processes continuing throughout the entire life span. Differential evolutionary conservation of the corresponding genomic regions implies that these regulatory processes, although beneficial in development, might be detrimental in aging. These results suggest a direct link between developmental regulation and expression changes taking place in aging.

Figure 1.

mRNA, protein, and miRNA expression changes during life span. (A–E) The first two principal components of mRNA, miRNA, and protein expression in human and rhesus macaque brains. The analysis was performed by singular value decomposition, using the “prcomp” function in the R “stats” package, with each gene scaled to unit variance before analysis. The numbers represent each individual's age in years. The proportion of variance explained by each principal component is shown in parentheses. For mRNA, arrows indicate pairs of technical replicates, and shades of blue represent two experimental batches. Individuals group according to their age, indicating a substantial influence of age on total expression variation. (F) Distribution of Pearson correlation coefficients between human and macaque expression time series, calculated for 3233 orthologous mRNAs (blue) or 98 orthologous miRNAs (red), showing significant expression change with age in human (Supplemental Text S1). The y-axis shows the relative frequency (the Gaussian kernel density estimate, calculated with the R “density” function) of genes showing a certain Pearson correlation coefficient.

Figure 2.

Major patterns of mRNA and miRNA changes with age. (A) Shows the average expression levels in eight coexpressed gene groups (see Methods). Dark blue points represent the mean expression level of all genes in a group per individual. The y-axis shows standardized expression levels, where each unit indicates one standard deviation difference from the mean. The x-axis shows age of individuals on (age)^¼ scale (i.e., the fourth root scale, which provided optimal resolution of expression changes during both developmental and adulthood periods). Bold vertical bars indicate the 25%–75% quantile range, and thin bars indicate the 2.5%–97.5% quantile range. Mean expression change with age within each group is summarized by spline curves. Green circles, green triangles, and purple crosses show the mean expression levels for the same genes among individuals from two published mRNA data sets (Lu et al. 2004; Somel et al. 2009) (GEO data set accession nos. GSE1572 and GSE11512), and the protein data set from the present study, respectively. (B) Major patterns of miRNA changes with age. Labels are as in A.

Figure 3.

Transition points of expression change with age for mRNA, miRNA, and proteins. (A) The age distribution of expression transition points determined on gene-by-gene basis. The y-axis shows the relative frequency (the Gaussian kernel density estimate, calculated with the R “density” function) of genes showing a certain transition point. The x-axis shows transition ages on the (age)^¼ scale. (Blue) Human mRNA; (green) macaque mRNA; (purple) a published human mRNA data set (Somel et al. 2009); (gray) human protein; (red) human miRNA; (orange) macaque miRNA. Only age-related genes following nonlinear trajectories and showing significant transition points are represented (Supplemental Table S3). (B) The transition point identification procedure illustrated using genes in groups 4 and 6 (as shown in Fig. 2A). The y-axes indicate mean normalized expression levels of genes in the group. The x-axes show individuals' ages in log₂ scale, allowing improved resolution of developmental changes (Methods; Supplemental Fig. S14). Blue points represent expression levels from the human mRNA data set; blue solid lines show spline curves fit to these data. Blue vertical lines show the transition points. Dotted blue lines show linear regression of expression on age before and after the transition point. Purple and brown points/lines represent mean expression levels/linear regression lines from two published data sets (Lu et al. 2004; Somel et al. 2009), respectively. Note that the results shown in A are calculated per gene, and in B using the means of gene groups.

Figure 4.

miRNA and TF regulation of expression changes with age. (A) Excess of negative correlations between miRNA/TF–target pairs. (Left panels) The y-axes show the relative frequency (the Gaussian kernel density estimate, calculated with the R “density” function) of Pearson correlation coefficients. Shown are correlations between expression profiles of age-related regulators (miRNA or TF) and their age-related target genes (Methods). Colored curves represent the distribution of regulator–target correlations for miRNA–mRNA (red), miRNA–protein (orange), and TF–mRNA (blue). Gray curves show the background distribution: correlations between regulators and non-targets (genes with no evidence of being targeted by the respective regulators). (Right panels) The difference between the kernel density distributions of regulator–target correlations and the background. The gray lines represent 100 simulation results, generated by randomly selecting the same number of background pairs, as regulator–target pairs. The bimodality of the correlation coefficient distributions is because we calculate correlations between age-related regulators and targets only; so, each pair shows some degree of correlation, positive or negative. We therefore test the excess of negative correlations for predicted miRNA–target pairs, relative to randomly paired age-related miRNA and mRNA. Similarly, we test the excess of strong positive and negative correlations for the predicted TF–target pairs (given the dual role of TFs as activator and/or repressor of transcription), relative to randomly paired age-related TFs and mRNA. (B) The proportion of expressed miRNA or TFs showing target enrichment among eight coexpressed gene groups (at HT, P < 0.05).

Figure 5.

miRNA and TF regulation in development and aging. (A) Excess of negative correlations among selected miRNA–target pairs in three coexpressed gene groups in human or macaque cortex. The colored bars show the proportions of negative correlations among miRNA with significant target enrichment within a gene group (at HT, P < 0.05) and their targets in that group at different correlation cutoffs. Hatched bars indicate the proportions of negative correlations among miRNA without target enrichment in a gene group and their targets in that group (the background). The asterisks indicate support for observed–background difference, calculated by bootstrapping the background set 1000 times; ***P < 0.001; **P < 0.01; *P < 0.05; ^o P < 0.10. Both observed and background correlations are calculated separately for developmental and aging periods. The names of identified putative regulatory miRNA are shown above each gene group. Genes in group 1 show limited expression change during aging; therefore, we do not estimate regulators for this group at this period. For macaque, regulators shown in the figure are predicted based on the macaque data and independently of human analysis results, using the same significance levels. Additionally, ∼80% of regulators predicted in humans show a tendency for negative correlation with their targets in macaques (Supplemental Table S5). (B) Excess of negative and positive correlations among TF–target pairs in three gene groups. The colored and hatched bars represent proportions of observed and background TF–target gene correlation pairs, as in B, but the x-axis shows the absolute Pearson correlation cutoff. The names of identified putative regulatory TFs are shown above each gene group. (C) A network of regulatory interactions identified in groups 4 and 6. Only part of the full network, listed in Supplemental Table S5, is shown. The represented genes are those containing the specific miRNA binding site and that show significant negative correlation with that miRNA's expression profile, either in development (green edges) or aging (blue edges). The figure was drawn using Cytoscape software (v 2.6.3).

Figure 6.

Functions, regulation, and specificity of coexpressed gene groups. Shown are mean expression profiles of selected genes within coexpressed gene groups, and their putative miRNA regulators. The empty triangles show mean standardized (z-transformed) human mRNA (blue) and miRNA (red) expression levels, while empty circles show mean standardized macaque mRNA (green) and miRNA (orange) expression levels (note the differences in timing of expression changes between human and macaque, for both mRNA and miRNA expression). The x-axis shows age of individuals on the (age)^¼ scale. The lines correspond to cubic spline curves. The depicted genes are associated with specific Gene Ontology functional terms significantly enriched within the given coexpressed group. For A, C, and E, the genes are further targeted by specific miRNAs. (A) miR-29a and its four cancer-related targets in group 1 (MMP2, TRAF4, COL4A2, COL4A1). (B) Seventeen genes involved in electron transport in group 5. (C) miR-222 and its target in group 6, REV1, involved in DNA damage repair. (D) Fifty-seven neuronal genes in group 8. (E) miR-34a and its seven target neuronal genes in group 4 (GREM2, CAMSAP1, TANC2, CALN1, RGMB, FKBP1B, RTN4RL1). (F) Cell-type specificity of gene groups. The y-axis shows the percentage of cell-type-specific genes among the eight coexpressed age-related gene groups (based on Cahoy et al. 2008; see Methods).

Figure 7.

Diminishing stabilizing selection pressure with age. (A) The symbols and fitted spline curves show the stabilizing selection scores (SSS) calculated for protein coding (blue and green), promoter (orange), and 3′-untranslated (UTR) (red) regions (Methods). The SSS indicate correlation between conservation values and standardized expression levels per individual, across 4084 age-related genes. Conservation scores are corrected for variation in mutation rates. The x-axis shows age of individuals on the (age)^¼ scale. Positive SSS indicate above-average correlation between expression levels and sequence conservation among genes, at a certain age. The dashed vertical line indicates 20 yr of age, when brain maturation is largely complete (de Graaf-Peters and Hadders-Algra 2006). (B) Same as A, but excluding genes possibly under positive selection (Methods). (C) Same as A, but only using genes with enriched expression levels in neurons. (D) Correlation between standardized expression levels and potential confounding factors across age-related genes: number of protein–protein interaction partners (blue), number of tissues (gray) or cell types (brown) a gene is expressed in (i.e., expression breath).