Similarities in Gene Expression Profiles during In Vitro Aging of Primary Human Embryonic Lung and Foreskin Fibroblasts

Replicative senescence is of fundamental importance for the process of cellular aging, since it is a property of most of our somatic cells. Here, we elucidated this process by comparing gene expression changes, measured by RNA-seq, in fibroblasts originating from two different tissues, embryonic lung (MRC-5) and foreskin (HFF), at five different time points during their transition into senescence. Although the expression patterns of both fibroblast cell lines can be clearly distinguished, the similar differential expression of an ensemble of genes was found to correlate well with their transition into senescence, with only a minority of genes being cell line specific. Clustering-based approaches further revealed common signatures between the cell lines. Investigation of the mRNA expression levels at various time points during the lifespan of either of the fibroblasts resulted in a number of monotonically up- and downregulated genes which clearly showed a novel strong link to aging and senescence related processes which might be functional. In terms of expression profiles of differentially expressed genes with age, common genes identified here have the potential to rule the transition into senescence of embryonic lung and foreskin fibroblasts irrespective of their different cellular origin.


Introduction
Cellular senescence is a terminal phase observed towards the end of a primary human fibroblast cell population after numerous cell divisions; it is considered to be the cellular aging process. Cellular senescence occurs either naturally or stress induced; that is, cells stop dividing after a finite number of cell divisions (termed "replicative senescence"), reaching the final cell cycle arrested state called the "Hayflick limit" [1]. The process of senescence is associated with a number of phenotypes; in general, the integrity and function of tissues decline, resulting in the body being susceptible to diseases associated with age [2,3]. Key factors driving cellular senescence are induced increase in Cyclin dependent kinase inhibitors (CDKIs) [4], oxidative stress [5], and DNA damage [6,7]. In senescence, despite their viability and active metabolism, cells are resistant to mitogenic or apoptotic stimuli [8,9]. On the one hand, cellular senescence results in irreversible growth arrest, limiting the proliferation of damaged cells susceptible to neoplastic transformation resulting in a decreased incidence of cancer. However, on the other hand, senescence results in in vivo aging, weakening the function and renewal of stem cells [10]. Markers are able to identify cellular senescence in vitro and in vivo: enlarged cell morphology, increase in amount of cellular debris, changes in chromatin structure, increase in Cyclin dependent kinase inhibitors (CDKIs) expression, presence of senescence associated secretory phenotype (SASP), and senescence associated ß-galactosidase (SA ß-Gal) [11][12][13]. DNA damage response and the p53-p21 and p16-pRb pathways are crucial for senescence induction [14], together with additional pathways including telomere uncapping, DNA damage (UV, ionizing radiation, and chemicals), cytoskeletal genes, the interferon pathway, nutrient imbalances, oncogenic activities, and oxidative stress [7,15,16]. In primates, 2 BioMed Research International the percentage of senescent skin fibroblasts increases with age in vivo [17]. Here, we therefore used primary human fibroblasts [9,13,18] as our model system.
Recently, we identified individual gene expression patterns during replicative senescence among five fibroblast cell lines of different cell origins [13]. In this study, we determined mRNA expression changes during different stages of their lifespan in two fibroblast cell lines of different cell origin. We analyzed the transcriptome, determined by RNAseq, at five separate population doublings (PDs) between young and senescent embryonic lung (MRC-5) and foreskin (HFF) fibroblasts. Using both molecular and systems biology approach, we studied the growth pattern of the two fibroblast cell lines in detail. By comparing fibroblasts from two different origins we were able to determine either mRNA changes specific for one of the cell lines or common transcriptomic patterns which underlie the process of replicative senescence.

Cell
Culture. Primary human fibroblast cells were cultured in Dulbeccos Modified Eagle's Low Glucose Medium (DMEM) with L-glutamine (PAA Laboratories, Pasching, Austria), supplemented with 10% fetal bovine serum (FBS) (PAA Laboratories) under normal air conditions in a 9.5% CO 2 atmosphere at 37 ∘ C. The cells were subcultured by removing the remaining medium followed by washing in 1x PBS (pH 7.4) (PAA Laboratories) and detachment using trypsin/EDTA (PAA Laboratories). Primary fibroblasts were subcultured in a 1 : 4 (= 2 population doublings (PDs)) or 1 : 2 (= 1 PD) ratio. For stock purposes, cryoconservation of the cell lines at various PDs was undertaken in cryoconserving medium (DMEM + 10% FBS + 5% DMSO). Cells were immediately frozen at −80 ∘ C and stored for two to three days. Afterwards, cells were transferred to liquid nitrogen for long time storage. Refreezing and rethawing was not performed to avoid premature senescence [20].
A vial of each of the two fibroblast cell lines (MRC-5 and HFF) was obtained and maintained in culture from an early PD. On obtaining enough stock on confluent growth of the fibroblasts in 75 cm 2 flasks, cells were subcultured into three separate 75 cm 2 flasks and were passaged until they were senescent in culture. At five different time points of the fibroblast's span in culture (MRC-5 = PDs 32, 42, 52, 62, and 72 and HFFs = PDs 16, 26, 46, 64, and 74), the total RNA was extracted and used for high-throughput sequencing.

Detection of Senescence
Associated ß-Galactosidase (SA ß-Gal). The SA -Gal assay was performed as described by [11] at each of the five PDs in both MRC-5 and HFF. Paired twosample type 2 Student's -tests assuming equal variances were done to examine the values obtained from SA ß-Gal assay for statistical significance [9].

RNA Extraction.
Total RNA was isolated using Qiazol (Qiagen) according to the manufacturer's protocol, with modifications as explained in [9].
2.6. Quantitative Real-Time PCR. Real-time PCR was performed using CFX384 thermocycler Biorad and Quantitect PCR system (Qiagen) as described earlier in [23]. Three reference genes (GAPDH, ACTB, and RAB10) were used for   [26] and the respective GTF gene annotation, obtained from the Ensembl database [27]. Gene counts were further processed using the R programming language [28] and normalized to RPKM values. RPKM values were computed using exon lengths provided by feature-Counts and the sum of all mapped reads per sample.

Sample Clustering and Analysis of Variance.
Spearman correlation between all samples was computed in order to examine the variance and the relationship of global gene expression across the samples, using genes with raw counts larger than zero. Correlation values were visualized using a heatmap ( Figure 1). Additionally, principal component analysis (PCA) was applied using the log 2 RPKM values for genes with raw counts larger than zero. Results were visualized in a three-dimensional scatterplot ( Figure 2 [29] and edgeR 3.4.2 [30] were used to identify differentially expressed genes. Both packages provide statistics for determining of differential expression in digital gene expression data using a model based on the negative binomial distribution. The nonnormalized gene counts have been used here, since both packages include internal normalization procedures. The resulting values were adjusted using Benjamini and Hochberg's approach for controlling the false discovery rate (FDR) [31]. Genes with an adjusted value < 0.05 found by both packages were assigned as differentially expressed. Since large sets of DEG were found more strict selection cutoffs have been used: adjusted value < 0.01 (by both packages) and absolute log 2 fold-changes > 1.

Comparison of RNA-seq with qRT-PCR and Protein
Expression. Correlation analysis was performed using all 15 samples (3 replicates for each of the 5 PDs) for MRC-5 and HFF, respectively. Spearman correlation coefficients were estimated using the RPKM values (RNA-seq data) and 2 −ΔCT values (qRT-PCR data). For comparison of RNA-seq with Western Blot data, only the first and the last PD were used. Log 2 fold-changes were calculated based on RPKM values (RNA-seq data), 2 −ΔΔCT ratios (rRT-PCR), and protein expression ratios (Western Blots).

Clustering of Expression Profiles.
Genes were clustered according to their temporal profiles using a fuzzy -means algorithm. We used the function cmeans from the package e1071 1.6-2 of the R programming language. Parameters were defined as = 1.2, iter.max = 500, d.obj.fun = 10 −8 . The number of trials for the fuzzy algorithm was set to 30. The optimum number of clusters was determined using a combination of several cluster validation indexes as described by [32]. See Supplemental Table 2 for detailed assignment of the genes to clusters.

Functional Enrichment
Analysis. Singular gene set enrichment analysis was performed using FungiFun2 [33] for selected sets of genes based on the clustering results. Although FungiFun2 is mainly suited for fungal gene enrichment analysis annotation for human genes is included as well and was recently updated. Default parameters were used while significant Gene Ontology (GO) terms and KEGG pathways were selected according to FDR corrected values < 0.05. Complete lists of GO-terms and KEGG pathways are available from Supplemental Table 3. The list of GO-terms was further summarized using TreeMaps of the REVIGO online tool [34]. Default parameters and GO term with adjusted values were used as input.
2.14. Monotonically Expressed Genes. In order to identify genes that change their expression levels monotonically with age, we calculated the Spearman correlation coefficient ( ) of each gene's temporal profile with the linearly increasing curve ( ) = . In order to incorporate the replicates at each time point, we repeated the calculations by randomly sampling over the replicates at each time point and by calculating an average correlation coefficient from the resampled curves afterwards. values ( ) were computed using the R base function cor.test. We used the calculated correlation coefficient ( ) of gene with the linear increasing curve as a criterion to split the genes into the following three groups: if ( ) > 0 and ( ) < 0.05, we considered a gene to be monotonically increasing with age, if ( ) < 0 and ( ) < 0.05, the gene was considered to be monotonically decreasing with age, and if ( ) > 0.05, the expression of the corresponding gene was considered nonuniformly [35]. See Supplemental Table 4 for detailed test results.

Functional Association Networks.
Gene symbols were used as input for functional association network creation using the STRING database [36], Cognoscente [37], and Gen-eMANIA [38] online tools. For STRING we used "multiple names" input and selected "Homo sapiens" as organism. In Cognoscente we selected "Human" and "Radius = 0 with intermediates" as input parameters in addition to the list of genes. The GeneMANIA network was created using default settings. The resulting networks are shown in Figure 9 and Supplemental Figure 6.

Results and Discussion
We studied the growth of two primary human fibroblast cell lines, MRC-5 and HFF, throughout their span in culture from an early PD until they achieved senescence at late PDs. Analysis of their growth behaviour (Supplemental Figure 1A) and their entry into senescence (Supplemental Figure 1B), measured by the induction of SA -Gal, revealed a cell line specific transition into senescence of these two fibroblasts (MRC-5 derived from embryonic lung and HFFs derived from human foreskin). Fibroblast cell line specific growth has been observed by us before [13,39]. Total RNA was extracted at five different time points of the fibroblasts span in culture and was subjected to high-throughput RNA sequencing (RNA-seq).

Global Expression Profiles Cluster according to Cell Line and Age.
Overall, the RNA-seq data of this study comprise 30 samples: 15 samples for each cell line (HFF and MRC-5), consisting of five different PDs, each with three biological replicates. For each sample, mapping and counting resulted in 56,299 raw gene count expression values (using Ensembl gene annotation). The largest group of these genes (21,226) belongs to the group of protein coding genes. For all the 30 samples, 19,237 genes have raw gene counts larger than zero; these genes were considered for further analysis.
First, we studied primary clustering of the global gene expression. We therefore created a heatmap showing the Spearman correlation for all 30 samples using the nonzero genes ( Figure 1). In this heatmap, both cell lines were clearly separated. In eight out of ten cases, the three values of the replicates were clustered together, showing the good quality of the data and low noise between the replicates. Next, we applied principal component analysis (PCA) to further investigate the effect of aging in the individual cell lines. Figure 2 shows the first three principal components which explain ∼97% of the variances in a three-dimensional plot. MRC-5 and HFF again were clearly separated (by PC2) and the effect of aging was covered by PC1 and PC3, with a larger separation between young and old HFF compared to MRC-5. Already at this global level, similarity between both cell lines is perceptible, since young and senescent samples are grouped concordantly.

MRC-5 and HFF Share Common Differentially Expressed
Genes Regulated by Aging. Differentially expressed genes (DEG) were identified by comparing all consecutive PDs as well as the first with the last PD in MRC-5 and HFF cells (10 comparisons; Table 1 . Figure 3 reveals that DEG were not specific for a certain PD comparison but recurred when later PDs were compared. MRC-5 and HFF shared a large fraction of DEG; only a minor fraction of DEG was identified uniquely in one of the cell lines (bars on the right in Figure 3). This indicates common processes which occur during aging in both cell lines rather than cell line specific changes. Most DEG were found when comparing the first with the last PD, leading to new DEG which had not previously been detected between consecutive PDs (orange and turquoise coloured bars in Figure 3). Both cell lines differed between the absolute number of DEG as well as the increased percentage of DEG for the first two transitions in HFF (PD 16 to PD 26; PD 26 to PD 46). In MRC-5, most of the changes seemed to occur at late PDs, while HFF cells indicated larger changes already after early PDs. This effect is also perceivable by the distances in the PCA plot ( Figure 2).

High Correlation of RNA-seq with qRT-PCR for Selected
DEG. For validation of the RNA-seq data, qRT-PCR was applied. Here, triplicates for all five PDs were measured.
Selection of genes was based on the comparison of the first with the last PD, using the strict DEG criteria ( < 0.01 and |log 2 fold-change| > 1), resulting in 2,117 DEG for MRC-5 and 4,651 DEG for HFF (5th and 10th bars in Figures 3(b) and 3(d)). We further filtered the intersection of those two gene sets (1,139)   showing again the similarity of gene expression changes in both cell lines. Overall 12 DEG were selected which either showed strong expression in the RNA-seq data (RPKM > 50; genes: EGR1, CCND1, CTSK, DKK3, IGFBP7, and TMEM47) or were proven to have an established role in cell cycle and senescence pathways (CCNA2, CCNB1, ID3, IGFBP2, MMP3, and WNT16). The minimal RPKM criterion was applied to ensure a strong expression signal in at least one condition for a set of genes. The expression profiles from both measurement techniques were then confronted using Spearman rank correlation in each individual cell line. The results showed high correlation coefficients indicating a good overlap of both measurement techniques and quality of highthroughput gene expression analysis ( Figure 4). In 17 out of 24 cases, correlation was larger than 50% (mean correlation of 63%), while only once negative correlation was found (EGR1 in MRC-5).

Consistent Changes of mRNA and Protein Expression.
Although mRNA expression changes are generally considered to consequently lead to corresponding changes in protein levels, correlation between both can be as little as 40%, as observed in large-scale proteome-and transcriptomeprofiling experiments [40]. We thus asked if the detected changes of strongly altered DEG correlated with corresponding protein expression levels. Triplicates of the first and last PD were selected for comparison. Gene selection was performed as described above, either by strong expression in RNA-seq (RPKM > 35) or by functional relation to cell cycle and senescence pathways. Overall, 28 DEG (16 downand 12 upregulated DEG) were selected (the genes mentioned above, validated by qRT-PCR, were included in this set). The results of this comparison showed consistent changes, in terms of their direction of regulation, between mRNA expression, measured by RNA-seq, and protein expression, measured by Western Blots, for all selected genes ( Figures  5 and 6). 44 out of the 56 protein fold-changes exhibited significant differences between young and old PDs.

Common Genes Ruling the Transition into Senescence in MRC-5 and HFF.
Then, we asked if common cellular markers are involved in the transition into senescence. We thus studied the genes most differentially expressed with age commonly in MRC-5 and HFF fibroblasts. We noticed that a large number of genes among the most differentially expressed genes belonged to the secretory phenotype (Figure 8(a), as explained in Section 3.7). The list of genes included CTSK, normally stimulated by inflammatory cytokines released after tissue injury [41], GRN, a previously functionally validated gene responsible for wound healing [42], CST3, associated with sarcopenia [43], and PERP, a p53 apoptosis effector, the mRNA expression level of which is upregulated in human mesenchymal stem cells [44]. We detected significant upregulation of IGFBP2 which was found upregulated with senescence in retinal pigment epithelial cells [45,46] and BJ fibroblasts [47]. ID1, ID3, CCNA2, and CCNB1 showed significant downregulation with age in our study for both human fibroblasts. Downregulation of ID1 and ID3 expression with senescence was detected in BJ foreskin, WS1 fetal skin, and LF1 lung human fibroblasts [48] and of CCNA2 in IMR-90 and WI-38 [49]. Targeting CCNB1 expression inhibits proliferation of breast cancer cells [50]. The list of most differentially expressed genes also included IGFBP7 and MMP3 which encode protein receptors predominantly located on the cell surface. Both IGFBP7 and MMP3 are upregulated with senescence in human melanocytes [51][52][53][54]. Recently we found that overexpression of recombinant IGFBP7 proteins induced premature senescence in early PD MRC-5 fibroblasts [13]. Among the genes significantly upregulated with age in both MRC-5 and HFFs we identified DKK3, having a role in Wnt signaling [55][56][57]. DKK3 has tumor suppressor activity in breast cancer patients [58] and in papillary thyroid carcinoma [57]. However, we had failed to demonstrate an induction of premature senescence in early PD HFFs on overexpression of recombinant DKK3 proteins [13]. Though not significantly differentially expressed with age in MRC-5 fibroblasts, one of the genes which were most significantly upregulated with age in HFFs was the SFRP4 gene, an antagonist for Wnt signalling [59]. SFRP4 acts as a tumor suppressor in gastric carcinoma [60] and epithelial ovarian cancer cell lines [61]. In a separate study, we functionally validated the expression of SFRP4 in early PD HFF and MRC-5 fibroblasts by treating them separately with human recombinant SFRP4 protein. This treatment resulted in premature senescence induction in HFFs but not in early PD MRC-5 fibroblasts [13].
Here, induction of SFRP4 mRNA expression was not detected by RNA-seq, explaining the lack of premature senescence induction in early PD MRC-5 fibroblasts. SFRP4 expression thus showed cell line specific differences.

Clustering of the Expression Profiles Shows Similar Pattern
in Both Cell Lines. We found many differentially expressed HFF RNA-seq HFF qRT-PCR HFF protein * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * genes commonly regulated in both MRC-5 and HFF. Next, we asked if both cell lines exhibit common temporal expression profiles rather than showing different effects for the same set of genes. Therefore, we applied fuzzy -means clustering comparing the expression profiles of both cell lines. We used 1,803 genes found to be differentially regulated between the four consecutive PDs and between the first and the last PD in both cell lines, according to the strict cutoffs as shown in grouped in clusters 1 and 4. Interestingly, most genes follow a monotonic profile (either up or down) while only few genes exhibit a parabolic-like shape. Clusters 1 and 3 show major changes in MRC-5 between the second to last and the last PD, while cluster 4 groups genes with large differences between the first and the second PD in HFF. This effect was already observed when comparing DEG between the consecutive PDs (see DEG section above).

Identification of Functional Categories Significantly Enriched for Genes with Common Expression
Profiles. Next, we deduced the main biological processes driven by the differentially expressed gene sets obtained from the cluster analysis. Using gene set enrichment analysis, for each of the five clusters significant GO categories and KEGG pathways could be identified. The results indicated a strong connection of upregulated genes (grouped in clusters 3 and 5) to "extracellular space" (GO:0005615) and "membrane" (GO:0016020) components (Figure 8(a)). Corresponding KEGG pathways, found for these genes, were for example, "ECM-receptor interaction" (hsa04512) and "ABC transporters" (hsa02010; see Supplemental Table 3). ABC proteins transport various molecules across extra-and intracellular membranes and are involved in aging and agerelated diseases [62]. Cluster 3 shows stronger upregulation at late PD for MRC-5 cells while HFF cells are upregulated more clearly in cluster 5 (see above). Comparing the GO-terms found for these single clusters, we found links of the stronger upregulation in MRC-5 with "integral component of plasma membrane" (GO:0005887) and the "Golgi apparatus" (GO:0005794), while "sarcolemma" (GO:0042383) and "nucleosome" (GO:0000786) were more specific for upregulation in HFF (Supplemental Figure 2). The structure of the secretion regulating Golgi complex is altered in senescent cells [63]. While our results indicate cell line specific differences during replicative senescence, the GO-term comparison revealed that in both cell lines many genes were similarly upregulated. A large set of GOterms associated with upregulated genes were related to the senescence associated secretory phenotype [64].
Downregulated genes (grouped in clusters 2 and 4) were associated to strongly enriched GO processes related to, for example, "cell cycle" (GO:0007049), "cell division" (GO:0051301) and "DNA replication" (GO:0006260) (Figure 8(b)). Here, differences between both cell lines are more obvious. After a slight initial gain, expression in MRC-5 cells declined strongly between PD 46 and PD 64. In HFFs, strong decline already started after PD32 without larger changes for late PD (Figure 7; cluster 4). Most of these cell cycle related genes, which account for the above-mentioned profiles, are related to the cellular component "nucleoplasm" (GO:0005654). Associated GO-terms for cluster 2, which depicts moderate downregulation, were more widespread and covered processes like "positive regulation of nitric oxide biosynthetic process" (GO:0045429), "endoderm formation" (GO:0001706), and "response to cAMP" (GO:0051591; Supplemental Figure 3).
Cluster 1 showed the largest differences between both cell lines. While, in MRC-5, genes are downregulated strongly at the last PD, no clear up-or downregulation is observed  for HFF. Significantly enriched GO-terms associated to these genes were, for example, "vasculogenesis" (GO:0001570), "response to lipopolysaccharide" (GO:0032496), and "cell adhesion" (GO:0007155) (Supplemental Figure 4).

Monotonically Regulated Genes in MRC-5 and HFF Are Connected in Functional Association
Networks. Since senescence is a continuous cellular process, it can be hypothesized that genes possessing key relevance for senescence change their expression values monotonically over time, while genes with irregular temporal expression patterns might be associated with response to environmental conditions, with the circadian rhythm or other processes. Amongst others, continuous increasing and decreasing profiles were found by the clustering analysis. In addition to this nonbiased approach, we intended to identify genes with a strong monotonic behaviour across the PDs investigated here. We calculated the Spearman correlation coefficient of each gene's temporal profile with a linearly increasing sequence. Replicates for each PD were incorporated by a random sampling approach. Subsequently, we classified genes into three classes according to their behaviour with age: (a) monotonically upregulated genes, (b) monotonically downregulated genes, and (c) nonuniformly regulated genes ( Table 2).
More monotonically up-and downregulated genes were found for HFFs compared to MRC-5 (888 versus 179). Only a small subset of these genes were commonly regulated in both cell lines (9 up and 14 down) but even less genes showed an opposite monotonic expression profiles (8; see Supplemental Figure 5 and Supplemental Table 4). The 23 commonly monotonically up-or downregulated genes were studied in more detail. Since in both cell lines the regulation of these genes strongly correlated with an increase of senescence, they might play an essential role in cellular aging and may rule common regulatory process. We used several online resources in order to find potential or validated interactions between these genes. The STRING database [36] only provides the interactions between four out of all the 23 genes (Supplemental Figure 6A). Using Cognoscente [37], 17 out of 23 genes were connected within one interaction graph (Supplemental Figure 6B). More interactions could be found using GeneMANIA [38], leading to a network which is widely connected by coexpression and common pathways like, for example, "epithelial cell proliferation" and "extracellular matrix organization" (Figure 9). Both of the latter tools integrate intermediate genes which were not in the input list. Hub genes in these networks included ATF7, MAF, UBC, and ELAVL which are interesting candidates for further studies. All the four of these genes were functionally associated with tumorigenesis. Members of the ubiquitin family including UBC have been associated with tumor progression [65]. In terms of ATF7, the activating transcription factor family is associated with cell proliferation and oncogenesis [66]. Both MAF and ELAV1 have been associated with oncogenesis and tumor progression [67,68]. Thus all the four genes had an association with cell proliferation. Then, we investigated the biological relevance of the monotonically up-and downregulated genes in both fibroblast cell lines. The list of monotonically downregulated genes included NLE1, AMMECR1, FIBCD1, ENPP2, TMTC4, ANPEP, MYC, EFNB3, HCLS1, FERMT1, FABP5, SPHK1, GOS2, and RPL36A. The genes monotonically upregulated included LRP10, TMCO3, CAV2, ADAMTS5, C5orf15, SDC2, ANKH, PCDHB16, and TGFB2. A number of genes in the above list have been functionally associated with proliferation.

Monotonically Downregulated
Genes. NLE plays a role in regulating the Notch activity and is involved in embryonic development in mammals by affecting the CDKN1A and Wnt pathways [69]. Forced expression of miR-26 inhibits the growth of stimulated breast cancer cells and tumor in xenograft models by reducing the mRNA expression levels of AMMECR1 and other genes [70]. AMMECR1 is associated with Alport syndrome, mental retardation, midface hypoplasia, and elliptocytosis [71]. FIBCD1 (fibrinogen C domain containing 1) binds to chitin of invading parasites [72]. FIBCD1 is primarily present in the gastrointestinal tract of humans; however, their presence in skin has been highly debated [73,74]. ENPP2 facilitates cell motility and progression and is related to the invasion of ductal breast carcinomas [75]. TMTC4 is a gene contributing to embryonic brain development; it interacts with Wntless, an integral Wnt regulator [76]. EFNB3, a member of the ephrin gene family, is associated with neural development [77]. ANPEP is a well-known marker for acute myeloid leukemia and tumor invasion; it has a regulatory role in angiogenesis [78,79]. FERMT1 is overexpressed in colon and lung carcinomas [80]. The MYC oncogene is associated with cell growth regulation by driving proliferation via upregulation of Cyclins and downregulation of p21 [81,82]. HCLS1 gene which is monotonically downregulated with age is associated with antigen receptor signaling and clonal expansion as well as deletion of lymphoid cells [83]. The FABP5 gene encodes the fatty acid binding protein in epidermal cells and is upregulated in psoriatic tissues [84]. SPHK1 has been previously associated with melanoma progression and angiogenesis [85,86]. The GOS2 gene promotes apoptosis by binding to BCL2, hence preventing the formation of protective BCL2-BAX; its mRNA and protein levels are downregulated in type 2 diabetic patients [87,88]. Thus, almost all genes, monotonically downregulated with age in both fibroblast cell lines, are associated with proliferation and cell survival.  Figure 9: Functional association network by GeneMANIA for 23 monotonically up-or downregulated genes. Input genes are depicted by either black, red, or blue color or large circles. Intermediate genes, added by the tool, are shown using small grey circles. Colored genes denote functions associated by the tools: blue = "epithelial cell proliferation"; red = " extracellular matrix organization." The edge colors indicate the type of interaction, as explained in the legend on the bottom right.

Monotonically Upregulated
Genes. LRP10, a negative regulator of Wnt signalling, was found monotonically upregulated with age [89]. CAV2 is a scaffolding protein within the caveolar membrane modulating cancer progression [90]. ADAMTS5 enables the destruction of aggrecan in patients with arthritic disease which is prevalent with aging [91]. The ANKH gene, associated with regulation of tissue calcification and in turn susceptibility to arthritis, is also monotonically upregulated with age in both fibroblast cell lines [92]. Syndecan-2 protein (SDC2) is upregulated in skin and lung tissues of patients suffering from (age-associated) systemic sclerosis and fibrosis [93,94]. mRNA expression of PCDHB16 is upregulated in patients with (age-associated) Alzheimer's disease [95]. TGFB2, also monotonically upregulated with age in fibroblasts, has suppressive effects on interleukin-2 dependent T cell proliferation and displays effector functions [96].
In summary, the genes, which we found here monotonically up-and downregulated with age in both fibroblast cell lines, have been studied before. In this study, we explicitly show for the first time the age-associated regulation of these genes in primary human fibroblast cells of two different origins. In a following study we will determine the protein expression of all age-related genes and functionally validate the expression of these genes.

Conclusion
We studied molecular aspects of cellular aging by determining the differential expression of genes during the aging of two primary human fibroblasts, MRC-5 and HFFs. RNAseq data analysis encompassed different levels, starting from the complete set of annotated and expressed genes, proceeding to different gene subsets and functional categories. Most of the detected changes were found to be common in both cell lines, as indicated by the large number of overlapping DEG and common expression profiles identified by clustering. We validated the expression patterns for selected genes, demonstrating an association of almost all most differentially expressed genes with proliferation or cell cycle arrest, consistent with previous senescence studies. Investigating expression changes across five consecutive PDs and comparing young with senescent cells enabled us to identify both monotonically up-and downregulated genes as well as the most differentially expressed genes. Both sets of genes strongly contributed to the transition into cellular senescence. Thus, we quantitatively describe similarities in gene expression profiles during the aging of two fibroblast cell lines of different origin.

Data Deposition
The RNA-seq data discussed in this publication have been deposited in NCBIs Gene Expression Omnibus and are accessible through GEO Series accession number GSE63577.