Senescence‐associated transcriptional derepression in subtelomeres is determined in a chromosome‐end‐specific manner

Abstract Aging is a continuous process leading to physiological deterioration with age. One of the factors contributing to aging is telomere shortening, causing alterations in the protein protective complex named shelterin and replicative senescence. Here, we address the question of the link between this telomere shortening and the transcriptional changes occurring in senescent cells. We found that in replicative senescent cells, the genes whose expression escaped repression are enriched in subtelomeres. The shelterin protein TRF2 and the nuclear lamina factor Lamin B1, both downregulated in senescent cells, are involved in the regulation of some but not all of these subtelomeric genes, suggesting complex mechanisms of transcriptional regulation. Indeed, the subtelomeres containing these derepressed genes are enriched in factors of polycomb repression (EZH2 and H3K27me3), insulation (CTCF and MAZ), and cohesion (RAD21 and SMC3) while being associated with the open A‐type chromatin compartment. These findings unveil that the subtelomere transcriptome associated with senescence is determined in a chromosome‐end‐specific manner according to the type of higher‐order chromatin structure.

However, the integration and time-dependent relationship between these aging pathways in different organs and major physiological functions as well as their roles in age-related diseases remain largely unrecognized (Roy et al., 2020). Various aging clocks have been recently described: transcription (Meyer & Schumacher, 2021), methylation (Horvath, 2013), and inflammation (Sayed, 2021), providing interesting tools to assay for biological age but we are still facing a "cause or consequence" challenge. Thus, there is a need to decipher the time-dependent and causality relationship existing between the different aging pathways.
The ends of chromosomes are protected from unwanted DNA damage response (DDR) by the shelterin protein complex that specifically binds to telomeric DNA. In humans, the shelterin is composed of six proteins: TRF1 and TRF2, binding the double-stranded DNA part of telomeres; POT1, binding the 3′-overhang; RAP1, binding to TRF2; and TIN2/TPP1, forming a protein bridge between TRF1/TRF2 and POT1 (de Lange, 2018;Ghilain et al., 2021).
In addition, telomeric chromatin is associated with the non-coding RNA TERRA transcribed from subtelomeric regions and can adopt peculiar conformations involving telomeric DNA looping (T-loop; Griffith et al., 1999). The formation of T-loops, which is facilitated by the shelterin subunit TRF2, contributes to the prevention of DDR checkpoint and repair activation (Benarroch-Popivker et al., 2016;Doksani et al., 2013;Sarek et al., 2019). The replication of telomeric DNA requires specific mechanisms to compensate for the inherent inability of the replication machinery to fully duplicate the extremities of the parental DNA molecule (Gilson & Géli, 2007). In many organisms, this involves a specialized reverse transcriptase, the telomerase (Blackburn et al., 2006).
Telomere shortening is programmed during normal development in humans and other vertebrate species. Indeed, in these organisms, the telomerase expression is downregulated at the end of embryogenesis in somatic tissues leading to telomere shortening at each cell division, eventually resulting in replicative senescence when a subset of telomeres becomes too short to ensure their anticheckpoint functions (Hayflick limit;Abdallah et al., 2009;Bodnar et al., 1998;Kaul et al., 2011). Throughout life, the rate of telomere shortening is paced as a function of the regeneration properties of the considered cell (Demanelis et al., 2020) and in response to external stressors (Jacome Burbano & Gilson, 2021). Telomere shortening can be a driver of human aging since germinal mutations of the telomerase and shelterin complexes, leading to critical telomere shortening and deprotection, cause rare progeroid syndromes such as Dyskeratosis Congenita (Armanios & Blackburn, 2012).
Moreover, short telomeres are associated in the general population with a broad spectrum of age-related pathologies, including cardiovascular and neurodegenerative diseases, type II diabetes, pulmonary fibrosis, and arthrosis (Blackburn et al., 2015;Martínez & Blasco, 2017).
A function of telomeres is to modulate gene expression. This was first discovered in yeast for the process leading to the repression of subtelomeric genes, named telomere position effect (TPE; Gottschling et al., 1990). In brief, in this organism, the major shelterin subunit, Rap1, serves to recruit a heterochromatin complex at telomeres, including the sirtuin Sir2, followed by its spreading toward the centromere (Kueng et al., 2013). The subtelomeric genes succumbing to this TPE silencing are determined in a chromosomeend-specific manner by specific combinations of insulators and proto-silencer sequences, leading to subtelomeric chromatin loops that regulate the heterochromatin spreading in a discontinuous manner (Fourel et al., 1999;Lebrun et al., 2003;Miele et al., 2009). Interestingly, the yeast subtelomeric genes whose expression is regulated by telomeres are related to stress responses and could thus contribute to an adaptive response in case of telomere dysfunction (Ai et al., 2002;Platt et al., 2013). Since these pioneer studies, the existence of TPE was reported in various organisms, including humans (Baur et al., 2001;Koering et al., 2002). Our understanding of the molecular mechanisms of human TPE remains fragmented with the involvement of TERRA (Arnoult et al., 2012), a shelterin subunit (TRF2) (Deng et al., 2009;Kim et al., 2016;Robin et al., 2014), a sirtuin enzyme (SirT6; Tennen et al., 2011), and long-range chromatin loops within the subtelomeric region (Robin et al., 2020;Wood et al., 2014). In particular, the heterochromatin nature of human telomeres remains largely indetermined since in many cell types the telomeres do not exhibit classical heterochromatin marks (Cubiles et al., 2018;Gauchier et al., 2019). Thus, despite striking similarities with yeast, no comprehensive model is established yet for human TPE. Telomeres can also modulate the expression of genes located in the interior of chromosomes by sequestering transcription factors at telomeres Maillet et al., 1996;Marcand et al., 1996). In the case of telomere shortening, their release throughout the nucleoplasm leads to genome-wide transcriptional regulation of targeted genes (Platt et al., 2013;Ye et al., 2014).
The link between telomere and gene expression led to the hypothesis that telomere shortening contributes to the transcriptional changes associated with aging Wright & Shay, 1992;Ye et al., 2014). This view was first supported in yeast replicative senescence (Platt et al., 2013) and recently in human aging transcriptomic studies revealing that telomere length correlates with age-related gene expression (Demanelis et al., 2020) and that genes upregulated with aging are enriched in subtelomeric regions (Dong et al., 2021).
In this work, we analyzed, by RNA sequencing, the transcriptional changes occurring at replicative senescence. We found a significant enrichment of derepressed genes at senescence in regions confined to two megabases from telomeres. Some but not all of them are sensitive to TRF2 and Lamin B1 dosage. Noteworthy, these genes are clustered in a subset of subtelomeres enriched in factors controlling chromatin higher-order structure such as CTCF, MAZ, and the cohesion complex. They are also enriched in chromatin marks of repression by the polycomb complex (EZH2 and H3K27me3) in contrast to the other subtelomeres rather enriched in marks of repression by HP1 (H3K9me3).

| Subtelomeric gene derepression during replicative senescence is restricted to a subset of subtelomeres
To investigate the transcriptional regulation occurring at replicative senescence (RS), we cultured human lung primary fibroblasts (MRC-5 cells) at 5% oxygen until they senesce (population doubling, PD, 71).
We defined a senescent culture when cells stop proliferating and are positive for senescence-associated β-galactosidase (SAβ-Gal) staining, and EdU is incorporated in <1% of the cells ( Figure S1A,B).
We harvested young (PD 30) and RS (PD 71) MRC-5 cells in triplicates, followed by RNA sequencing (RNA-seq) and differential expression analysis. This allowed us to identify 1051 differentially expressed genes (DEGs). When plotted according to the distance to telomeres, the upregulated DEGs found at the first 2 Mb from the telomeres were positively enriched (fold enrichment = 3.6 and 3.8 with p = 7.86 e −15 and 7.13 e −14 , respectively, hypergeometric test, Figure 1a). On the contrary, no specific enrichment was observed in subtelomeric regions for the downregulated DEGs ( Figure 1a).
We hypothesized that the DEG upregulation could result from a senescence-associated process of either derepression or specific transactivation. To decipher between these two possibilities, we compared the expression levels of the DEGs and no-DEGs in young and senescent cells by analyzing the RNA-seq counts. We used the DESeq2 normalized counts of genes located in the first 2 Mb from telomeres, and we separated them into upregulated and downregulated genes ( Figure 1b). As expected, this analysis revealed a significant increase in read counts in senescent versus young MRC-5 of the upregulated DEG group (p < 0.0001, Kruskal-Wallis test, Figure 1b).
Notably, there were significantly fewer reads in the young upregulated DEG group as compared to the no-DEG group (p < 0.0001, Kruskal-Wallis test), indicating that the upregulated DEGs are repressed in young cells compared to the no-DEGs where no change is evident. Therefore, since the upregulated DEGs are less expressed in young cells than the no-DEGs, their upregulation in senescent cells is better explained by a senescence-associated mechanism of derepression.
Regarding downregulated DEGs (Figure 1b), the read counts in the senescent DEG group were significantly lower compared to the young DEG group (p = 0.0151, Kruskal-Wallis test). This indicates that downregulated DEGs in senescent cells, at least in subtelomeric regions, result from senescence-associated transcriptional repression.
Next, we investigated whether the increase of subtelomeric transcription at senescence is a global feature of subtelomeres. We looked for positive fold enrichment of upregulated genes chromosome by chromosome by defining bins of 2 Mb. The enrichment was calculated considering gene density (represented by the grey dots, Figure 1c; see the method section for further details). We found seven chromosome ends showing significant positive enrichment of upregulated genes at the first 2 Mb from the telomeres (chromosome arms 8q, 9q, 11p, 16p, 16q, 19p, and 22q, Figure 1c).
On the contrary, only one chromosome arm (5p) was enriched for downregulated genes ( Figure S2).

| The shelterin protein TRF2 contributes to the senescence-associated subtelomeric gene expression
We asked whether gene expression in subtelomeres can be regulated by the shelterin protein TRF2, which is downregulated as cells approach senescence (Fujita et al., 2010;Mendez-Bermudez et al., 2022) and involved in subtelomeric gene expression (Robin et al., 2014(Robin et al., , 2020. To this end, we infected pre-senescent MRC-5 cells with a TRF2 expressing lentivirus or an empty control and extracted RNA when the cells reached senescence. In the case of the TRF2-expressing cells, senescence was reached after 75 PD in culture, while the control cells transduced with an empty vector stopped dividing at PD 71. In parallel, we ectopically expressed TRF2 in young MRC-5 cells (PD 30) for 6 days before RNA extraction. In both, young and RS cells, the levels of TRF2 were high but constant ( Figure S3A). After RNA sequencing and differential expression analysis, we identified 1366 DEGs when we compared senescent vs young MRC-5 overexpressing TRF2 and 1051 DEGs in senescent versus young MRC-5 transduced with an empty vector ( Figure 2a). By crossing the two datasets, we obtained three groups of DEGs ( Figure 2a) which are dependent upon senescence but that behaved differently according to the ectopic expression of TRF2. The intersection between the two datasets identified 856 genes that are not influenced by the TRF2 levels. We called them TRF2-independent DEGs. We found 510 DEGs detected only when TRF2 is overexpressed; thus, we named this group as TRF2 high -dependent DEGs. Finally, we found 195 DEGs detected only in the empty control group. We assumed that the expression regulation of this group of genes is influenced by the natural decrease of the TRF2 protein level seen in senescent cells and we named this group TRF2 low -dependent DEGs. We confirmed the results obtained from our RNA-seq for some of these DEGs by RT-qPCR ( Figure S3B). Similar to what we found for the whole transcriptome of senescent cells (Figure 1a), enrichment at subtelomeres was found for the upregulated TRF2-independent DEGs (fold enrichment = 3.0 and 4.3 with p = 7.57 e −8 and 1.83 e −13 , hypergeometric test) and TRF2 lowdependent DEGs (fold enrichment = 5.3 and 2.3 with p = 2.87 e −9 and 0.031, hypergeometric test) but not for the TRF2 high -dependent DEGs ( Figure 2b). These results suggest that TRF2 downregulation upon telomere shortening is responsible for the derepression of some of the subtelomeric genes. Interestingly, out of the 195 TRF2 low -dependent DEGs, 27 are subtelomeric and from those, 26 were upregulated DEGs, showing that the naturally low levels of TRF2 at senescence lead to an upregulation of subtelomeric genes.
Next, we asked whether TRF2 levels also regulate the expression of the TRF2 low -dependent DEGs in a non-senescent context, i.e., in young cells of the same cell line. Indeed, downregulating TRF2 in Gene ontology analyses using Reactome reveal that the TRF2independent DEGs belong to pathways related to cell cycle arrest and senescence ( Figure S3D). The other two groups, TRF2 low and TRF2 high -dependent DEGs, did not show a particular pathway enrichment probably due to the low number of DEGs present in those groups.

| Decreased Lamin B1 expression leads to subtelomeric gene derepression
Lamin B1 is an important regulator of chromatin architecture and a component of the nuclear lamina whose expression decreased at senescence (Freund et al., 2012;Shimi et al., 2011). Together with the evidence of interactions between TRF2 and Lamin B1 (Pennarun et al., 2021), the fact that senescence-associated TRF2 We concluded that the Lamin B1 downregulation in young fibroblasts recapitulates in part the transcriptome of senescent cells, including the subtelomere-specific enrichment of derepressed genes.
Noteworthy, the majority of TRF2 low -dependent DEGs are also regulated by Lamin B1 in young cells, suggesting commonalities between TRF2 and Lamin B1 to regulate their expression.

| Specific subtelomeric gene derepression is not restricted to replicative senescence
We then asked whether the subtelomeric gene expression profile identified here in the context of "pure" replicative senescence (i.e., cells grown at 5% oxygen to avoid a combined effect with exogenous oxidative stress) is also observed when senescence is induced by We asked which chromosome ends were enriched for upregulated DEGs in these datasets. We found enrichment in three chromosome ends (9q, 11p, and 16p) for the RS signature and five (3q, 8q, 9q, 16q, and 20q) for the OIS signature (Table S1). Interestingly, most of those chromosome ends overlap with the seven chromosome ends described for our RNA-seq with MRC-5, although not necessarily the same genes are differentially expressed.
It was previously shown that age-dependent DEGs in human tissues were enriched in subtelomeres (Dong et al., 2021). Thus, we asked whether the chromosome ends described in this work (MRC-5 RS cells) were also enriched in age-dependent DEGs. We F I G U R E 1 Distribution of DEGs at senescence. (a) Distribution of DEGs in senescent (PD 71) versus young (PD 30) MRC-5 cells. Upregulated (left) and downregulated (right) genes are shown. The line plots illustrate the number of DEGs related to their distance to telomeres in 1 Mb intervals. Fold enrichment (right Y-axis) is only shown if there is a significant p-value and the size of the circles corresponds to the p-value of the enrichment obtained by a hypergeometric test. (b) Dot plot of DESeq2 normalized counts from the first 2 Mb of telomeres. Counts from non-differentially expressed genes (no DEGs) and from DEGs in the two conditions (young and senescent) are plotted. Upregulated DEG counts are plotted in blue and downregulated DEG counts are plotted in red. Data represent median with interquartile range of three biological replicates. Statistical analyses were performed using the Kruskal-Wallis test, (*p < 0.05; **p < 0.001; ***p < 0.0001; ****p < 0.00001). (c) Upregulated genes in the first 2 Mb from telomeres. The grey empty dots represent all protein-coding genes and pseudogenes. The black dots represent upregulated genes (senescent versus young MRC-5 cells) from this work. The red arrows mark the start of the telomere. The grey boxes delimitate the subtelomeres that are significantly enriched in upregulated genes.
found that five out of the seven subtelomeres were enriched for age-dependent genes (Table S1)  and SMC3, which were also highly correlated with these chromosome ends, are associated with insulation and chromatin looping (Table S2 and Figure 4a,b).
Besides, we found a correlation of several histone marks such as H3K27me3, H4K20me1, H3K36me3, H3K4me3, H3K79me2, and H3K9ac. Amongst those, the polycomb H3K27me3 factor appeared to be an important regulatory element of upregulated DEG subtelomeres as also the histone methyltransferase enzyme Enhancer of Zeste Homolog 2 (EZH2), which catalyzes the addition of methyl groups to histone H3 at lysine 27, was present in our correlation analysis (Table S2 and (Figure 5a and Figure S6) and highly correlated with H3K27me3 histone mark (chr1p, chr16p, chr19p and chr20q).
In agreement with ENCODE data, we found that those DEGs were We further showed that the derepression of some of these genes is dependent upon the decreased levels of TRF2 and Lamin B1 occurring at senescence. Moreover, knocking down TERF2 and LMNB1 genes in young cells is sufficient to derepress their expression.
Since TRF2 is known to mediate the formation of chromatin loops between telomeres and subtelomeres by the interaction with interstitial telomeric sequences (ITS; Kim et al., 2016;Robin et al., 2014Robin et al., , 2020, it is likely that the higher-order organization of subtelomeres controlled by TRF2 plays a role in the senescence-associated gene derepression observed at subtelomeres. One cannot rule out that other telomere factors, potentially able to interact with ITS, such as TRF1 (Simonet et al., 2011) or TERRA (Feretzaki et al., 2020), could mediate chromatin loops affecting subtelomeric gene expression.
These results raised the question of whether this subset of subtelomeres has peculiar chromatin composition. Indeed, as a result of a large screening of ChIP-seq data from primary young fibroblasts, we found that these subtelomeres are specifically enriched with loops (Rao et al., 2017;Xiao et al., 2021). It is also known that cohesin and polycomb proteins interact to maintain chromatin organization (Dorsett & Kassis, 2014;Schaaf et al., 2013). Moreover, the expression of many of those factors decreases upon senescence (Ito et al., 2018).
Clustering analyses showed that the subtelomeres with a high proportion of upregulated DEGs are preferentially associated with the A-type chromatin compartment and become more open upon OIS induction. One limitation of this analysis is the heterogeneity of cell lines and in most ChIP-seq cases the lack of corresponding RNA-seq data. Nevertheless, our analyses unveil that a subset of subtelomeres has specific chromatin features associated with gene upregulation at senescence which could render them more susceptible to gene expression at senescence.
All this raises the question of a role of telomere DNA length in the regulation of subtelomeric gene transcription during replicative senescence. Previously, the expression of the subtelomeric gene, ISG15, has been directly associated with telomere length (Lou et al., 2009). However, we did not find ISG15 in our RNA-seq results nor in the published transcriptome dataset of replicative senescence we used in our analyses. Arguing against a direct role of telomere length is the fact that we also found an enrichment of upregulated subtelomeric genes in OIS, a situation where the average telomere length is unaltered (Suram et al., 2012). However, in the case of OIS, one cannot rule out a telomere effect on subtelomeric gene expression caused by DNA replication stress leading to stochastic telomere loss, which are telomere events that cannot be detected by measuring the mean telomere length (Suram et al., 2012).
Other non-exclusive possibilities are an alteration at senescence of telomere-subtelomeric chromatin loops mediated by factors such as TRF2 (Kim et al., 2016;Robin et al., 2014Robin et al., , 2020, the expression of which is known to decrease upon p53 activation (Fujita et al., 2010;Mendez-Bermudez et al., 2022) or a subtelomeric enrichment of Polycomb-repressed genes, as shown in this work, being upregulated as a consequence of the general attenuation of polycomb repression occurring during senescence (Bracken et al., 2007;Jacobs et al., 1999).
Altogether, these results unveil the existence of a subset of subtelomeres with specific chromatin profiles enriched in genes prone to be derepressed during senescence. This highlights the importance of subtelomeric regions in the senescence process.

| Cell lines
MRC-5 human primary lung fibroblasts were obtained from ATCC.

| EdU proliferation assay
To measure the proliferation of MRC-5 cells, the click-iT EdU Alexa Fluor 647 Imaging Kit (Thermo Fisher Scientific) was used. EdU was incubated for 14 h at a final concentration of 10 μM. Cells were imaged using fluorescence microscopy.

| Chromatin immunoprecipitation (ChIP)
Around 5 million cells were cross-linked for 10 min with 1% formaldehyde and washed with cold PBS. The cells were centrifuged, and the pellet re-suspended in cell lysis buffer ( Cells were visualized using phase-contrast microscopy and the percentage of SAβ -Gal positive cells was calculated.

| Statistical analysis
Statistical analysis was performed using the GraphPad Prism software v9. Student's t-test and the corresponding nonparametric Mann-Whitney test were used for qPCR analysis. The Kruskal-Wallis test was used for DESeq2 counts analysis. Differences were considered statistically significant when the p-value was <0.05.
The p-value of the enrichment analysis was based on the cumulative distribution function of the hypergeometric distribution using the function phyper from the stats R package. Fold enrichment (FE) was calculated using the following equation: where k represents the number of successes, s is the sample size, M is the number of successes in the population, and N is the size of the population.
Regarding the enrichment analysis in the distribution of the DEGs in subtelomeres, the distance to the closest telomere was calculated for each gene. Chromosomes were then divided into 2 Mb non-overlapping intervals according to the distance from a telomere, and DEGs were grouped in those intervals. The enrichment analysis was calculated as described above, with s as the total number of genes of interest, M as the total number of genes in the interval, and N as the total number of genes in the genome.
The correlation analysis between subtelomere ends with upregulated DEGs vs the rest of subtelomeres and different ENCODE factors was performed by calculating the length of the overlap between the coordinates of DEGs and the ones of the ChIP-seq peaks. This overlap was expressed as a percentage followed by a Spearman's rank correlation coefficient computed using the R package lares.
A cluster heatmap (Figure 3b) was generated by calculating the percentage of upregulated genes present in the different chromosome ends. For each data set, the percentage was calculated as the number of upregulated genes in one chromosome end divided by the total number of upregulated genes in all chromosome ends.
Heatmap of C was generated using the coordinates of different ChIP-seq peaks from a selection of young primary fibroblast available in ENCODE, ATAC-seq and Hi-C data. The coordinates were used to calculate the length they covered across each subtelomere and then standardized using the scale function in R. Gilson.

ACK N OWLED G M ENTS
The IRCAN's Molecular and Cellular Core Imaging (PICMI) is sup-

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
RNA-seq data were deposited in the Gene Expression Omnibus database under accession numbers GSE180406 and GSE160503.