Multiple time-series expression trajectories imply dynamic functional changes during cellular senescence.

Cellular senescence is a dynamic process driven by epigenetic and genetic changes. Although some transcriptomic signatures of senescent cells have been discovered, how these senescence-related signals change over time remains largely unclear. Here, we profiled the transcriptome dynamics of human dermal fibroblast (HDF) cells in successive stages of growth from proliferation to senescence. Based on time-series expression profile analysis, we discovered four trajectories (C1, C2, C3, C4) that are dynamically expressed as senescence progresses. While some genes were continuously up-regulated (C4) or down-regulated (C2) with aging, other genes did not change linearly with cell proliferation, but remained stable until entering the senescent state (C1, C3). Further analysis revealed that the four modes were enriched in different biological pathways, including regulation of cellular senescence. These findings provide a new perspective on understanding the dynamic regulatory mechanism of cellular senescence.

Introduction

Cellular senescence is a state of cell growth arrest induced by various factors, and can impact complex organismal processes and disease development, such as cancer [1], atherosclerosis [2], and Alzheimer disease [3], [4]. Cellular senescence is also closely related to organism senescence [5], [6], [7]. Various factors and signaling pathways are involved in intracellular regulation during the proliferative to replicative senescence (RS) stages, including the p16INK4a/Rb and p53/p21Cip1 pathways [8]. In addition to canonical aging markers like p16 and p21 [9], [10], many senescence-associated markers have been identified from single or combined transcriptome datasets based on paired analysis of young and senescent cells [11].

Gene expression is a temporal process. For instance, transcription factors (TFs) can function through distinct trajectories to initiate or repress transcription under different conditions [12], [13]. Previous studies have shown that cellular senescence is a highly dynamic multistep process, during which cells exhibit a range of phenotypically diverse cellular states [14]. Mounting evidence emphasizes the significance of time-series gene expression data in understanding developmental processes and disease progression [15], [16], [17]. However, systematic analysis of the transcriptomic profiles of cells at continuous time points has not yet been reported. To fill this knowledge gap, we used gene expression profiling of replicative senescent cells at distinct time points to reveal intrinsic expression trajectories during cellular RS. An in-depth understanding of the molecular mechanisms underpinning progressive aging may provide new perspectives and strategies for studies on cell senescence.

Materials and methods

Cell strains and culture

The human dermal fibroblast (HDF) cell line was isolated from an adult and cultured in Dulbecco’s Modified Eagle Medium (DMEM)/High glucose medium (C11995500BT, Gibco) with 10 % fetal bovine serum (35–076-CV, Gibco) and 1 % penicillin/streptomycin (15140–122, Gibco) in a 37 ℃/5% CO2 incubator. For RS, cells were propagated in culture from passages 10 to 34 and collected for transcriptome analysis every three generations from passage 13.

Senescence-associated beta-galactosidase (SA-β-Gal) staining

SA-β-Gal staining was conducted using a Senescence Cells Histochemical Staining kit (CS0030, Sigma-Aldrich) according to the manufacturer’s instructions. Cells were visualized using a Nikon eclipse Ti inverted microscope.

Cell proliferation assay

The HDF cells were seeded into 96-well plates at adensity of 8,000. The DNA synthesis was measured using the Cell Light EdU Apollo567 In Vitro Kit (C10338-1, RiboBio) following the manufacturer’s protocols. In brief, cells were labeled with 10 µM EdU for 24 h, then fixed in 4 % paraformaldehyde for 30 min, then stained with Apollo 567 and Hoechst 33342. The cells were imaged by Cytation 5 Cell Imaging Multi-Mode Reader (BioTek).

RNA isolation and sequencing

Total RNA samples were extracted using the TRIzol method (Invitrogen). The polyA-enriched RNA-sequencing (RNA-seq) libraries were prepared for sequencing using the Illumina HiSeq 4000 platform.

Public datasets

Raw RNA-seq data of the MRC5 (human lung fibroblasts) and HFF (human foreskin fibroblasts) cell lines [18] were collected from the ‘‘GEO Repository’’ (https://www.ncbi.nlm.nih.gov/geo/) (Accession number GSE63577). The fastq files were downloaded using the SRA Toolkit v2.10.8. The RNA-seq data of human foreskin fibroblasts (HAC2) subjected to ionizing radiation were obtained from the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress) (ID code: E-MTAB-5403).

Bioinformatics analysis

FastQC v0.11.9 (https://github.com/s-andrews/FastQC) and fastp v0.20.1 [19] were used for quality control of all raw data, including public datasets and our own. The clean data were then aligned to the human reference genome (hg19) using Hisat2 v2.1.0 [20]. The FeatureCounts [21] command in Subread v1.6.0 was used to count reads of genes. The raw count matrix was normalized as Transcripts Per Million (TPM). Only genes annotated as protein-coding were included in analysis. To identify significant differential expression changes with cellular senescence, we used the lm function in R 4.1.2 to fit the linear regression model. Significance was determined a priori at a Benjamini-Hochberg (BH)-corrected p-value of < 0.05.

Time-series analysis

The gene expression time-series data were clustered using the fuzzy c-means algorithm in the Mfuzz package [22] of the R platform, and genes with consistent expression changes were grouped into a cluster.

Functional enrichment analysis

Differentially expressed genes (DEGs) were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using the clusterProfiler v4.2.1 package [23] and visualized using the “ggplot2” package in R. Functional enrichment analysis of multiple gene lists in the four patterns was conducted using Metascape [24] (https://metascape.org). Protein-protein interaction network (PPI) analysis was carried out using the PPI analysis function embedded in Metascape, with the hub gene network identified using the Molecular Complex Detection (MCODE) algorithm. Network visualization was generated through Cytoscape (v3.9.0) [25].

Results

Transcriptome time-series expression profiles of HDF cells during RS

To better understand the molecular changes during RS, we conducted RNA-seq on HDF cell samples spanning eight time points, from proliferation (passage 13) to senescence (passage 34), with two biological replicates per passage (Fig. 1). We first verified the senescent level of cells by senescence-associated beta-galactosidase (SA-β-Gal) staining cell percentages and EdU-staining cell ratio (Fig. S1A and S1B). Next, we examined the expression patterns of classical senescence markers (p16/CDKN2A and p21/CDKN1A) to confirm the senescence state of the HDF cells. As expected, the mRNA expression levels of p16 and p21 increased significantly with cellular senescence in the HDF cells (p16: p = 3.55e-07; p21: p = 0.027) (Fig. 2A, 2B), and p16 protein expression was also high in the cells at passage 34 (p = 0.0066) (Fig. 2C). The protein level of p16 and p21 was validated by western blotting assay (Fig. S1C). In addition, we also analyzed the expression changes of the senescence-associated secretory phenotype (SASP) factors (i.e., CCL2, CXCL8, IL1B, IGFBP5, IGFBP7, IL-6, ANGPTL4 and VEGFA) in old and young HDF cells. Results showed that the expression of all SASP genes except one was significantly up-regulated in senescent HDF cells (P34) (Fig. S1D). As the loss of heterochromatin is an important feature of cellular senescence, we analyzed and found the heterochromatin regulator gene CBX5 (HP1-α) was significantly down-regulated along with senescence in both mRNA and protein level (Fig. S1E).

Fig. 1

Pipeline for transcriptome sequencing analysis.

Fig. 2

Transcriptomic changes in HDF cells with senescence. (A, B, C) Expression profiles of classical senescence markers at the mRNA and protein level. The lm test in R was used to examine differences between passages. (D) Volcano plot of DEGs in senescent cells. Red dots represent up-regulated genes, and blue dots represent down-regulated genes. (E) Functional enrichment analysis of DEGs. Nine KEGG pathway terms (adjusted p < 0.05) are shown. (F, G) Presentation of time-series analysis of DEGs. (F) Four gene clusters were obtained from 3113 DEGs based on Euclidean distance and c-means objective function [22]. The four trajectories were described as “stable first and then decreasing” (“C1”), “continuously decreasing” (“C2”), “stable first and then increasing” (“C3”), and “continuously increasing” (“C4”). (G) Statistics of number of genes in four clusters. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To detect changes in gene expression during proliferation and senescence, we identified significant DEGs using linear regression (BH-corrected p < 0.05). In total, 3 113 DEGs were obtained, including 1 300 up-regulated and 1 813 down-regulated genes (Fig. 2D and Table S1). Enrichment analyses indicated that these DEGs were enriched in adherens junction, endocytosis, and insulin resistance processes (hsa04520, hsa04144, and hsa04931) (BH-adjusted p < 0.05) (Fig. 2E and Table S2), which play important roles in cellular senescence [26], [27], [28], [29]. To explore the patterns of genetic changes during RS, we performed time-series analysis of the expression profiles of DEGs across eight passages using the Mfuzz R package. Interestingly, the DEGs were grouped into four expression trajectories, with the up-regulated and down-regulated genes clustered into two patterns, respectively (Fig. 2F, G). Although previous research has reported on transcriptomic signatures at different time points in irradiation-induced senescence (IRIS) [11], we depicted the four gene expression trajectories from proliferation to senescence.

Dynamic changes in different functions during RS

To elucidate the underlying functions of the four patterns, we conducted functional enrichment analysis of genes in the four clusters. Comparing functional enrichment among the four expression trajectories, we found that genes in the different clusters were significantly and specifically enriched in distinct functional pathways (p < 0.01) (Fig. 3 and Table S3). Genes in the “C1” pattern cluster were specifically enriched in spliceosome, lysine degradation, base excision repair, Notch signaling pathway, and VEGF (vascular endothelial growth factor) signaling pathway (hsa03040, hsa00310, hsa03410, hsa04330, and hsa04370), indicating that these functions were stable in the cell proliferation phase but decreased with growth arrest. In contrast, the C2-specific functions, including ribosome, sphingolipid signaling pathway, and NOD-like receptor signaling pathway (hsa03010, hsa04071, and hsa04621), decreased continuously with RS. Among the up-regulated genes, the HIF-1 signaling process (C4-specific function) (hsa04066) increased continuously with cellular aging, whereas glycan biosynthesis (hsa00510, hsa00514) and Wnt signaling process (hsa04310) (C3-specific function) only increased during the senescent stage. Of note, almost all pathways mentioned above play functional roles in the regulation of cellular senescence [30], [31], [32], [33]. In particular, the ribosome biogenesis [34], [35], base excision repair [36], [37], and VEGF signaling pathways [38], [39] are essential for senescence progression. Taken together, these results suggest that different pathways function in distinct stages during cellular senescence.

Fig. 3

Functional enrichment analysis of genes in four clusters. Enrichment analysis of multi-gene list was conducted using Metascape. Pathway terms with p < 0.01, minimum gene count of 3, and enrichment factor > 1.5 were retained.

Gene expression trajectories in different cell lines and senescent types

To determine whether the four expression patterns were confined to the HDF cell line or also existed in other cells, we collected public datasets of other human fibroblasts (HFF and MRC5 cells) subjected to RS [18]. Interestingly, results showed that the four expression patterns were validated in at least one other cell type, indicating that these trajectories were stable in distinct cell lines during RS (Fig. 4A and Fig. S2). For each trajectory, we compared genes in the three cell lines, and obtained 437 genes (C1: 126 genes; C2: 207 genes; C3: 61 genes; C4: 43 genes) (Fig. 4A and Table S4).

Fig. 4

Time-series analysis of different cell lines and types of senescence. (A) Conservation of four gene clusters in different cell lines. Left panel: Four gene expression patterns in three cell lines (HDF, HFF, MRC5). Right panel: Heatmap of 437 overlapping genes in three cell lines in four gene clusters. Figure shows fold-change of genes in three cell types relative to proliferating cells. (B) Four expression patterns in HAC2 cells subjected to ionizing radiation.

We also explored whether the four expression trajectories were stable in cells subjected to different senescence-inducing stimuli. First, we conducted a time-series analysis of the RNA-seq dataset at different time points in HAC2 cells after ionizing radiation [11]. Surprisingly, the four expression trajectories also existed in IRIS (Fig. 4B and Fig. S3). Next, we confirmed whether the dynamic functions were consistent in the different types of cellular senescence. Unexpectedly, enrichment analysis revealed that most functional pathways of genes in the four trajectories differed from those in RS (Fig. S4 and Table S5). Nevertheless, the ribosome pathway (hsa03010; C2), valine, leucine, and isoleucine degradation process (hsa00280; C2), and mitophagy function (hsa04137; C3) consistently changed in both RS and IRIS. Among genes of the mitophay pathway, PINK1 is the only gene shared in RS and IRIS. PINK1 has proved to be a master regulator of mitophagy, and its silencing could suppress mitophagy [40]. Prior studied have emphasized the importance of homeostasis of mitophagy in the regulation of cellular senescence, and either too high or too low is detrimental to cell growth [41]. The lysosome-related function (hsa04142) was stable at the early stage and then increased with RS (C3) but was continuously enhanced in IRIS (C4). Of the 437 genes identified in the four trajectories in RS, 62 showed consistent expression changes in IRIS (Fig. 5).

Fig. 5

Key genes in four gene clusters independent of cell type and senescence-inducing type. Genes meeting MCODE criteria are marked in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We next identified key functional genes in the four expression trajectories in the HDF cells using PPI enrichment analysis and MCODE algorithm. Using default parameters, the physical interactions in STRING [42] (physical score > 0.132) and BioGrid [43] were retained, resulting in 116 genes involved in 228 interaction edges after MCODE analysis (Fig. S5), with BCAT2 and TAB1 conserved in the HAC2 cells during IRIS.

In conclusion, the four expression trajectories were stable in the different cell lines during RS progression, but only a few functional pathways were retained in cells during IRIS.

Discussion

Most biological processes change dynamically over time, and analysis of time-series experiments can provide important information regarding how certain functional pathways respond to various signals. Here, we found that gene clusters with distinct functions exhibited changes in expression in distinct trajectories. Further enrichment analysis showed that certain cellular senescence-related processes did not change linearly with growth, but remained stable during the proliferative phase, with linear down-regulation or up-regulation not occurring until the senescence phase.

We also detected DEGs with cellular senescence, 33 of the DEGs identified here are involved in the 55 core senescence-associated genes identified by Hernandez-Segura et al. [11] (Fig. 2 and Table S1). As cellular senescence is a dynamic process [14], [44], we discovered four DEG expression trajectories in the HDF cells across sequential passages. Intriguingly, we found that some pathways did not change continuously with growth time (C1, C3). Despite the heterogeneity of senescence in different cell types, the four expression trajectories were also stable in HFF and MRC5 cells. Coincidentally, previous research has shown that certain pathways exhibit time-point specific responses to ionizing radiation [11]. In our study, the four expression trajectories were also found during IRIS, but the genes were different from those involved in RS. Differences between replicative and radiation-induced senescence have been investigated previously [45], [46]. Consistently, we found that pathways related with cell cycle and base excision repair (hsa04110 and hsa03410) were more stringently regulated in RS other than IRIS. Notably, we found that ribosome and amino acid catabolism functions continued to weaken with cellular senescence in both RS and IRIS progression. Indeed, previous studies have suggested that diminished ribosome biogenesis could drive cell into senescence [30].

Taken together, we identified four gene expression trajectories during cell proliferation and senescence. Clarifying these four trajectories should shed new light on critical biological processes ranging from proliferation to senescence. Our study also highlights the importance of time points in cell growth when exploring the regulatory mechanisms that underpin cellular senescence. Cellular senescence is a complex process with heterogeneity in the regulatory mechanisms of different cell types. We only included fibroblast-based transcriptome data in our study, and thus further research on the dynamics of senescence in other cell types is needed.

CRediT authorship contribution statement

Ming-Xia Ge: Data curation, Formal analysis, Visualization, Writing – original draft, Writing – review & editing. Qin Yu: Investigation, Resources, Validation. Gong-Hua Li: Data curation. Li-Qin Yang: Resources. Yonghan He: Supervision, Funding acquisition. Ji Li: Resources. Qing-Peng Kong: Conceptualization, Visualization, Supervision, Project administration, Funding acquisition, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.