Analysis and pharmacological modulation of senescence in human epithelial stem cells

Abstract Human epithelial stem cells (ESCs) are characterized by long‐term regenerative properties, much dependent on the tissue of origin and varying during their lifespan. We analysed such variables in cultures of ESCs isolated from the skin, conjunctiva, limbus and oral mucosa of healthy donors and patients affected by ectrodactyly‐ectodermal dysplasia‐clefting syndrome, a rare genetic disorder caused by mutations in the p63 gene. We cultured cells until exhaustion in the presence or in the absence of DAPT (γ‐secretase inhibitor; N‐[N‐(3, 5‐difluorophenacetyl)‐L‐alanyl]‐S‐phenylglycine T‐butyl ester). All cells were able to differentiate in vitro but exhibited variable self‐renewal potential. In particular, cells carrying p63 mutations stopped prematurely, compared with controls. Importantly, administration of DAPT significantly extended the replicative properties of all stem cells under examination. RNA sequencing analysis revealed that distinct sets of genes were up‐ or down‐regulated during their lifetime, thus allowing to identify druggable gene networks and off‐the‐shelf compounds potentially dealing with epithelial stem cell senescence. These data will expand our knowledge on the genetic bases of senescence and potentially pave the way to the pharmacological modulation of ageing in epithelial stem cells.


| INTRODUC TI ON
The homeostasis of the epithelial tissues is ensured by the presence of adult stem cells (SCs), which are endowed with two important features: self-renewal and unlimited potency. 1 Recent evidence has demonstrated that in vivo SCs reside in specific regions of the tissue and remain quiescent until activation is required either due to the normal need to maintain tissues or in response to diseases and injuries. 2 The structure and localization of such niches are tissue-specific and depend on tissue turnover and regenerative potential. The ability of SCs to survive and retain their proliferative potential throughout their lifespan does not necessarily imply that they have an endless capacity to divide, undergoing constant self-renewal. As tissues have different developmental needs and cellular turnover rates, the in vivo self-renewal frequencies of SCs are different. Interestingly, the behaviour of the different epithelia reflects the physiological role of the tissues they belong to. The epidermis undergoes keratinization, a process in which the epidermal cells progressively mature from basal cells with proliferative potential to the lifeless, flattened squames of the stratum corneum to generate a functional barrier to the external environment. The incessant exposition to thermal shocks, pathogens and chemical agents prompts the skin epithelium to undergo a continuous replacement of cells, which involves cell proliferation and differentiation that culminate in the desquamations of keratinized flakes. The homeostasis is guaranteed by epithelial stem cells (ESCs) located in the interfollicular epidermis (IFE), with a turnover time of 28 days. 3 The oral mucosa, that protects the oral cavity from bacteria and from the mechanical stress during mastication, has similar functions. It is possible to distinguish the different areas of the mucosa according to differences in keratinization. The gingivae and hard palate are keratinized, while the floor of the mouth, the buccal regions located under the surface of the tongue, and the sulci are not keratinized. The oral mucosa has higher numbers of cell layers compared with the epidermis ( Figure 1A). [4][5][6] Quiescent cells are located in the basal layer, while proliferation is restricted to the parabasal and adjacent suprabasal layer. [6][7][8][9][10][11][12] The ocular surface, which is composed of the cornea centrally and the sclera peripherally, shares the same ectodermal origin.
Located in the outermost part of the eye, the cornea is vulnerable to damages caused by burns, abrasions, contact lenses, alterations in tear production, infections and other diseases. 13 The corneal epithelium consists of a single basal layer and of stratified squamous epithelial cells and is essential to guarantee the structural uniformity and the transparency of the tissue. During corneal homeostasis, human limbal epithelial stem cells (H-LESCs) proliferate and generate transient amplifying cells (TA) that divide, differentiate and migrate to the centre of the cornea to regenerate the epithelium. 13 H-LESCs are present in the basal layer of the limbus, in specific niches called Palisades of Vogt. 14 Unlike the H-LESCs, the precise location of the human conjunctival stem cells (H-CESCs) is still controversial. Although previous studies suggested that H-CESCs are uniformly distributed throughout the whole ocular surface, 15,16 recent reports have demonstrated that they are located in the fornix and/or in the bulbar conjunctiva. 17,18 In fact, the fornix may provide greater physical protection, intraepithelial mucous crypts, vasculature and immune cells, typical features of the SC niches.
If SCs are exhausted too quickly, or if genetic defects alter their proliferative potential, tissue atrophy and premature ageing can arise.
For example, mutations in the p63 gene, known to cause ectrodactylyectodermal dysplasia-clefting (EEC) syndrome, induce a rapid exhaustion of clonogenic and self-renewal potential of ESCs, resulting in accelerated ageing. [19][20][21] Conversely, mutations that promote frequent SC divisions without an appropriate differentiation balance can result in abnormal tissue development and even cancer. 22 In this study, we demonstrate that the frequency and timing of SC divisions are tightly regulated and specific for each tissue, as well as preserved in vitro to ensure a defined lifelong maintenance of the SCs population. In addition, the treatment of cells with the γ- S-phenylglycine T-butyl ester) resulted in the extension of the lifespan of ESCs. Finally, we identified several genes strongly related to the senescence process, by means of RNA sequencing (RNAseq) and gene ontology (GO) analyses, reconstructed the protein networks involved and identified potential ageing modulating drugs.

| Cell culture
Primary human keratinocytes were isolated from fresh skin, conjunctiva, oral mucosa and limbus biopsies (N = 12) of healthy donors, aged Clonal analysis was performed to obtain holoclones, meroclones and paraclones, as previously described. 23 DAPT was administered as previously described. 24

| Evaluation of the relative telomere length
Cytogenetic analysis of cells was performed according to standard procedures. 25 Fluorescence in situ hybridization (FISH) was done on chromosome preparations according to manufacturer's recommendations and as previously reported. 26 Commercially available all-telomere probes were applied (Telomere PNA Kit, Agilent, Santa Clara, CA, USA).
Pictures were taken with a Zeiss Axioplan Imaging microscope (Jena, Germany and ISIS software (MetaSystems, Altlussheim, Germany).
Telomere FISH signals were quantified on a per-cell basis using the open-source plugin "Telometer" program (version 3.0.6, available at: https://demar zolab.patho logy.jhmi.edu/telom eter/). FISH images were converted to 16-bit grayscale TIFF images using the ImageJ program (available at: https://imagej.nih.gov/ij/). The Differences between means were determined through a Multiple Range Test using Bonferroni's analysis. Statistical significance was set at p < .05.

| RNA isolation, library preparation and sequencing
Total RNA was isolated from cultured keratinocytes at each passage and subsequently purified with the RNeasy Clean-Up Kit  Overall, the preprocessing scheme selected 3807 isoforms and 4824 genes, which will undergo further analysis. To assess the quality of the preprocessed data set, we computed the coefficient of FPKM variation (CV) for each gene and isoform.
Low-quality bases were trimmed using erne-filter 28 ; residual adapter sequences were removed using cutadapt. 29 Cleaned, trimmed reads were aligned against the UCSC hg19 reference human genome using tophat2. 30 Transcript quantification based on the genome annotation was performed using cufflinks. 31 Differential expression was performed using cuffdiff. 32 Temporal evolution of expression was assessed by measuring the ratio of isoform expression in FPKM to the expression of GAPDH in FPKM.
Expression at the first time-point is set to 0 and the expression at following time is represented as the log2 variation of the expression relative to time 0. Genes presenting a pattern of constant up-or down-regulation over time were selected as interesting candidates and plotted. Temporal evolution of ΔNp63α was also plotted.

| Gene ontology analysis
Gene ontology (GO) enrichment analysis was performed using the Enrichment analysis tool available from the GO consortium, [33][34][35] to identify biological processes significantly enriched in genes either up-regulated or down-regulated. To this end, the PANTHER Overrepresentation Test (Released 20,200,728) was performed using the GO Ontology database DOI: 10.5281/zenodo.4033054 Released 2020-09-10. Genes and isoforms either up-regulated or down-regulated were analysed using the Homo sapiens reference list and the GO biological process complete Annotation Data Set, using the Fisher's exact test with FDR correction.

| Protein association analysis and druggable targets identification
Associations between up-and down-regulated genes were analysed using the String suite 36 and the following settings: network type, full network; meaning of network edges, confidence; active interaction sources, all; medium confidence. Druggable gene products were identified using the Genecards suite 37 and by literature mining.

| Quantitative real-time PCR
Total RNAs from primary cell cultures were extracted and purified using the RNeasy Micro kit (Qiagen, Milan, Italy). The purity of the RNA preparation was verified by measuring its absorbance ratio at 260/280 nm. 500 ng of RNA was used to synthesize the cDNA with random hexanucleotide primers and MoMULV reverse transcriptase (Fisher Scientific, Milan,) at 42°C for 1 h. 50 ng of cDNA was amplified in an AB7900 real-time PCR detection system (Fisher Scientific, Milan, Italy) using TaqMan™ Universal Master Mix II (Fisher Scientific, Milan, Italy), in a total volume of 20 μl. For the absolute ΔNα-tp63 quantification (Ab-qPCR), the level of expression of the target gene was normalized to GAPDH (Glyceraldehyde-3-Phosphate Dehydrogenase).
For the relative gene expression analysis (Rel-qPCR), the difference in relative target gene expression was performed using the 2 − ΔΔCt method. GAPDH was used as internal control gene. The efficiency of the target amplification (ΔNp63α) and of the reference amplification (GAPDH) were measured and found approximately equal.

| Statistical analysis
Differences between groups were analysed by Student's t-test, while multiple comparisons were performed by one-way or two-way anova as indicated in the text. All statistical analyses were performed using Graphpad Prism 9.0. For data shown in Figure 6, analysis of covariance was performed. The test fits independent linear models for each time series; the coefficients of the models are then compared using standard methods to identify the genes that have statistically different expression dynamics.

| In vitro characterization of human primary epithelial keratinocytes reflects native and specific differences of the respective tissues of origin
We compared the proliferative potential of ESCs from different sources in vitro. To this end, human primary epithelial keratinocytes derived from skin, limbus/cornea, conjunctiva and oral mucosa were isolated ( Figure 1A), cultured in vitro and serially propagated until exhaustion, with CFE being quantified at each passage. We also ana-

| Karyotype analysis, telomere length and mitochondrial activity throughout DAPT induced proliferation
To assess any possible transformation of cells following DAPT administration, we performed a karyotype analysis on EEC-OMESCs at different passages in culture (pIV, pVII, pXVII). The cytogenetic assay confirmed that no numerical or structural chromosomal abnormalities were observed compared with cells cultivated without DAPT (i.e., cells at passage pI, the initial and, therefore, standard of reference for lifespan analysis) ( Figure 5A-D).
Q-FISH analysis of the telomeres of EEC-OMESCs showed increased telomere length in cells grown with DAPT (p IV, pV), compared with the same cells cultivated without DAPT ( Figure 5E).
These results correlate well with the extended lifespan observed in SCs treated with continuous administration of DAPT, thus confirming that elongation of telomeres is linked to increased cellular half-life. 38 Since the senescence process in vitro is likely to involve both telomeres and mitochondria, 39 we evaluated the mitochondrial activity in EEC-OMESCS after DAPT administration by means of a MitoTracker assay, and analysed the signal of probes, which passively diffuse across the plasma membrane and accumulate in active mitochondria. As shown in Figure S2 Figure S2F.
Because a low mitochondrial activity has been previously reported to be linked to a higher SC potential in ex vivo experiments, 40 our results seem to confirm that the administration of DAPT enriches the SC population, as shown by the low mitochondrial activity of treated cell, and delayed premature senescence.  (Table S1), whereas 30 genes and 8 isoforms were consistently expressed at higher levels in the first passages, thus suggesting a potential down-regulation with senescence (Table S2). We further performed a GO enrichment analysis on the up-regulated and down-regulated gene lists reported in Tables S1 and S2, to identify biological processes involved in ESC senescence. We could identify 27 biological processes associated with down-regulated genes ( Figure 7A) and 6 associated with up-regulated genes ( Figure 7B).

| Comparative epithelial stem cell transcriptomics analysis
Our results indicate that SC senescence is associated to a downregulation of cellular products mainly associated with "negative regulation of gene silencing", "prophase", "mitotic prophase", "histone H3-K27 trimethylation" and "regulation of chromatin condensation" (all with an enrichment >95x). On the other hand, up-regulated genes were mainly linked to inflammation and immune response.

| Analysis of marker expression during lifespan
We extrapolated the expression profiles of a panel of genes involved in the maintenance of telomeric integrity (RAD50, F I G U R E 3 DAPT administration extends the lifespan of stem cells. Cell cultures derived from holoclones and meroclones obtained from donor limbus and EEC-OMESCs were treated with DAPT. The number of passages in the absence (red columns) or presence (blue columns) of DAPT (A). Significant increase in the number of passages after DAPT treatment of EEC-OMESCs (B). Clonogenic cells, final cell number and ∆Np63α expression quantified by qPCR from holoclones (C-E), meroclones (F-H) and R279 H-OMESCs (I-K) either untreated or treated with DAPT. Data are shown as mean ± standard deviation along with the results of statistical significance as calculated by ordinary one-way anova with multiple comparisons (*p < .05; **p < .005; ***p < .0005; ****p < .0001). EEC-OMESCs: ectrodactyly-ectodermal dysplasia-clefting oral mucosa epithelial stem cells; DAPT: N-[N- (3, 5- seems to be strongly down-regulated in LH cultures, as expected, while the initial expression levels are maintained during the different passages in cells treated with DAPT, thus confirming that the senescence process is delayed in this culture. 48 Overall, these data confirm that following treatment with DAPT, ageing process of LH cultures seems to slow down.

F I G U R E 7
Gene ontology (GO) enrichment analysis of deregulated genes in epithelial SCs. Genes either down-regulated (A) or upregulated (B) during SC lifespans were used to perform GO enrichment analysis as described in the materials and methods section. All significantly enriched biological processes are shown on the left: the size of each slice is proportional to the enrichment of each biological process with respect to the reference list, and the p value is indicated. All other relevant information is reported in the right panels F I G U R E 8 Expression profiles of senescence markers. The expression profiles of a panel of genes involved in the maintenance of telomeric integrity and/or associated with differentiation and senescence in keratinocyte stem cells were obtained from Limbal Holoclone cultures either untreated (LH) or treated with DAPT (LH-DAPT) by means of RNAseq analysis. Normalized expression of the indicated markers is plotted against the culture passage. The time series relative to LH and LH-DAPT is compared using an analysis of covariance test (p < .05). LH, limbal holoclones

| Identification of druggable protein networks involved in epithelial stem cell senescence
Because our data showed that it is possible to modulate the senescence of ESCs pharmacologically by continuous administration of DAPT, we sought to identify additional drugs capable of further modulating such process. To this end, we reconstructed a protein networks based on the RNAseq data and identified potential drugs able to modulate the senescence process based on the information available on the Genecard suite ( Figure 9). Our results suggest that five down-regulated and four up-regulated gene products are targetable by 19 different approved or experimental drugs. Among down-regulated gene products, CDK1 and TOP2A, associated with 10 and 7 other factors, respectively, appeared the most central druggable factor ( Figure 9A). Among druggable up-regulated factors, STAT1 was associated with three additional gene products ( Figure 9B).

| DISCUSS ION
The functional decline of SCs has been associated with many age- The ability of DAPT to delay cell senescence, and therefore extend the replicative properties, was also tested on oral mucosa ESCs derived from three EEC patients (R279H-OMESCs). It is known that cells from these patients are affected by premature senescence.
Interestingly, the addition of DAPT increased the number of cell passages in vitro from 7 ± 2 to 18 ± 2 i.e., about three times more than the normal length of the lifespan of these cell cultures. The analysis of the karyotypes of EEC-OMESCs treated with DAPT did not reveal any numerical or structural chromosomal abnormalities. However, both the telomere length and the number of passages in culture were increased, thus confirming that the elongation of telomeres is associated with an increased half-life of the cells (27). Our preliminary findings would therefore indicate that the use of DAPT could slow down the senescence of epithelial SCs from patients with EEC syndrome, thus suggesting interesting pharmacological opportunities for the treatment of these patients.
Our data show a differential longevity of the analysed epithelia in vitro, greater for the skin and lower for the cornea (Figure 2), which is consistent with the extension and the physiological function of the tissue. 51 Accordingly, the limbus has a smaller number of SCs than its neighbouring conjunctiva, which requires more SCs for normal renewal as having a larger extension. While the niches of the cornea are located in a specific and delimited area (the limbus) in which the palisades of Vogt are organized to protect the SCs, the distribution of the niches in the skin encompass the entire basal layer. This is likely due to the epithelium of the skin being constantly subjected to stress and continuously renewed. 52 In our previous studies, [19][20][21] we provided evidence that a major factor in the pathophysiology of EEC syndrome is the premature de- • CLU, a molecular chaperone that increases stress resistance and significantly extends lifespan, when overexpressed 58 ; • SEPW, a selenoprotein involved in redox-related processes, muscle growth and differentiation, and in the protection of neurons from oxidative stress during neuronal development 59 ; • STAT1, a signal transducer and transcription factor playing an essential role in response to interferon (IFN) signalling and regulating many cellular processes, such as proliferation, differentiation and cell death. 60,61 Among the down-regulated genes, we found the following ones to be interesting: Finally, numerous gene involved in regulatory processes important for development, survival, proliferation and differentiation were found up-(KIAA1191, S100A9, SEPW1 and TMBIM1) or downregulated (CSRP2, IFI16 and PSIP1).

| CON CLUS IONS
In the current study, we provided some insights into the self-renewal and senescence features of epithelial tissues. We show that the in vitro proliferative capacities of epithelial SCs from skin, cornea, conjunctiva and oral mucosa tissues, are different and lineage-specific and reflect the different physiological needs and characteristics of each tissue. We provide evidence that lifespan may represent a valid model for studying ageing and senescence features in vitro and that it can be extended with the help of suitable pharmacological inhibitors.
Moreover, we report a list of genes, the expression of which strongly changes during SC lifespan. Additionally, we show that DAPT treatment can slow the senescence process by extending the replicative capacity of SCs, further confirming the importance of Notch signalling in ageing. 66 Finally, we identified a number of approved or experimen-

ACK N OWLED G EM ENTS
We thank Dr. Patrizia Nespeca, Dr. Greta Bernardo and Dr Giulia Masi for their support in the expression analysis of SC markers, Claudia Breda for her support with cell cultures and Dr. Angelo Migliorati for the contribution in the graphic elaboration of Figure 1.

CO N FLI C T O F I NTE R E S T
The authors confirm that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.