Replicative senescence of HUVECs
We prepared senescent HUVECs, which ceased to proliferate, in the presence of a complete culture medium for comprehensive transcriptome analysis. After 60 days of culture, four independent HUVECs (C1–4) stopped dividing, which showed that the population doubling (PD) rate was <0.1 per day (Fig. 1A). Approximately 70%–80% of the cells showed SA-β-Gal activity, and the cell size was 3.4-fold greater on day 60 than it was on day 22 (Fig. 1B,C). The expression of CDKN2A was induced 10-fold greater on days 73 and 74 than it was on days 5–7 (Fig. 1D). These cellular phenotypes were characteristic of replicative senescent cells, consistent with a previous report (Salama et al., 2014).
Comparative transcriptomic analysis of senescent HUVECs and fibroblasts
To investigate transcriptome responses during endothelial cellular senescence, we conducted RNA-seq analysis of two independent HUVEC clones (C1 and C2), with four different culture periods (5, 18/20, 60/61, and 74 days) (Fig. 2A). First, to compare the gene expression signature with the previous data set of senescence-related genes, we conducted hypergeometric test-based enrichment analysis of senescence-related gene sets extracted from three different databases, namely, Molecular Signatures Database (MSigDB), REACTOME, and Human Ageing Genomic Resources (HAGR) [17-19]. Genes showing differential expression between days 5 and 74 (fold change > 2) were significantly enriched among the known senescence-related gene sets (Fig. 2B). The most significantly enriched gene set of “FRIDMAN_SENESCENCE_UP” was previously reported to be up-regulated in senescent cells (C1, P < 1 × 10−9; C2, P < 1 × 10−4) . The gene set of the p53-mediated pathway and senescence-associated secretory phenotype was enriched among genes both up- and down-regulated in senescent cells (up-regulated genes: C1 and C2, P < 0.005; down-regulated genes: C1, P < 0.05; C2, P < 0.0005) and up-regulated in senescent C1 and C2 clones (C1, P < 1 × 10−4; C2, P < 0.05). Other senescence-related gene sets did not show significant enrichment with genes showing expression changes during HUVEC senescence. These results suggest that the signature expression pattern of some genes is common to senescent HUVECs and fibroblasts, and the remaining genes are involved in a distinct senescence-related pathway in HUVECs.
Next, we examined the distribution of differentially expressed genes (DEGs) and the alternatively spliced genes (ASGs) between young (day 5) and senescent (day 74) HUVECs. The distributions of fragments per kilobase of exon model per million reads mapped (FPKM) in young and senescent HUVECs are shown in Fig. 2C. ASGs were quantified by MISO analysis, with variation in percent spliced-in (ΔPSI) (Bayes score > 10) as shown in Fig. 2D. Because AS occurs in various biological processes, especially cellular differentiation , we analyzed DEGs and ASGs during neuronal differentiation from induced pluripotent stem cells (iPSCs) for 10 days . The results showed a three-fold change in FPKM values, thus providing 1222, 1520, and 3227 genes corresponding to C1, C2, and iPSCs, respectively (Fig. 2E). The magnitude of DEGs during senescence was statistically greater than that during neuronal differentiation from iPSCs (mean values: C1, 3.91; C2, 3.80; iPSCs, 3.64 as a scale of log2), as determined by Wilcoxon tests (C1, P = 0.04; C2, P = 0.0004). Additionally, 637, 451, and 1561 genes were identified in C1, C2, and iPSC data sets, respectively, using the thresholds of ΔPSI > 0.1 and Bayes factor > 10, and the magnitude of ΔPSI in senescent clones was smaller than that in the iPSC data set (mean values: C1, 0.25; C2, 0.27; iPSCs, 0.32), as determined by the Wilcoxon test (C1, P < 1 × 10−22; C2, P < 1 × 10−9) (Fig. 2F). These results suggest that the distribution of DEGs was greater than that of the neuronal differentiation of iPSCs and the pattern of ASGs in senescent clones was smaller than that in iPSCs. This motivated us to conduct a more detailed computational analysis of gene expression and AS.
Key regulators indicating monotonic changes in FPKM during HUVEC senescence
To examine the biological significance of responsive transcripts in senescent HUVECs, we conducted Weighted Correlation Network Analysis (WGCNA) of the FPKM transcript data. We identified four dominant FPKM profile modules (turquoise, blue, brown, and yellow) containing 60.6% of the clustered peak events among 20 different modules (Figs 3A, S1, and S2). These four dominant modules showed a monotonic increase or decrease in responses with the culture period, whereas the other modules showed nonmonotonic patterns in a small number of genes, implying that nonmonotonic patterns detected secondary effects during senescence (Fig. 3A). Next, we conducted a functional enrichment analysis of each set of four dominant modules to identify the pathways affected by senescence. Only modules with decreasing FPKM (turquoise and blue) showed significant enrichment of biological processes such as cell cycle, DNA repair, metabolism, transcription, translation, splicing, and chromosome organization, as determined using the hypergeometric test (P < 1 × 10−4) (Fig. 3B, Table S1). Among these, metabolism and chromosome organization processes were down-regulated in both the turquoise and blue modules, implying that key molecules for cellular senescence are included in these processes.
Next, we analyzed the regulatory network of each of the four dominant modules by determining which genes were involved in the regulation of genes in each module. Genes with high connectivity (>0.8) were defined as hub genes; these genes were considered as key genes in a co-expression network . A hub gene has the potential to affect many genes in its network. We identified 963, 247, 163, and 179 hub genes in the turquoise, blue, brown, and yellow modules, respectively (Fig. S3). To predict the transcriptional regulatory network of hub genes, we conducted a computational analysis using two different databases: the ingenuity pathway analysis (IPA) upstream database and transcription factor binding site (TFBS) database [24, 25]. The IPA upstream regulator analysis revealed potential upstream regulators of genes in each module, and upstream gene analysis using TFBSs in MSigDB revealed upstream gene candidates for each module (Fig. 3C). Canonical senescence-related molecules, identified as upstream regulators such as E2F1, E2F4, RB1, CDKN1A, CDKN2A, and MYC, were included in the turquoise module. Particularly, MYC included in the blue module was identified as an up-regulator in both turquoise and blue modules, suggesting that MYC is a key regulator of not only HUVEC senescence but also fibroblast senescence. Interestingly, the sirtuin (SIRT) family, SIRT2 and SIRT6, was identified as an up-regulator of genes in the blue module. SIRT2 was included in the yellow module as a gene showing a monotonic increase in expression. SIRT1 is a well-known senescence-related SIRT family protein; however, the role of SIRT2 in senescence remains unclear. Taken together, WGCNA and upstream pathway analysis suggest that canonical and non-canonical senescent pathways, such as MYC and SIRT2 signaling pathways, regulate HUVEC senescence.
Predominant changes in AFE- and ALE-type splicing during HUVEC senescence
To examine how HUVEC senescence affects AS, we quantitated the AS events between senescent HUVECs at early (day 5) and late (days 60–74) timepoints during culture (Fig. 4A). We classified the splicing events as alternative 3ʹ splice sites (A3SSs), alternative 5ʹ splice sites (A5SSs), AFEs, ALEs, mutually exclusive exons (MXEs), retained introns (RIs), and skipped exons (SEs) (Fig. 4A,B). The proportion of AFE- and ALE-type splicing increased after HUVEC senescence, and the fold increase in the number of AFEs and ALEs was greater than that in the number of SEs in both C1 and C2 clones during senescence (AFE: 3.1- and 1.5-fold; ALE: 1.8- and 3.1-fold; SE: 1.6- and 1.2-fold, respectively). Among the alternatively spliced events between day 5 and day 74 in HUVEC culture, 74 splicing events overlapped between C1 and C2 clones in a statistically significant manner, as determined by Fisher’s exact test (P = 8.96 × 10−72), suggesting that the identified ASGs depend on HUVEC senescence, not reflecting the difference between C1 and C2 (Fig. 4C). Some of the overlapped splicing events were detected in three senescence-related genes (EFEMP1, FLT1, and TCF3) (Fig. 4D, Table S2). A3SS-, ALE-, and SE-type splicing was altered during senescence in EFEMP1, FLT1, and TCF3, respectively. Among the ASGs, except senescence-related genes, the ΔPSI distribution for the AFE-type splicing of ACACA (acetyl-CoA carboxylase-α [ACCα]) was significantly altered in both C1 and C2 clones during senescence (Bayes factor > 1012) (Fig. 4E). Multiple promoters regulate ACCα mRNA, registered as NM_198836 and NM_198834, in the presence of thyroid hormone, glucagon, and medium-chain fatty acids  (Fig. 4F). Quantitative PCR analysis revealed that the expression of NM_198836 was decreased compared with that of NM_198834 in a time-dependent manner during senescence, implying that AFE-type splicing of ACCα mRNA is involved in HUVEC senescence (Fig. 4G). Taken together, splicing alterations in senescence- and non-senescence-related genes were induced during HUVEC senescence, and AFE- and ALE-type splicing was predominantly altered compared with other types of splicing.
Microexon alterations during HUVEC senescence
Besides canonical splicing analysis, we examined whether microexons are induced by HUVEC senescence using VAST-TOOLS, which detects splicing alterations at all hypothetical splice junctions formed by the usage of annotated and unannotated splice sites . The results of VAST-TOOLS analysis were similar to those of MISO analysis; however, more RI events were detected because of the increased sensitivity for unannotated splice sites (Fig. 5A). We identified several microexon alterations during HUVEC senescence, although to a lesser extent (Table S3). Among these, two microexons in PRUNE2 and PSAP, each containing 9 nt, were altered in both C1 and C2 clones during senescence. To confirm these results, we amplified the microexon-containing region from each isoform by reverse transcription PCR (RT-PCR). Besides C1 and C2 clones, we examined clones C3 and C4. Consistent with VAST-TOOLS results, microexons in PRUNE2 decreased and microexons in PSAP increased the PSI value, with statistical significance or tendency, during HUVEC senescence in all clones (Fig. 5B). Each altered microexon contained the CRAL-TRIO lipid-binding domain and saposin B-type domain in PRUNE2 and PSAP, respectively, implying that microexon alteration affects the function of the encoded protein.
Furthermore, to examine whether readthrough transcripts were induced by HUVEC senescence, we conducted an RNA-seq analysis using deFuse . The readthrough transcript refers to a conjoined gene arising from the upstream transcript to the downstream transcript of partner genes with poly(A) sites [29, 30]. However, RNA-seq data showed that the number of conjoined genes did not change, indicating that HUVEC senescence does not involve conjoined genes (Fig. 5C). Together, our results suggest that a limited number of microexon splicing alterations are induced during HUVEC senescence.