Affymetrix Mouse Genome 430 2.0 oligonucleotide arrays, which include probe sets for more than 39,000 transcripts selected from GenBank, dbEST and RefSeq. Mouse Genome 430 2.0 array annotation was obtained from release 31 of the NetAffx database [67]. Raw data have been deposited in the NCBI Gene Expression Omnibus database and are accessible through accession number GSE35322. Statistical analysis of aging effects in human and mouse skin Quality of array hybridizations was evaluated according to standard metrics, including degradation scores based upon 39 to 59 trends in expression measurements, average background and the percentage of transcripts with detectable expression levels [70]. Expression scores for all arrays were calculated using robust multichip average (RMA) [70], and both human and mouse skin datasets were adjusted for batch effects using Empirical Bayes methods [71]. In both mice and humans, functions from the Limma linear modeling package were used to estimate model coefficients associated with aging effects [72], and FDR-adjusted p-values were calculated using the Benjamini-Hochberg method [73]. For mouse tail skin data, aging effects were evaluated by a two-sample treatment comparison between young and old mice [72]. With respect to human sun-protected skin, expression values were available across a continuous span of ages. The effects of aging on gene expression were thus evaluated using iteratively reweighted least squares robust regression [74]. This approach estimates model coefficients by iteratively minimizing a loss function that assigns reduced weight to larger residuals, thereby attenuating the influence of outlying observations [74]. To implement this method, we used the ‘‘robust’’ option in the ‘‘lmFit’’ function from the Limma linear modeling package (available from Bioconductor) [75], which provides an interface to the MASS package ‘‘rlm’’ function with default settings [76]. Default settings for the ‘‘rlm’’ function initialize the procedure by obtaining ordinary least squares coefficient estimates, with subsequent refitting and reweighting of residuals using the ‘‘Huber function’’ with coefficient value of 1.345 [74]. Weights for individual observations were not provided during the initiation step, such that all observations were weighted equally when obtaining the initial least square coefficient estimates [76]. Both human and mouse analyses were performed using RMA- derived scores in which gene expression is quantified on a log2 rather than absolute scale. In human regression analyses, we expected that expression estimates on the log2 scale would have a linear association with individual age. For each of the 54,675 probe sets represented on the Human Genome U133 Plus 2.0 array, we evaluated whether including a quadratic term in our models led to a significant decrease in model deviance (i.e., F-test for lack of fit) [74], but did not obtain significant evidence in favor of this type of non-linear association between log2-based gene expression and individual age (FDR-adjusted P$0.99 for both sexes). This supported adequacy of the linear regression approach in this context, and we expected that any opposing trends for some genes would be mitigated by our use of robust regression analysis [74,76]. We note, however, that we have estimated 40-year fold- change in humans and 2-year fold-change in mice on an absolute scale (e.g., see Figures S1 and S10). For these fold-change estimates, we have assumed that gene expression changes with age are linear on an absolute scale. Since this assumption may not be valid for all genes, our fold-change estimates should be interpreted as a summary measure, which is accurate in cases for which gene expression on an absolute scale changes in linear fashion with age, but represents only the average rate of change across the lifespan Mechanisms of Transcriptome Aging in cases for which gene expression on an absolute scale changes in non-linear fashion with age. Overrepresentation of gene ontology (GO) terms within specific sets of genes was evaluated using a conditional version of Fisher’s exact test [70,77]. With this approach, significance of a GO term is evaluated based upon a filtered set of associated genes, which excludes any gene associated with both that term and a more specific child term that is already significantly overrepresented. The main advantage of this approach is that results are decorrelated with less redundancy within a set of GO terms identified as significantly overrepresented with respect to genes increased or decreased by aging. We note, however, that significant GO terms are identified based upon p-values haven’t been corrected for multiple hypothesis testing, since an appropri- ate p-value correction has not yet been established in this context [70]. In all GO analyses (Figures S2 and S11), overrepresentation of terms among age-regulated genes was evaluated with respect to a background of genes with detectable expression in the majority of human or mouse skin samples (12,500 human genes; 14,681 mouse genes). Presence or absence of transcripts in human or mouse skin samples was evaluated based upon p-values generated from the Wilcoxon signed rank test as implemented in the MAS 5.0 algorithm [70,78]. Meta-analysis of age-associated gene expression patterns Meta-analysis of age-associated gene expression patterns was carried out using data available from either the Gene Expression Omnibus or ArrayExpress databases (accession numbers for human data: GSE876, GSE5086, GSE9419, GSE9103, GSE11906, GSE11375, GSE11882, GSE16028, GSE19743; accession numbers for mouse data: GSE3150, GSE3253, GSE4332, GSE4786, GSE6591, GSE8150, GSE10000, GSE10965, GSE11291, GSE11667, GSE13799, GSE21716, GSE22317, E-MEXP-839, E- MEXP-1504; Tables S1 and S2). Additionally, we have utilized a complete dataset provided as supplemental material from a previous study [45]. All datasets used the same Affymetrix oligonucleotide array platforms that we have utilized to generate expression data in human and mouse skin (i.e., the Human Genome U133 Plus 2.0 and Mouse Genome 430 2.0 array platforms). Statistical analysis of aging effects was carried out using two-sample comparisons (young versus old) or robust linear regression with age as a single predictor variable (see above). We have not incorporated cDNA array data from the AGEMAP project [24] into our analysis because the AGEMAP cDNA platform provides only partial coverage of the mouse genome (i.e., 8932 mouse genes), because AGEMAP cDNA array data have already been meta-analyzed in a previous study [11], and because differences between Affymetrix oligonucleotide and cDNA array platforms render direct comparisons difficult to interpret [79–81]. Calculation of human and mouse inflammation profiles Procedures followed for calculation of inflammation profiles have been described in a previous publication [31]. To generate human profiles (Figure 1), a database of 1382 Affymetrix Human Genome U133 Plus 2.0 array hybridizations was manually assembled using annotated CEL files downloaded from Gene Expression Omnibus. In total, the 1382 CEL files belonged to 309 cell population treatments, where each treatment contained two or more replicate CEL files associated with a given cell population (average of 1382/309 = 4.47 replicates per treatment). Signature transcripts associated with each treatment were identified by comparison to an independent reference set of 61 samples (representing treatment A in Figure 2 from [31]), consisting of array hybridizations involving a heterogeneous collection of PLoS ONE | www.plosone.org 14 March 2012 | Volume 7 | Issue 3 | e33204