Gene expression profiling of aging reveals activation of a p53-mediated transcriptional program

General

Among the biological classes of transcripts with an age effect, there were several gene families whose mRNA levels were significantly increased or decreased with age [see additional file 1]. There were 8 procollagen genes (11 probe sets) that were altered with aging; 7 of these showed higher levels in the young muscle (col1a1, col1a2, col3a1, col5a1, col6a1, col6a2, col6a3) while only one procollagen gene was increased with aging (col17a1). The keratins, important cytoskeletal structural proteins, were represented by 8 genes (9 probe sets) all of which displayed increased transcript levels in older muscle. Three members (alpha, beta and delta) of the CCAAT/enhancer-binding protein (CEBP) family of transcription factors had increased levels of transcripts in the older mice. The encoded CEBP alpha, beta and delta proteins are important in the regulation of genes involved in immune and inflammatory responses, cell cycle arrest and body weight homeostasis (for a review, see [18]). Because many genes are represented more than once on the MOE430A Affymetrix microarray, several genes had multiple hits on our list of genes with an age effect in skeletal muscle. Among those genes with four or more independent listings were: fibromodulin (6 probe sets, all with reduced expression in the old), transferrin receptor (5 probe sets, all with reduced expression in the old), immunoglobulin kappa chain variable 8 (4 probe sets, all with increased expression in the old), and myristoylated alanine rich protein kinase C substrate (marcks) (4 probe sets, all with reduced expression in the old).

Transcript and Protein Levels of p53-Related Genes Altered with Aging

To confirm the age-related changes in expression observed in p21 and GADD45a using microarrays, we probed young, old and old CR mice using RT-PCR (Figure 1A). We also examined the expression levels of the p53 gene itself, along with the expression levels of another important tumor suppressor, p16. Both p53 and p16 transcripts were considered absent (below levels of background hybridization) in all animals according to the Affymetrix software, but we were able to detect expression using the more sensitive method of RT-PCR. As is shown in Figure 1A, transcript levels for these genes increased significantly (Student's t-test, P < 0.05) with age. Furthermore, old animals undergoing CR demonstrated reduced levels for these genes compared to age-matched mice fed normal diets (only the difference in p53 levels was considered significant at P < 0.05). To examine how the expression of p21, p16, GADD45a and p53 is altered over the course of the adult mouse life span, we used RT-PCR to measure transcript levels for these genes in gastrocnemius muscle at 5 month intervals from 5 to 30 months of age (Figure 1B–E). With the exception of GADD45a (Figure 1D), the increase in expression of these genes does not show a linear pattern with age. p16 and GADD45a expression increases significantly only in the last third of the lifespan, whereas expression of p53 and its transcriptional target p21 occur as early as 10 months of age.

To determine whether the p53-mediated transcriptional program observed in older mouse skeletal muscle translated into higher levels of p53-associated proteins being expressed, we probed young, old and old CR muscle for p53-transcriptional targets p21 and GADD45a. Western immunoblot confirmed a significant increase in p21 protein (4.6 fold, P < 0.005) and GADD45a (1.6 fold, P < 0.05) in older mouse skeletal muscle (Figure 2). Levels for these proteins were reduced in the old CR mice as compared to old mice on the control diet (approximately 57% and 77% of old levels for p21 and GADD45a respectively). An attempt was made with several antibodies to measure p53 protein levels (both phosphorylated and non-phosphorylated isoforms) in GN muscle, but the signal from the Western blots was either absent or too low to accurately quantitate p53 in any age group (data not shown).

Transcript Levels of p21, p16 and GADD45a in Mouse Models of Increased Oxidative Stress

To test whether the changes in expression of the p53-related genes could be attributed to cumulative oxidative damage with age, we compared transcript levels for p21, p16, GADD45a and p53 in C57BL/6 mice with a heterozygous knockout for either glutathione peroxidase 4 (GPX4^+/-) (5 mo old) or MnSOD (SOD2^+/-) (~9 mo old) mutations to age-matched wild-type mice. GPX4^+/-mice have approximately half the amount of mRNA and protein levels for this essential Se-dependent antioxidant protein as compared to wild type mice [19]. Mouse embryonic fibroblasts (MEFs) derived from mice heterozygous for the gpx4 gene show increased evidence of lipid and DNA oxidation, and are sensitive to gamma irradiation [19, 20]. SOD2^+/- mice show a ~50% reduction in MnSOD protein content and MnSOD activity in gatrocnemius tissue [21] and have increased oxidative damage as demonstrated by significantly elevated levels of 8-oxo-2-deoxyguanosine (8oxodG) in both nuclear and mitochondrial DNA in all surveyed tissues compared to wild type mice [22]. RT-PCR analysis of p21, p16, GADD45a and p53, in GPX4^+/- mice showed a significant increase (P < 0.005) in p53 expression compared to the age-matched wild-type mice. We did observe a ~2 fold increase in p16, although it was not statistically significant (P < 0.12). Transcript levels of these genes were not significantly different between SOD2^+/- mice and age-matched controls [see additional file 2]. Thus, constitutive oxidative damage does not always trigger the expression of the genes examined. Possibly, different types of DNA damage in GPX4-/+ and SOD2+/- mice explain the observed differences.

Increased p53 Levels in Older Muscle is Associated with Increased Expression of Genes Involved in p53-Mediated Apoptosis

A critical function of p53 is the induction of apoptosis in cells that have undergone DNA damage or exhibit uncontrolled cellular growth [23]. To investigate whether increased p53 activity with age is associated with transcription profiles indicative of apoptosis in old gastrocnemius muscle, we examined mRNA levels of four genes known to be key transducers of p53-mediated apoptosis: Bcl-2 binding component 3 (bbc3/puma) [24], phorbol-12-myristate-13-acetate-induced protein 1 (pmaip1/noxa) [25], tumor necrosis factor receptor superfamily, member 10b (tnfrsf10b/killer/dr5) [26] and Bcl-2-related ovarian killer protein (bok) [27]. We observed a significant increase in mRNA levels with age for these genes in mouse gastrocnemius muscle and a significant decrease in levels due to caloric restriction in old, age-matched controls (Figure 3A). The expression levels for all of these p53-mediators of apoptosis were also examined in the gastrocnemius muscle over the course of the mouse adult lifespan (Figure 3B–E), and increased levels of transcript were usually observed in the first half of the lifespan. We also compared the levels for these proapoptotic genes in the gastrocnemius tissue of GPX4+/- and the SOD+/- mice to age-matched wild type mice using RT-PCR. Only the levels for tnfrsf10b were considered significantly higher (P < 0.005) in the GPX4+/- mice while expression of these genes was not significantly different in the SOD2+/- mice [see additional file 3]. These findings indicate that p53-dependent mediators of apoptosis are regulated with age and that the mechanism of activation is independent of endogenous levels of oxidative stress.

Increase in mRNA Levels of p53 and p53-Associated Genes in Multiple Tissues of Old C57Bl6/J Mice

Animal treatments

Male C57BL/6NHsd mice were purchased from Harlan Sprague-Dawley at 6–7 weeks of age and individually housed for the remainder of their lifespan. Upon arrival to our facility, mice were fed one of two life-long dietary regimens: normal diet (ND) and caloric restriction (CR). Mice on ND were provided acidic water ad libitum and fed 84 kcal/wk which is ~15% less than the average ad libitum diet. C57BL/6NHsd mice on this dietary regiment typically have an average lifespan of 30 mo. in our colony. Mice on CR were fed a diet that resulted in a 26% difference in caloric intake to the ND animals. The diet fed to CR mice was enriched in protein, vitamins and minerals such that CR and ND mice were fed nearly identical amounts of these components. Detailed information on diet composition and feeding for ND and CR mouse diets have previously been described [49]. Five replicate mice were used for each diet and age group for microarray analysis. Dissected gastrocnemius muscle from the mice were flash-frozen in liquid nitrogen and stored at -80°C until processed for microarray or Western analysis.

Tissue preparation and oligonucleotide array hybridization

RNA was extracted from frozen mouse tissue, transcribed to cDNA which was then used as a template to produce biotin-labeled cRNA as previously described [50]. A total of 10 μg of fragmented cRNA was used for hybridization to each MOE430A array (Affymetrix, Santa Clara, CA). After hybridization, the arrays were washed and stained with signal amplification protocol using antibody. The arrays were scanned twice and the averaged images were used for further analysis.

Data analysis

The complete set of gene expression data has been deposited in the Gene Expression Omnibus database [51] accession # [GSE6323]. Robust Multi-Array Analysis (RMA) was used to pre-process and normalize the raw Affymetrix GeneChip data [52]. RMA is available in Bioconductor [53], an open source and open development software project for the analysis and comprehension of genomic data based primarily on the R programming language [54]. RMA fits a log-scale linear model to derive a single, normalized, summary score of expression for each probe on each array. RMA has been shown to provide greater sensitivity and specificity in detection of differentially expressed genes than the Affymetrix analysis software MAS 5.0 [52, 54]. Because RMA generally provides attenuated estimates of fold-change [55], the values reported in the expression ratios for the array data are likely minimal estimates. MAS5.0 was used to compute absent/present calls. Any gene considered absent across all 15 arrays was removed from the analysis. This left 10,892 probes (out of 22,690).

Statistical analysis

To identify genes differentially expressed over time, we used an empirical Bayes hierarchical modeling approach called EBarrays [56, 57], also implemented in Bioconductor [54]. To briefly illustrate EBarrays, first consider the simple case of two biological conditions. For a given gene g, there are then two underlying levels of expression, one in each condition. These are the quantities of interest that are approximated by the expression measurements. Of most interest is the relationship between these means. In the case of equivalent expression (EE), the two means denoted by μg1 and μg2 are equal; μg1 ≠ μg2 denotes differential expression (DE). We refer to these equality and inequality relationships as expression patterns.

EBarrays assigns a probability to each gene for each pattern. To do this, the probability distribution governing the data must be described. A Gamma-Gamma (GG) model is often used. For that model, we assume the underlying expression levels are approximately distributed as Inverse Gamma; and, for a given gene g in condition k (in our simple example, k = 1,2), measured intensities arise from a Gamma distribution with mean μgk. The Inverse Gamma distribution here describes variation in the true underlying intensity levels of different genes (some genes have low expression, others high) and the Gamma distribution describes within gene variability due to biological and technical error variations. Similarly, a log-Normal Normal (LNN) model can also be used, where the underlying expression levels are assumed Gaussian, and log intensities for gene g in condition k are considered as Normally distributed with mean prescribed by μgk. Since we never know which genes are in which patterns, the full data distribution is given by a mixture: there is some probability p that the data for a given gene arise from pattern DE and a (1-p) chance that the data are EE.

For both cases, model parameters are determined by fitting the model to the full set of array data (details are given in [57]). In this way, the approach utilizes information across genes and arrays to optimize model fit, and is thus more efficient than a number of methods that make gene inferences one gene at a time. Once model parameters are obtained, probabilities for each pattern are calculated for every gene. These probabilities are referred to as posterior probabilities, since they are calculated after observing a set of data. For a given gene and expression pattern, the posterior probability quantifies evidence that the gene is in that pattern. Posterior probability thresholds can be set to control the false discovery rate (FDR) at a desired level [58]. For example, defining DE genes to be those with a posterior probability of DE exceeding 0.95 controls the FDR at 0.05.

This approach, described for two conditions, extends naturally to multiple conditions. Whereas for two biological conditions there are only the EE and DE patterns to consider, for multiple conditions there are many patterns. Given the three conditions here (young, old, and old CR), 5 expression patterns are possible (μ1 = μ2 = μ3; μ1 ≠ μ2 = μ3; μ1 = μ2 ≠ μ3; μ1 = μ3≠μ2; μ1≠μ2≠μ3). Model fit proceeds as in the two-pattern case (the details of the multiple pattern case are also given in [57]). Posterior probabilities for each of the 5 patterns were calculated for 10,892 genes under both GG and LNN assumptions. The appropriateness of model assumptions was checked using diagnostic plots. In most cases, results were similar for both models and only results from the LNN model are reported.

Real-time (RT)-PCR

Total RNA was extracted from frozen tissue using the TRIZOL reagent (Life Technologies, Grand Island, NY). Polyadenylate [poly(A)+] RNA (mRNA) was then purified from the total RNA with oligo-dT-linked oligotex resin (Qiagen, Valencia, CA). Isolated mRNA was reverse transcribed to cDNA and amplified using the TaqMan EZ RT-PCR Kit (Applied Biosystems, Foster City, CA) (5 ng mRNA per 25 μl RT-PCR reaction volume). All primers and fluorescent probes for assayed genes, with the exception of the custom-made p21 and p16, were purchased from Assays-on-Demand Gene Expression probes (Applied Biosystems, Foster City, CA). The gene symbol and Applied Biosystem assay number for all Assays-on-Demand primer/probe sets used for RT-PCR are as follows: gadd45 (Mm00432802_m1), p53 (Mm00441964_g1), bbc3 (Mm00519268_m1), pmaip1 (Mm00451763_m1), tnfrsf10b (Mm00457866_m1), bok (Mm00479422_m1), tbp (Mm00446973_m1), and gapd (Mm99999915_g1). The sequence for the primers and probes for p21 and p16 are as follows: for p21, forward 5'-GGCAGACCAGCCTGACAGAT-3', reverse 5'-TTCAGGGTTTTCTCTTGCAGAAG-3', and probe 5'-TCTATCACTCCAAGGCA-3'; for p16, forward 5'- CCCAACGCCCCGAACT-3', reverse 5'- GCAGAAGAGCTGCTACGTGAA-3', and probe 5'-TTCGGTCGTACCCCGATTCAGGTG-3'. The sequence for the RT-PCR probe and primer for p16 was provided to us by Norman Sharpless's research group. The expression of TATA box binding protein (tbp) and glyceraldehyde-3-phosphate dehydrogenase (gapd), which showed approximately the same expression levels as measured on the microarrays for all animals regardless of age or diet, was used to normalize all RT-PCR results. RT-PCR was performed and quantified using the ABI PRISM 7000 Sequence Detection System (Applied Biosystems, Foster City, CA).

Western immunoblot

Tissues homogenate was prepared from frozen in modified RIPA buffer (Tris-HCL 50 mM pH 7.4, NP-40 1%, Na-deoxycholate 0.25%, NaCl 150 mM EDTA 1 mM) containing protease and phosphatase inhibitors (Sigma). One hundred micrograms of each tissue extract were separated on 12.5% SDS PAGE gel (BioRad) and transferred to PVDF membrane (Pierce). Membranes were blocked in 5% milk in Tris buffered saline 0.1% Tween 20 (TBST) for 1 hour. Membranes were incubated in primary antibody overnight at 4°C with agitation in either 5% milk or 5% BSA in TBST at the indicated dilutions (GADD45A Santa Cruz Biotechnology Inc. 1:500; p21/WAF1 Ab-6 Calbiochem 1:4000; Tubulin Sigma 1:2000). After 3 × 5 min washes, the membranes were incubated in HRP conjugated secondary antibody (Vector labs, 1:10,000 dilution) for 1 hour at room temperature. Membranes were washed as before, incubated in SuperSignal West Pico Chemiluminescent Substrate (Pierce) and exposed to film (Pierce). Band density was calculated using NIH Image J software.