Age, gender and UV-exposition related effects on gene expression in in vivo aged short term cultivated human dermal fibroblasts.

Ageing, the progressive functional decline of virtually all tissues, affects numerous living organisms. Main phenotypic alterations of human skin during the ageing process include reduced skin thickness and elasticity which are related to extracellular matrix proteins. Dermal fibroblasts, the main source of extracellular fibrillar proteins, exhibit complex alterations during in vivo ageing and any of these are likely to be accompanied or caused by changes in gene expression. We investigated gene expression of short term cultivated in vivo aged human dermal fibroblasts using RNA-seq. Therefore, fibroblast samples derived from unaffected skin were obtained from 30 human donors. The donors were grouped by gender and age (Young: 19 to 25 years, Middle: 36 to 45 years, Old: 60 to 66 years). Two samples were taken from each donor, one from a sun-exposed and one from a sun-unexposed site. In our data, no consistently changed gene expression associated with donor age can be asserted. Instead, highly correlated expression of a small number of genes associated with transforming growth factor beta signalling was observed. Also, known gene expression alterations of in vivo aged dermal fibroblasts seem to be non-detectable in cultured fibroblasts.

Introduction

Various biological effects have been related to ageing including accumulation of DNA damage and reactive oxygen species (ROS), metabolic alterations (especially energy metabolism) and cellular senescence [1]. Many of these ageing related effects, like accumulation of mutations in somatic and mitochondrial DNA and telomere attrition are not directly linked to altered mRNA expression. But still, transcriptomic profiling, by outlining many physiologic effects in parallel, provides sensible information on global uniformly age associated effects.

For whole transcriptome analysis, RNA-seq is an established platform [2]. For each step of analysis (alignment to genome, differential expression analysis, functional classification), a variety of standard technologies are available. Gene expression in normal cells [3] and age-related alterations of gene expression [4] have been found to be very tissue specific. Age-related gene expression changes are also species specific [4].

Ageing skin

During the ageing process, human skin undergoes characteristic morphological and functional changes, for example, reduced epidermal thickness, flattening of dermo epidermal junction [5] and the reduction of fibrillar collagen content [6-9]. In young skin, thick collagen fiber bundles are present with little open space. Fibroblasts appear orientated along collagen bundles. In old skin, collagen fibers are more disorientated with present empty space and fibroblasts show little orientation along fibroblast bundles [7]. Additionally, aged skin contains increased amount of fragmented collagen [9]. Also, the in vitro growth capacity of aged fibroblasts is reduced [5, 7].

Homeostatis of dermal extracellular matrix

Fibroblasts produce the dermal collagen matrix consisting of 80–90% of Type I collagen and 10–15% Type III collagen. For both types of collagen, a linear age-related decrease of 29% over a 49-year period in cultured fibroblasts has been reported [10]. Gene expression of Type I procollagen has shown to be reduced by 75% in fibroblasts from direct dermis extracts [7, 8].

The content of matrix metalloproteinase 1 (MMP1) is elevated in aged upper dermis and in aged fibroblasts [11]. Additionally, an increased content of fragmented collagen, the product of collagen degradation by MMP1 can be found in aged skin [12].

Regulation of collagen production by TGF-β

Type I procollagen production is mainly regulated by the transforming growth factor beta (TGF-β) / SMAD signalling pathway [8]. Connective Tissue Growth Factor (CTFG) is an important downstream regulator of TGF-β/ R-SMAD mediated reduction of collagen and levels of mRNA TGF-β1 and CTGF decrease by 70% and 57% respectively in in vivo aged dermal fibroblasts [8]. Regulation of Type I collagen expression by TGF-β is mediated via Type II TGF-β receptor (TβRII) and SMAD3 in dermal fibroblasts [13] and expression of SMAD3 is reduced in aged skin [14]. A universal negative feedback regulation mechanism of TGF-β signalling is mediated via SMAD7 [15] which has also been shown to be active in dermal fibroblasts [8]. Age-related dermal reduction of collagen content therefore is a consequence of concomitant diminished expression of TGF-β, SMAD and CTGF in dermal fibroblasts.

In vivo interaction of dermal fibroblasts with extracellular matrix

Extracellular mechanical forces, mediated by Integrin receptors (e.g. Integrin α2β1) induce prominent alterations of cellular function (e.g. gene expression) in dermal fibroblasts [11, 16]. Presence of fragmented collagen increases Integrin α2β1 mRNA which in turn causes strong up-regulation of MMP-1 mRNA in dermal fibroblasts. Additionally, presence of reactive oxygen species (ROS) is markedly increased by contact with extracellular fragmented collagen. Based on these findings, a model of self-perpetuating age dependent collagen fragmentation has been proposed [11].

Increased ROS also down-regulate TGFBR2 (TGF-β receptor II) thereby impairing the TGF-β/ R-SMAD pathway [13]. In aged skin, dermal fibroblasts have lower cellular contact area to collagen fibers resulting in less mechanical stimulation [7]. Therefore, in vivo ageing effects of dermal fibroblasts depend on at least two different mechanisms:

Cell intrinsic senescence

Altered mechanical interaction with extracellular matrix

Genes associated with known cellular alterations

Known alterations associated with ageing skin as well as cellular senescence mechanisms draw attention to genes evidently related therewith. For skin, known gene groups with altered expression include collagens, metalloproteinases and TGF-β signalling pathway associated genes. With cellular senescence, functionally associated genes are TGF-β signalling pathway and cyclin dependent kinases [17]. Additionally, for senescent cells, commonly utilised reporter genes are p16INK4a(CDKN2A) and p21Cip1 (CDKN1A, WAF1) [18]. The indicator SA-βgal is mainly used for histochemical detection [19, 20]. Also, proteins secreted as part of Senescence-Associated Secretory Phenotype (SASP), for example, interleukins, proteinases or some chemokines are of interest [19, 21, 22].

Photoageing

Photoageing, the most important causative of extrinsic skin ageing [6], is mainly caused by UVA (320–400 nm) exposition [23, 24]. Characteristic for UV aged skin are deep wrinkles and reduced stiffness [6].

In UV radiated skin various molecular effects are detected, for example, reduced TGF-β signalling causing diminished collagen production [25].

Principle findings

We describe that consistent age-related alterations in gene expression are not detectable in short term cultured fibroblasts. Using an analysis approach based on monotone alignment depth ratios, we filtered out 42 genes with consistently increasing or deceasing alignment depth. In a subset of 9 TGF-β signalling related genes (ATOH8, SNAI1, ID3, SPHK1, ID1, PRRX2, SMAD7, FAM83G, SERTAD1) high pairwise correlated gene expression was found (correlation coefficients >0.8).

Materials and methods

Sample donors

Sixty 4 mm punch biopsies of healthy skin were taken from 30 human donors. From each donor one sample was derived from a sun-exposed site (neck/shoulder) and one sample from a sun-protected (buttock/gluteal) site. Fifteen donors were female and fifteen donors were male. Samples were analysed histologically to verify presence of characteristic age-related or UV-exposition related alterations. Sample donors were assigned to three groups according to their age: Young (18–25 years), Middle (35–49 years) and Old (60–67 years). Six samples (two samples from each age group) were excluded from further analysis due for serious disturbing effects identified by analysis of DNA 6-mer spectra using the Bioconductor package seqTools [26]. Samples from 27 subjects (13 female and 14 male) were included into the analysis (in total 54 samples, 18 samples per age-group). The study was approved by the Ethical Committee of the Medical Faculty of the University of Düsseldorf (# 3361) in 2011. Informed written consent was obtained from all donors before sample acquisition.

Fibroblast cultures

From each sample a part was fixed in 3.4% formaldehyde and prepared for histology using standard procedures. Isolation and primary culture of the cells were performed as previously described [27]: 4-mm punch biopsies were washed in 70% ethanol, followed by sterile phosphate-buffered saline. In brief, skin pieces were incubated with dispase (10 mg/ml in phosphate-buffered saline, sterile filtered) at 37°C and 5% CO2 for 2 hours to remove the epidermis. Dermal pieces were dried for 10 min under the sterile bench, then culture medium was added. Cells were propagated further at 37°C, 5% CO2 and in high glucose and L-glutamine containing Dulbecco’s modified Eagle’s medium (DMEM; Biochrom) supplemented with 10% fetal calf serum (FCS Superior; Biochrom) and 1% antibiotics (PAN Biotech) with the medium changed every 2–3 days. The fibroblasts were sub-cultured by treatment with 0.05% trypsin for 5 min. Mycoplasma and virus contamination was excluded by the Multiplex Cell Contamination Test (Multiplexion). Fibroblasts started to migrate out of the dermal pieces after approximately 1–2 weeks. Fibroblasts for RNA-seq analysis were taken from passages 2–7.

RNA-sequencing

For alignment and all subsequent analysis, Human genomic sequence (GRCh38) and annotation data (release 82) were downloaded from Ensembl [28] and BioMart [29, 30]. cDNA libraries were synthesised using TruSeq RNA SamplePrep kit (Illumina) according to the manufacturer’s protocol. One microgram of total RNA was used for poly(A) RNA enrichment. The samples were amplificated on 9 Illumina flow cells (v1.5) and sequenced on a Illumina HiSeq 2000 sequencer using TruSeq SBS kits v1. From each lane, the resulting 101-nt sequence reads were converted to Fastq by CASAVA (1.8.2). Subsequent alignments were calculated on unprocessed Fastq files. Alignments were calculated using bowtie2 (2.2.5) [31], TopHat (2.0.14) [32] and STAR (2.4.1d modified) [33]. The raw Fastq files are available under ArrayExpress accession E-MTAB-4652 (ENA study ERP015294).

Analysis of alignment counts

Analysis of alignment counts were executed using R / Bioconductor software [34]. Read counts were calculated from BAM files using two different approaches: summarizeOverlaps (GenomicAlignments) [35] and readExpSet(rbamtools) [36]. For standard differential expression (DE) we used Quasi-likelihood F-Tests from the edgeR (3.12.0) framework [37] (EQLF). P-values were corrected for multiple testing using Benjamini Hochberg procedure. Adjusted p-values (FDR levels) below 0.1 were considered significant. Analysis of gene set enrichment (GSEA) for GO terms and KEGG pathways was done using functions topGO and topKEGG from package limma (Bioconductor) [38]. For GESA, p-values (P.DE), below 0.05 were considered significant.

Estimation of mRNA expression level

Gene expression levels were estimated with normalised read counts derived from summarizeOverlaps (GenomicAlignments). The read counts are presented as counts-per-million (CPM) values. We do not use model-based CPM estimates (e.g. obtained using cpm (edgeR)) because we want demonstrate the variability of CPM values within groups.

Alignment gap-sites

Analysis of read counts in this study utilises our R/Bioconductor analysis pipeline specialised on identification of splicing events consisting of the R packages rbamtools [39], refGenome [40] and spliceSites [41]. The central surrogate of splicing events in read alignment data are gap-sites: Inner borders of gapped read alignments, sharing coordinates with all other reads covering the same splicing event. In this analysis, read numbers on gap-sites are used for identification of splicing events present in alignment data and for estimation of gene expression abundance.

Gene filtering using monotone alignment depth (AD) ratios (MALDR approach)

For identification of age-related differential gene expression, a sensible approach should focus on small differences which monotonically increase with proceeding age surrounded by considerable inter-individual variation. Besides standard differential expression (DE) procedures, we therefore implemented a filter method based on monotone alignment depth relations (MALDR). The initial step in this procedure uses read counts on gap-sites. Gap-sites identified from TopHat alignment data using readExpSet (spliceSites) were first filtered for presence of at least one read in each sample and then filtered for annotated splice sites in Ensembl 82. The method consists of the following steps:

Count aligns on annotated gap-sites using readExpSet (spliceSites).

Perform edgeR differential expression analysis (Quasi likelihood F-test) on gap-site alignment counts and on two age groups: Young and Old.

Filter out genes for which DE analysis resulted in FDR below 0.1 on at least one splice site.

Count alignment depth on complete genomic regions for all filtered genes and all samples.

First cut of genomic regions: Intronic regions (defined by alignment gap sites) are removed from alignment counts resulting in a restricted genetic region.

Combine alignment depth data for each group using mean and smooth mean group alignment depth data using loess regression.

Second cut of genomic regions: Regions with low read coverage (below 2% of maximal alignment depth within gene) are removed from alignment counts.

Determine regions with monotone alignment depth relation (Young

Determine fraction of restricted genetic region where relation between all adjacent groups is >1.2.

Genes, for which fraction of restricted genetic region is >99% of restricted genetic region are considered significantly differential expressed.

We present MALDR results in increasing fraction of restricted genetic region order containing no information on effect size or baseline expression level, a property shared with ordering along p-values. We utilised MALDR approach for a more sensible detection of age-related alterations in gene expression level.

Testing for differences in different qualities

The samples analyzed in this study are grouped by three donor qualities: Age, gender and sample location. As we expected to find small group differences on a small set of differential expressed genes, we conducted three parallel test procedures (one for each quality) where each test ignores a heterogeneity with respect to other qualities. With this strategy, we maintained the full sample size for each quality. The hypothesis of small numbers of differentially expressed genes was confirmed by the results of this study where the number of significant different expressed genes is below 3 ‰ in all comparisons.

Results

Histological analysis

We show two representative histological images of biopsies from two female donors aged 24 and 60 years. Thereby we can show the characteristical age-dependent loss of rete ridges and thinning of the epidermis including a potentially reduced stratum corneum and a less dense dermal collagen structure (Fig 1). In addition, sun-exposed sites from the same donor (shoulder versus gluteal) are characterised by an increased number of melanocytes and in particular in the skin of the older donors an even more profound disturbance of the dermal collagen structure. More histological sections are shown in supplemental material (see S2 File).

Fig 1

Representative histological sections from two female donors.

A: Female donor, 24 years old (Young age group), gluteal (sun-protected) location. B: Female donor, 60 years old (Old age group), gluteal (sun-protected) location.

Analysis of age related differential gene expression

Standard differential expression analysis for age group

We performed differential expression analysis restricted to two age groups (Young and Old). Each group consists of 18 samples. Standard differential expression analysis of gene-wise read counts with EQLF (egdeR, Quasi-likelihood F-Test) resulted in no significant differential expressed gene (lowest FDR-value: 0.23).

Differential expression analysis for age group MALDR approach

We therefore re-evaluated the aligned reads using MALDR approach. From 54 samples, we identified 1,000,380 different unique gap-sites from which 100,040 reside on Ensembl 82 annotated splice sites. From testing number of alignments on each gap-site using edgeR QLF test we obtained 790 genes with at least one significant gap-site. Subsequent filtering using the beforehand described MALDR approach resulted in 42 genes with age dependent variation in gene expression. For these 42 genes, we examined alignment depth data:

CPM values (as estimates for alteration of gene expression level) derived from summarizeOverlaps (shown in S3 File).

Gene-wise maximum of readExpSet derived counts of gapped alignments crossing gap-sites (shown in S4 File).

Results from both methods have high similarity. They show that no consistent age dependent progress is present (in accordance with results from EQLF framework). Rather, the differences identified by MALDR base on different degrees of inter-individual variation in different age groups which is further illustrated in the following paragraph.

Exemplified results for gene ID1

We show the loess smoothed alignment depth data for gene ID1 (Inhibitor of DNA binding 1) after intronic and low coverage regions have been removed. A clear monotone ratio effect can be seen for age-group mean values of alignment depth (Fig 2). Plotting CPM values from individual samples (separated by gender and location) reveals that still for the majority of individuals, no systematic tendency is present (Fig 3). Obviously, the observed gene expression differences identified by MALDR are due to the presence of three out of 27 individuals (a male and a female at age of 25 and a female at age of 36) with markedly elevated transcript levels for ID1 (assigned to age groups Young and Middle). These three individuals generate a larger variation and increased mean values in the Young and Middle aged groups. We later on will explore, whether other genes with the same expression pattern can be identified in this data set and examine correlation coefficients for CPM values therefore.

Fig 2

Align depth estimates for gene ID1.

The figure displays alignment depth in absolute numbers. Three lines estimate mean alignment depth for each age group (y = Young, m = Middle, o = Old).

Fig 3

Gene-wise CPM values for ID1.

Gene-wise counts per million (CPM) values derived from summarizeOverlaps.

MALDR derived result set

In order to differentiate the effect in MALDR filtered genes from a global, by differential expression analysis identified change we call these 42 genes age MAR (Monotone Alignment depth Ratio) genes. The complete table of age MAR genes is available in supplemental material S1 File. The number of age MAR genes is small compared to what has been reported for age related differential expression of genes in other tissues for example, in artery (3082 genes [4]) or blood (3287 genes [4] or 1497 genes [42]). But it is in the same order of magnitude as reported for skin (12 genes [4]) and dermal fibroblasts (104 genes [43]; In contrast, Glass et. al. describe 1672 differentially expressed genes in skin [44]).

Comparison age MAR genes from this study with age related differential expressed genes identified in the recently published GTEx study [4] resulted in no intersecting gene. Changing the FDR limit to 0.1 increased the number of age MAR genes in skin to 55 genes; we finally intersected our list herewith.

Functional characterisation of age MAR genes with GO and KEGG

Standard analysis for significantly enriched GO terms resulted in 609 terms (shown in S1 File). Thereof, only five GO terms were associated with fibroblast physiology (e.g. GO:0017134 fibroblast growth factor binding) or extracellular matrix (e.g. GO:0010715 regulation of extracellular matrix disassembly). Analysis for enriched KEGG pathways resulted in 12 enriched terms (shown in Table 1). Only for TGF-β signalling (path:hsa04350) and Rap1 signalling pathway (path:hsa04015) a relation to fibroblast physiology is plausible. The other associations are likely to be incidental.

Table 1

Enriched KEGG pathways in age MAR genes.

Pathway	N	Age MAR genes
TGF-beta signaling pathway	67	ID1, ID3, SMAD7
Signaling pathways regulating pluripotency of stem cells	96	ID1, ID3
Hippo signaling pathway	112	ID1, SMAD7
Rap1 signaling pathway	137	FGF13, ID1
Metabolic pathways	905	(4 genes)
Pentose phosphate pathway	23	PRPS1
Propanoate metabolism	28	ACSS1
Arginine and proline metabolism	33	CKB
Sphingolipid metabolism	36	SPHK1
VEGF signaling pathway	44	SPHK1
Arrhythmogenic right ventricular cardiomyopathy	47	GJA1
Melanoma	47	FGF13

Pathway: Name of KEGG pathway, N: Total number of genes in pathway, Age MAR genes: Number of selected genes in pathway

Functional characterisation of age MAR genes based on literature

We therefore explored putative physiologic roles of the 42 age MAR genes directly from literature by searching PubMed and online resources (HGNC, http://www.genenames.org/ and GeneCards, http://www.genecards.org/). Outlines of functionality with reference to literature is listed for each age MAR gene in the supplementary material (see S2 File). For every gene, we determined whether association with cellular senescence, cancer progression (or cell cycling), embryonal development (or tissue differentiaton), is described in the literature because of the close functional relationships between senescence, tissue differentiation, proliferation and cancer. The categories and associations of age MAR genes are shown in Table 2.

Table 2

Functional characterisation of age MAR genes.

Function	Associated genes
Cellular senescence(11 genes)	ATOH8, ID3, ID1, ERRFI1, MEG3, STC1, HSPB7, SMAD7, FAM83G, CKB,SERTAD1
Cancer progression / Cell cycle(20 genes)	PODXL, SNAI1, ID3, SPHK1, ID1,ERRFI1, SEPT5, MEG3, STC1, PRRX2,SMAD7, FAM83G, DDR1, EVA1A,FGFRL1, FILIP1L, ENC1, SERTAD1,ADGRL4, KCNC4
Differential expressed on cancer(11 genes	PODXL, SPHK1, ERRFI1, MEG3, PRRX2,DDR1, EVA1A, FILIP1L, ENC1, ROBO1,KCNC4
Development / Differentiation(10 genes)	ATOH8, PODXL, SNAI1, ID3, ID1,TRNP1, PRRX2, FAM83G, ENC1, SERTAD1
Cell migration / Adhesion(5 genes)	PODXL, ERRFI1, DDR1, FILIP1L, ROBO1
Cellular filaments / Skeleton(4 genes)	PODXL, SEPT5, HSPB7, ENC1
Vascularisation / Endothel(5 genes)	PODXL, FILIP1L, ADGRL4, ROBO1,KCNC4
Apoptosis / Autophagy(6 genes)	SPHK1, HSPB7, PRRX2, EVA1A, FILIP1L,GJA1

Functional associations of age MAR genes based on available literature. Genes possibly associate with multiple categories.

Correlation of gene expression between functional related genes

From our list of 42 age MAR genes, 8 genes (ATOH8, SNAI1, ID3, SPHK1, ID1, CNN1, PRRX2, ROBO1) are regulated by TGF-β. SMAD7 is part of TGF-β signalling and two genes (FAM83G, SERTAD1) interact with SMAD proteins. lncRNA MEG3 regulates TGF-β pathway genes. Altogether, TGF-β signal related genes comprise more than 25% of our age MAR genes. Analysis of correlation coefficients of age MAR gene CPM values revealed that for a set of 9 genes (ATOH8, SNAI1, ID3, SPHK1, ID1, PRRX2, SMAD7, FAM83G, SERTAD1), all pairwise correlation coefficients are greater than 0.8 (detailed analysis in supplemental material S2 File).

In five genes (ATOH8, ID3, ID1, SMAD7, FAM83G) gene expression in our samples show a similar expression pattern: CPM values (derived from summarizeOverlaps) for these genes are markedly elevated in the same three individuals (and they are therefore classified as age MAR genes). Fig 4 shows a panel of scatterplots of CPM values for genes ATOH8, ID3, ID1, SMAD7 and FAM83G revealing this relation. Highly correlated gene expression of SMAD7, ID1 and ATOH8 has indeed been described in mouse liver cells where the gene expression correlates with iron content of liver cells and is regulated by BMP6 [45].

Fig 4

CPM values for five genes.

CPM (counts per million) values derived from summarizeOverlaps for Genes ATOH8, ID3, ID1, SMAD7 and FAM83G for all 54 samples.

Analysis of gender and location related differential gene expression

Due to absence of age-related differential expressed genes in EQRF analysis, we regard the population as homogeneous with respect to donor age and age group and analysed differential expression of genes for gender and sample location (UV exposition) using EQLF.

Analysis of gender related DE genes

We analysed differential expression between samples derived from female and male sample donors (sex biased gene expression). Testing with EQLF revealed 168 gender related differential expressed genes (gender DE genes). The complete gene list is contained in S1 File (Gender_EQLF_DE_genes). Therefrom, 95 (57%) genes were higher expressed in male individuals (male-biased) and 73 (43%) genes were higher expressed in female individuals (female-biased). The chromosomal distribution of gender DE genes is shown in Fig 5. Twofold over-representation is restricted to Y chromosome for male-biased genes and is present on chromosomes 2, 10, 11, 15, 17, 19 and X for female-biased genes. Functional annotation with GO and KEGG reveals predominantly metabolic processes over-represented (e.g. GO:0006694 and path:hsa00100, both steroid biosynthesis).

Fig 5

Chromosomal distribution of gender DE genes.

Raw Number of significant gender DE genes per chromosome. On chromosome 18, no gene was differential expressed.

Comparison of gender DE genes with other tissues

In order to compare the list of gender DE genes with results from other tissues, we downloaded lists of differential expressed genes from two mircroarray studies: A study on liver samples revealing 1,018 gender DE genes [46] and a study using peripheral blood samples revealing 649 gender DE genes. In Fig 6, a Venn diagram shows the number of shared genes between three tissues: fibroblasts, liver and peripheral blood. All study overlaps are below 10% of gender DE genes of each sample.

Fig 6

Number of gender related DE genes in different tissues.

Comparison of gender related differential expressed genes. Gender DE genes in liver were described by [46]. Gender DE genes in peripheral blood were described by [47].

Analysis of location (UV exposition) related DE

From each donor, two samples were taken from different locations: One sample from the gluteal region (UV protected) and one sample from the shoulder (UV exposed). Differential expression analysis using EQLF results in 56 location related differential expressed (location DE) genes. The complete gene list is contained in S1 File (Location_EQLF_DE_genes).

Analysis of gene expression with respect to fibroblast physiology

Collagen Type I gene expression

From collagen Type I it is known, that in vivo gene expression of dermal fibroblasts progressively decreases with age. Therefore, the fact that in our samples, collagen Type I is not identified as differential expressed by EQLF and not filtered by MALDR approach probably requires further exploration. Age-group-wise alignment depth relations as shown in Fig 7 indicate, that collagen gene expression consistently is lowest in the Middle aged group explaining, that monotone alignment depth relations are absent for this gene. Alignment count values (CPM) show how this age group relation arises. Fig 8 shows that individuals expressing high amounts of COL1A1 (CPM >25,000) follow different patterns in female and male subjects. For male subjects, all COL1A1 high expressing individuals are >40 years old. For female subjects, one participant is 65 years old and all other are less than 30 years old.

Fig 7

Alignment depth for gene COL1A1.

Align depth in genetic region COL1A1 after cutting out intronic regions and regions with low alignment depth. Group-wise mean alignment depth values have been smoothed using loess regression.

Fig 8

CPM values for gene COL1A1.

Gene-wise counts per million (CPM) values directly derived from summarizeOverlaps.

Genes associated with cellular senescence

We tabulated gene expression level (logCPM) and changes in gene expression (logFC) for genes known to be associated with cellular senescence and SASP (see S1 File). Although, for some genes logFC values considerable different from 1 are present (-0.7 for TBX3 or 2.75 for CCL13) though changes are not significant (all FDR values are 1).

Discussion

Known effects of in vivo aged fibroblasts

In our data, consistent age-related changes gene expression for COL1A1, TGFB1, CTGF (CCN2) and MMP1 are absent. Therefore, the major known age-related physiological alterations in in vivo aged dermal fibroblasts are not present in our data. The apparently arising conflict is resolved by the fact, that changed gene expression for these genes have been shown to be undetectable in cultured dermal fibroblasts. For example, gene expression of procollagen Type I, TGF-β and CTGF does not differ between in vivo young and aged dermal fibroblasts [8]. Expression of MMP-1 is increased in aged dermal fibroblasts in vivo, but no difference has been observed in untreated cultured fibroblasts [11].

Fibroblasts interact very intensively with the extracellular fibrillar matrix (via integrin receptors) [7, 11] and physical forces regulate many important biological processes. Reduced mechanical forces down regulate TGFBR2 thereby impairing the TGF-β/ R-SMAD pathway. As consequence gene expression for Type I collagen, fibronectin and CTGF is reduced [48]. Removal of fibroblasts from tissue environment into cell culture is therefore likely to have profound impact on regulatory networks.

Gene groups in differential expressed genes

We filtered out 42 age-related differential expressed genes using MALDR approach. Therefrom, a subset of 8 genes show highly correlated gene expression all of which are functionally related to TGF-β signalling pathway. The conclusions drawn from this observation are weakened by the fact that identified age-related alteration is mainly due to high expression of these genes in few individuals.

Effects of photoageing

UV radiation induces a variety of changes in gene expression in the skin. Expression of Type I procollagen is reduces in human dermis and this reduction is mediated via TGF-β signalling pathway [25]. mRNA levels of JUN (c-jun) and FOS (c-Fos), members of transcription factor AP-1 family, are increased while the gene expression of SMAD 7 is up-regulated [12, 49]. Expression of matrix metalloproteinases (MMPs) is increased [50, 51]. These effects are not reflected in our location DE genes, implicating that chronic effects of UV radiation may be different from acute effects.

Common age-related effects on gene expression

There are numerous studies on age-related physiological alterations and especially on gene expression. Age differences in gene expression show large variation between tissues as well as between species [4, 52].

Consistently age-related differential expressed genes

A central objective for high-throughput studies on samples from differently aged individuals is whether a global (i.e. not tissue specific) age related gene expression signature can be identified. Although there have been results suggesting that these gene sets may exist [53], already early studies reported that no such gene sets are present [54]. Results by Glass et. al. [44] suggest the same conclusion (by finding only one common age MAR gene in three tissues).

Functional enrichment analysis

Functional enrichment analysis of differential expressed genes commonly identifies mitochondrial pathways, RNA metabolism, and immune response pathways [4, 42, 53].

Enrichment of these pathways can not be concluded from our data. KEGG pathway analysis on our data is based on a very small set of genes (i.e. one gene for most pathways) and therefore statistical results are not reliable.

Experimental validation

The functional background of correlated gene expression in five genes (ATOH8, ID3, ID1, SMAD7, FAM83G) related to cellular senescence could further be explored using gene knockdown experiments, for example in the TGF-β pathway. As recently has been shown that reduction of mechanical forces leads to downregulation of TGF-β pathway [48], further cultivation likely reproduces the effect of cultivation and results could not directly be related to skin ageing.

A second question is whether the results from this transcriptomic study translate into protein levels especially when possibly insufficient correlations are considered. In a recent study, it was shown that predictability of protein levels from transcript level could be significantly enhanced when gene-specific RNA-to-protein (RTP) conversion factors are used leading to median pearson correlation coefficients of 0.93 [55]. As gene-specific RTP factors would not influence the results of gene-wise testing, it can be assumed that the major conclusions of this study also apply to protein levels.

Concluding remarks

Ageing from evolutionary theory perspective

Evolutionary mechanisms select out morbidity and mortality before reproductive (and upbringing) ages but not beyond (selection shadow) [56]. Since ageing is not under evolutionary selection [57], single genes propagating (or preventing) ageing are unlikely to exist [56]. In contrast, since absence of selection evokes variation, a spectrum of ageing phenotypes with considerable inter- and intra-individual discrepancy can be expected. Human lifespan, a related property, is a complex trait where associated genomic loci have small effect size [58-60].

A component of evolutionary ageing models is antagonistic pleiotropy: Properties with beneficial effects in early life (which are under selection) and late detrimental effects (in the selection shadow) [56]. The tight connection between cellular differentiation, senescence, cancer and ageing now provides a model for antagonistic pleiotropy where cellular senescence is essential for embryonic tissue differentiation and stress response [61] but deleterious when organisms age or when cancer cells escape from senescence (a model which similarly has already been proposed [62, 63]).

Conclusion

Limitations of current study

Sample size

There is considerable variation in biological ageing processes, between individuals, tissues (and species). For capturing systematic effects as well as underlying variations, large sample sizes are needed [64]. Early size estimates were in the range of 16–63 samples per group [54, 64]. As tissue ageing appears gradually over time, experimental settings should be able to detect small fold changes (>2 [64]) which further increases sample size requirements [65]. Therefore, this study is under-powered, especially regarding two additional sources of variation (apart from age: Gender and UV exposition).

Gene tables.

Spreadsheet file containing the following tables:

Age_MALDR_DE_genes: Listed age-related differential expressed (age MAR) genes.

Age_MALDR_GO: Listed enriched GO terms for age MAR genes.

Age_MALDR_KEGG: Listed enriched KEGG pathways for age MAR genes.

Gender_EQLF_DE_genes: Listed gender related differential expressed (gender DE) genes.

Gender_EQLF_GO: Listed enriched GO terms for gender DE genes.

Gender_EQLF_KEGG: Listed enriched KEGG pathways for gender DE genes.

Location_EQLF_DE_genes: Listed location related differential expressed (gender DE) genes.

Location_EQLF_GO: Listed enriched GO terms for location DE genes.

Location_EQLF_KEGG: Listed enriched KEGG pathways for location DE genes.

Literature_selected_genes: Listed DE values for genes associated with senescence or fibroblast ageing.

(XLSX)

Click here for additional data file.

Supplemental material.

Background information on methodology and additional results.

(pdf)

Click here for additional data file.

Age MAR CPM values.

PDF file with gene expression data for the 42 age-related DE genes (derived from edgeR QLF test and ReadExpSet).

(pdf)

Click here for additional data file.

PDF file with summarizeOverlaps derived CPM data for the 42 age-related DE genes.

(pdf)

Click here for additional data file.