Proteome-wide analysis reveals an age-associated cellular phenotype of in situ aged human fibroblasts

Proteome-wide characterization of in situ aged and ex vivo cultured dermal fibroblasts

First, we analyzed proteins expressed in ex vivo cultured human dermal fibroblasts by a proteomic approach using quantitative mass spectrometry (MS). Using a label-free MS approach, we identified a total of 2409 proteins analyzing 15 primary fibroblast cultures obtained from each five young (20-30), middle aged (40-50) and old (60-70) female donors. Due to high reproducibility of our label-free MS analysis we were able to consider 1607 proteins for quantitative analysis. Statistical analysis revealed 43 proteins with a significant age-associated expression (p ≤ 0.05), whereof 20 and 23 proteins exhibit positive and negative correlation with age, respectively (Tab. 1 and Tab. S1). However, principal component analysis revealed an almost homogeneous protein expression pattern of the three age groups with a low separation within the aged samples (Fig. S1). By detailed analysis of the proteome data we revealed that over 47% (755 proteins) of quantified proteins exhibited a constant abundance across all age groups and thus likely account for this homogeneous protein expression pattern (Fig. S2). By means of enrichment analysis we found that fibroblasts maintain protein abundance for biological processes such as ‘translation’, ‘metabolism of proteins’ and ‘metabolism of RNA’ during aging (Tab. S2). To reveal proteomic changes beyond the housekeeping processes, we analyzed the proteome data considering proteins’ abundances. Therefore, we ranked the quantified proteins into four abundance classes over a dynamic abundance range of approximately five orders of magnitude (Fig. 1). This approach revealed that groups of proteins linked to specific biological processes demonstrated coordinated abundance changes during aging. For instance, a group of five proteins (MTCO2, NDUFA5, NDUFA9, NDUFA10, NDUFS6) linked to the biological process ‘mitochondrial ATP synthesis coupled electron transport’ appeared among the semi-abundant proteins (class III) for the young and middle age group, but among the low-abundant proteins in the old aged donors (class IV; Tab. S3-5). For the other abundance classes (class I-II) we found no age-associated changes. Due to observed age-associated changs of mitochondrial proteins we performed a detailed analysis of proteins related to mitochondria. Cluster analysis revealed 23 age-associated mitochondrial proteins (p ≤ 0.1; 14 negative and 9 positive). Under the proteins with negative age-correlation the biological processes ‘ATP synthesis’ and ‘cellular respiration’ were significantly enriched. The biological processes ‘response to stress’ and ‘catabolism’ were enriched for proteins with a positive age-correlation (Fig. S3; Tab. S6). Proteome-wide characterization revealed that in situ aged fibroblasts exhibit a moderate age-associated cellular phenotype with a large number of proteins exhibiting a constant abundance across the three age groups.

Transcriptome-wide characterization of in situ aged and ex vivo cultured dermal fibroblasts

Next, we considered gene expression data from the same in situ aged fibroblasts already published by Kalfalah and colleagues [10]. It was reported that about 17,000 transcripts were successfully quantified, whereof a total of 137 genes exhibited a significant age-associated differential expression (Tab. S7). Comparison between the proteome and transcriptome data revealed that for almost all proteins (98%) gene expression data were available, but on the contrary the coverage for the 137 age-associated gene by proteome data was quite low (Fig. 2A). Only eight of the 137 differentially changed genes were successfully quantified at the proteomic level and exhibited no agreement with the respective mRNA data concerning age-associated differential expression. Due to intraspecies comparison with an 87% and 98% correlation within the proteome and transcriptome data, respectively, we can exclude technical problems for the minor overlap of the expression data (Fig. S3). The interspecies comparison with about 30% correlation between proteins and mRNAs data was in the expected range (Fig. 2B; Fig. S4). Further detailed analysis revealed that over 63% (10665 quantified mRNAs) of the analyzed genes exhibited constant expression across the different age groups (Fig. S5). In contrast, for all of the 43 age-associated proteins we obtained gene expression data and found out that 77% (33 proteins) of the corresponding transcripts were not changed during aging. Comparison of these age-associated proteins with transcriptome data revealed that age-associated proteins were not differentially regulated at the transcript level.

MicroRNA profiling of in situ aged and ex vivo cultured dermal fibroblasts

Then, we were interested whether age-associated changes in miRNA expression might be related to the proteome data and transcriptome data listed above. Microarray-based miRNA profiling revealed 12 miRNAs showing evidence of age-related differential expression (Tab. 2) [11]. To determine if these miRNAs play a role in mRNA abundance, we compared their targets with the mRNAs from Kalfalah and colleagues [10]. This comparison revealed 164 mRNAs as candidate targets of the 12 candidate miRNAs (Fig. 3; Tab. S8). However, none of the targeted mRNAs were found to be significantly altered with age. Moreover 59% (96 mRNAs) exhibited a constant expression during aging. At the protein level we detected three proteins with age-associated alteration that possibly may be explained by miRNA mediated regulation (Tab. 2; Fig. 4). DNAJC3, DDX3X and CALM1 exhibited an anti-correlation with miR-107, miR-29c as well as with miR181a and miR-409-3p, respectively.

Age-associated biological processes

With the integration of data from proteome, transcriptome and miRNA analysis we assessed a broad spectrum of candidate genes/proteins involved in aging processes. Next, we were interested to reveal the biological processes which allow us to predict aging processes at the functional level. Kalfalah and colleagues reported on 117 age-related biological processes based on the transcriptome data (Tab. S9) [10]. From the proteome data, we deduced for the 43 age-associated proteins 71 unique biological processes as being significantly enriched (p ≤ 0.01; Tab. S10). We found eight main categories of GO terms: ‘proteostasis’, ‘response to stimuli and stress’, ‘development and differentiation’, ‘cell organization’, ‘cell communica-tion and signal transduction’, ‘RNA metabolism’, ‘cell cycle and proliferation’ and ‘cell death’ (Fig. 5). Remarkably, comparison of biological processes related to the identified age-associated mRNA and protein candidates revealed ‘cell cycle’, ‘development’, ‘transcription’, ‘translation’, ‘regulation of actin cytoskeleton’, ‘proteasome’, ‘RNA and protein metabolism’, ‘signaling’, ‘RNA biosynthesis’, ‘oxidative phosphorylation’, ‘respiratory electron transport’ and ‘signal transduction’ as highly enriched in both transcriptome and proteome data. Interestingly, we revealed most of these processes by the enrichment analysis of 164 target mRNAs of the 12 identified age-associated candidate miRNAs (Tab. S11). Although not confirmed at the protein level or by mRNA alteration, enrichment analysis showed that identical age-associated biological processes such as ‘signal transduction’, ‘cell communication’, ‘cell death’, ‘RNA and protein metabolism’, ‘cell communication’, ‘response to stimuli’, ‘development’, ‘proliferation’, ‘cellular component organization’ and ‘RNA and protein biosynthesis’ were highly enriched. The comparisons between different data sets by the deduced biological processes resulted in a significant overlap and led to a biological processes-oriented picture of aging in fibroblasts not expected from the single gene/protein level.

Age-associated loss of proliferative potential

The limited proliferating potential of in vitro cultured cells has been widely investigated in the context of cellular aging [4]. We identified six proteins known to be involved in cell cycle regulation as showing age-associated expression changes in dermal fibroblasts. These included reduced expression of the two subunits TCP1 and CCT6A of the T-complex protein 1 that plays an important role in the transition from S to M phase [17]. Another protein with age-associated decreased expression was CAPNS1 which is the regulatory subunit of the calcium-regulated non-lysosomal thiol-protease. Knock-down of CAPNS1 inhibited proliferation and cell growth [18]. Similarly, DNAJC3, an inhibitor of the kinase PKR, showed reduced expression in aged fibroblasts. Reduced expression of DNAJC3 resulted in increased activity of PKR and phosphorylation of EIF2A, whereby protein translation and cell growth is down-regulated [19]. The protein DDX3X, a multifunctional ATP-dependent RNA helicase involved in several steps of gene expression, also demonstrated lower expression in fibroblasts from aged versus young donors. DDX3X is important for translation of Cyclin E1 mRNA and thereby may stimulate cell cycle progression [20]. In contrast THBS1 protein abundance increased with age in cultured fibroblasts. Previous data indicate that THBS1 may cause cell cycle arrest of endothelial cells by increasing the amounts of p21 and unphosphorylated RB1 [21]. Collectively, our findings on aging-associated expression of cell cycle regulatory proteins are in line with a role of reduced proliferation capacity of aged fibroblast discussed in skin aging and provide candidate proteins for further functional validation of this concept.

Age-associated features of geroconversion

Cell isolation and culture

Cell isolation and culture were performed as previously described [10]

Sample preparation for MS-analysis

Human dermal fibroblasts were cultivated under standardized conditions until passage three as described above. Cells were harvested, homogenized and lysed in urea-buffer (2 M thiourea, 7 M urea, 30 mM Tris-HCL, pH 8.0). Subsequently, total protein amount was determined using BCA protein assay (Thermo Scientific Pierce, Rockford, USA) and 10 μg of each sample were proteolytic digested with trypsin (Promega, Mannheim, Germany). Here, the proteins were incubated with 10 mM dithiothreitol at 56 °C for 45 min in Amicon Ultra filters (3 kDa cut-off; Millipore, Billerica, USA). Subsequently, 0.55 M iodoacetamide was added to the mixture and incubated for 30 min at room temperature and light protected. The mixture was centrifuged for 30 min at 14,000 x g and 4 °C. Next, 50 mM ammonium hydrogen carbonate was added and again centrifuged for 60 min at 14,000 x g and 4 °C. Proteins were recovered and 1:50 trypsin was added. Digestion was performed at 37 °C for 12 h. Afterwards, 0.1% trifluoroactic acid (TFA) was added to stop the digestion.

Mass spectrometric analysis

Mass spectrometric analyses of proteomes were carried out using highly reproducible and stable LC-MS/MS system and a label-free approach for quantification. Peptides of each sample were analyzed with a nano-HPLC/ESI-MS system composed of an RSLCnano HPLC and a LTQ Orbitrap Velos (Thermo Fisher Scientific, Bremen, Germany) mass spectrometer equipped with a nano-electrospray ion source. Each sample was loaded onto a trap C₁₈ trapping column (2 cm × 100 μm × 5μm, 100 Å, Thermo Fisher Scientific) and desalted with 0.1% TFA for 10 minutes. Peptides were eluted from the trap, separated with an analytical column (Acclaim PepMap RSLC C₁₈; 25 cm × 75 μm × 2 μm, 100 Å, Thermo Fisher Scientific) for 120 minutes with a flow rate of 300 nL/min, and sprayed into the MS. The Orbitrap parameters were as follows: spray voltage, 1.4 kV; ion transfer tube temperature, 275 °C; collision gas, helium; collision gas pressure, 1.3 mTorr; normalized collision energy for MS/MS, 35%. MS-spectra were acquired in the Orbitrap with a mass range of 300-1500 m/z and a resolution of 60,000. The Orbitrap was operated in a TOP4 data-dependent mode to automatically switch between MS and MS/MS acquisition. MS/MS-spectra were acquired applying gas-phase fractionation. Therefore, five m/z areas were used for precursor selection: 300-450 m/z, 450-550 m/z, 550-650 m/z, 650-800 m/z and 800-1500 m/z (for more information see Fig. S6). Polysiloxane (445.120030 Th) was used as lock mass. For MS/MS, ions were isolated with an isolation width of 2 m/z and fragmented using collision-induced dissociation. MS/MS-spectra were acquired in the linear ion trap in centroid mode. Target ions selected for MS/MS were dynamically excluded for 45 sec. The ion selection thresholds were 500 counts for MS/MS. An activation Q of 0.25 and activation time of 10 ms was applied in MS/MS acquisitions.

Protein identification and quantification

For protein identification, the Proteome Discoverer (version 1.3, Thermo Fisher Scientific) and MASCOT search engine (version 2.4.1, Matrix Science, London, UK) were used. MS/MS-spectra were searched against the UniProtKB database (date 04/03/2013). Search parameters were as follows: enzyme, trypsin; missed cleavage sites, 2; taxonomy, homo sapiens; precursor mass tolerance, 10 ppm; fragment mass tolerance, 0.4 Da; oxidation of methionine as dynamic modification; carbamido-methylation of cysteine as fixed modification. The false discovery rate was set to 1% (p ≤ 0.01). Label-free quantification of proteins was carried out using Progenesis LC-MS (Nonlinear Dynamics, Newcastle upon Tyne, UK) as describes previously [50]. A minimum of 2 unique peptides were required for quantification. Only proteins identified in 14 of 15 samples were regarded for quantification. Statistical analysis for the proteomes was performed using two approaches: multiple t-tests (ANOVA) and Pearson correlation. For ANOVA analysis a significance threshold of 5% (p ≤ 0.05) was applied. Pearson correlation was performed without any grouping with a significance threshold of 5% (p ≤ 0.05). Stability of protein and mRNA abundance was determined as described elsewhere [51]. Therefore, protein/mRNA intensities were logarithmized (log10) and for each protein/mRNA the absolute mean difference between the age-groups (‘mean young’ – ‘mean middle’; ‘mean young’ – ‘mean old’; ‘mean middle’ – ‘mean old’) was calculated. Afterwards, the ‘two one-sided test’ for equivalence (TOST) was applied to test for similarity of the groups (e.g. H₁: |‘mean young’ – ‘mean middle’| < equivalence range (ɛ)) [52] [51] [50] [50] [48]. As equivalence range, the three-fold standard deviation (ɛ = 0.3) was used. Proteins/mRNAs whose p-value was significant (p ≤ 0.05) in all of the three pairwise comparisons were marked as unchanged.

Functional annotation, network and enrichment analysis of proteins

Gene set analysis of unchanged proteins was carried out using Consensus Path DB [53]. The entire identified proteome was used as background list. Enrichment was applied on Gene Ontology, Reactome, Wikipathways and Pathway Interaction Database biological processes. We discarded gene sets that were redundant, had ≤ 5 members or a p-value above 0.01. Enrichment analyses of abundance classes and mitochondria related proteins were carried out using DAVID [54] and the Gene Ontology Biological Process categorization. We discarded processes that had a p-value above 0.05. Network and enrichment analysis was carried out using Cytoscape environment [55] and ClueGo plug-in [56]. Parameters for network analyses were applied as follows: Ontology source, Gene Ontology biological processes or cellular components (date 12/10/2013); statistical test, enrichment/depletion (two-sided hypergeometric test); p-value correction, Benjamini-Hochberg; p-value restriction, p ≤ 0.01; network specificity, ‘medium’; GO term restriction ‘min level = 4′, ‘max level = 8′, min percentage = 4.0; use GO term fusion; use GO term grouping; GO term connection restriction, kappa score ≥ 0.5. The entire identified proteome was used as reference set.

mRNA/protein correlation analysis

Correlation analysis between measured mRNA-expression levels and protein abundance was performed using R ( http://www.r-project.org). The mRNA intensities where log2 scaled, the lowest 25% based on the median intensities as well as probes with no corresponding EntrezID where discarded. In the case of multiple probes corresponding to a single Entrez GeneID the mean intensity was used for further analysis. The maximum mRNA intensity was set to 100%. For each identified protein the sum intensities were divided by the underlying number of quantified peptides. The protein data was log2 scaled and the maximum value was set to 100%.

miRNA network and enrichment analysis

The age-associated miRNAs reported by Röck et al. [11] were used to build up a miRNA-gene regulatory network as follows: We built up an initial network for these miRNAs based on target predictions, which were derived from TargetScan [57]. We then integrated gene expression measurements in order to find relationships between miRNAs and genes, which can be deduced by the data. We fitted a linear model with the gene expression data as response and the expression data of its predicted targeting miRNAs as predictor variables. We then performed a multiple linear regression analysis for each gene with elastic net penalty in order to select for miRNAs that have an influence on the gene expression data. For this purpose, we used the glmnet package for the R [58] and introduced a negativity constraint to select only for negative miRNA effects. The penalty parameter was determined by 10-fold cross-validation. We finally obtained two networks by applying this procedure to the mRNA and the protein expression data, respectively.