# Single-cell transcriptomics reveals age-resistant maintenance of cell identities, stem cell compartments and differentiation trajectories in long-lived naked mole-rats skin

Aleksandra Savina, Thierry Jaffredo, Frederic Saldmann, Chris Faulkes, Philippe Moguelet, Christine Leroy, Delphine Del Marmol, Patrice Codogno, Lucy Foucher, Antoine Zalc, et al.

► To cite this version:

Aleksandra Savina, Thierry Jaffredo, Frederic Saldmann, Chris Faulkes, Philippe Moguelet, et al. Single-cell transcriptomics reveals age-resistant maintenance of cell identities, stem cell compartments and differentiation trajectories in long-lived naked mole-rats skin. Aging (Albany N.Y.), 2023, 14 (9), pp.3728-3756. 10.18632/aging.204054. hal-03977899

HAL Id: hal-03977899

https://cnrs.hal.science/hal-03977899

Submitted on 7 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Single-cell transcriptomics reveals age-resistant maintenance of cell identities, stem cell compartments and differentiation trajectories in long-lived naked mole-rats skin

Aleksandra Savina<fnref n="1" />, Thierry Jaffredo<fnref n="2" />, Frederic Saldmann<fnref n="3" />, Chris G. Faulkes<fnref n="4" />, Philippe Moguelet<fnref n="5" />, Christine Leroy<fnref n="6" />, Delphine Del Marmol<fnref n="7" />, Patrice Codogno<fnref n="6" />, Lucy Foucher<fnref n="8" />, Antoine Zalc<fnref n="1" />, Mélanie Viltard<fnref n="3" />, Gérard Friedlander<fnref n="6" />, Selim Aractingi<fnref n="1" />,<fnref n="9" />,*, Romain H. Fontaine<fnref n="1" />,*

**Correspondence to:** Romain H. Fontaine; email: romain fontaine@inserm.fr

**Keywords:** naked mole-rat, skin, stem cells, wound healing, aging

**Received:** December 1, 2021

**Accepted:** April 25, 2022

**Published:** May 4, 2022

**Copyright:** © 2022 Savina et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/3.0/) (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

## ABSTRACT

Naked mole-rats (NMR) are subterranean rodents characterized by an unusual longevity coupled with an unexplained resistance to aging. In the present study, we performed extensive *in situ* analysis and single-cell RNA-sequencing comparing young and older animals. At variance with other species, NMR exhibited a striking stability of skin compartments and cell types, which remained stable over time without aging-associated changes. Remarkably, the number of stem cells was constant throughout aging. We found three classical cellular states defining a unique keratinocyte differentiation trajectory that were not altered after pseudo-temporal reconstruction. Epidermal gene expression did not change with aging either. Langerhans cell clusters were conserved, and only a higher basal stem cell expression of *Igfbp3* was found in aged animals. In accordance, NMR skin healing closure was similar in young and older animals. Altogether, these results indicate that NMR skin is characterized by peculiar genetic and cellular features, different from those previously demonstrated for mice and humans. The remarkable stability of the aging NMR skin transcriptome likely reflects unaltered homeostasis and resilience.

## INTRODUCTION

Skin acts as an essential barrier and protects organism from external threats, preventing fluid loss, stabilizing body temperature, and relaying sensory information to

the brain. Maintaining skin homeostasis is essential, as alterations in skin functions can cause various deleterious conditions ranging from fluid loss to more severe diseases such as infections or UV-induced cancers. Because the skin is frequently wounded,

<footnote n="1">¹Université Paris Cité, CNRS, INSERM, Institut Cochin, Paris, France</footnote>

<footnote n="2">²Institut de Biologie Paris Seine (IBPS), Laboratoire de Biologie du Développement, Sorbonne Université, CNRS, INSERM, Paris, France</footnote>

<footnote n="3">³Fondation pour la Recherche en Physiologie, Brussels, Belgium</footnote>

<footnote n="4">⁴Queen Mary University of London, School of Biological and Chemical Sciences, London, United Kingdom</footnote>

<footnote n="5">⁵Service d'Anatomie et Cytologie Pathologiques, Hôpital Tenon, Assistance Publique-Hôpitaux de Paris, Paris, France</footnote>

<footnote n="6">⁶Université Paris Cité, CNRS, INSERM, Institut Necker-Enfants Malades, Paris, France</footnote>

<footnote n="7">⁷Université de Namur ASBL, Namur, Belgium</footnote>

<footnote n="8">⁸Ecole Nationale Vétérinaire d'Alfort, Centre de Recherche Biomédicale, Maisons-Alfort, France</footnote>

<footnote n="9">⁹Service de Dermatologie, Hôpital Cochin, Assistance Publique-Hôpitaux de Paris, France</footnote>

<footnote>*Senior co-author</footnote>
maintenance of an adequate protection requests rapid and finely tuned cutaneous repair steps. In addition to epidermal stem cells committed to the formation, differentiation and regeneration of the epidermis, the skin harbors an elaborate network of immune cells including subsets of T cells but also a substantial component of myeloid-derived cells such as dendritic cells and Langerhans cells (LC). All these specialized immune cells within the skin provides the first line of defense, protecting the body against infections [1, 2].

Constantly exposed to both internal and environmental stresses such as UV radiation or air pollutants, the skin ages and undergo profound changes in its appearance and functions. Indeed, aged skin undergoes gradual structural and functional degeneration, leading to thinning of epidermal and dermal layers, loss of elasticity, wrinkling and dryness. These changes are responsible for delayed wounding, more frequent infections, pruritus, enhanced allergen/irritant penetrations with variable degree of dermatitis and eventually carcinogenesis. In rodents and humans, these phenomena have been partially attributed to loss or lineage skewing of keratinocytes stem cells and immune cells, and/or the regulation of their niches, altering normal homeostasis and tissue repair [3-7].

Aging is often described as an inevitable and irreversible process. However, one species seems to defy the natural laws of aging: the naked mole-rat (NMR, *Heterocephalus glaber*). The NMR are small poikilothermic and hairless rodents native to East Africa, where they live underground in eusocial colonies. These mouse-sized rodents live almost five times longer than expected on the basis of body size, with a maximum lifespan exceeding 37 years in captivity and up to 17 years in their natural habitat [8-11]. Despite being the longest-lived rodent, NMR do not show any increase in age-specific hazard of mortality in defiance of Gompertz's law [12] and all of the classical signs of aging such as decreased fertility, muscle atrophy, bone loss, changes in body composition or metabolism seem to be mostly absent in these animals [13]. Moreover, the incidence of age-related diseases such as metabolic disorders, neurodegenerative diseases and even cancers are absent or extremely low in NMR [14, 15]. The NMR genome has been sequenced in 2011 revealing a unique genetic profile and molecular adaptations correlated with all the characteristics described above: resistance to cancer, poikilothermia, absence of pelage, resistance to low levels of oxygen, as well as impaired vision, taste and circadian rhythms. 93% of the NMR genome has a synteny with the human, rat or mouse genome [16]. Thereby, the NMR, as an example of successful slow-

aging mammal, emerges as a highly attractive organism for the study of skin aging. Investigation of the NMR skin has revealed unique adaptations to thermoregulation and their extreme environmental challenges. The skin morphology and histology of the NMR have been studied, however NMR skin biology during aging is still unexplored [17-20].

We used single-cell RNA-sequencing (scRNA-seq), to obtain the unbiased molecular RNA profile of the NMR epidermal cell populations. By profiling 10,000+ cells from skin epidermis in young and older NMR, we found that epidermal compartments and cell populations, especially the stem cells pool, remained unaffected despite aging. *Igf/bp3*, expressed by keratinocyte stem cells and known to play a major role during epidermal homeostasis, was found upregulated in older animals, contrary to what is observed in other species. In addition, functional skin healing experiments revealed that the wound closure rate was strictly comparable between young and older animals.

## RESULTS

### Histomorphology of young versus middle-aged NMR skin

First, to investigate NMR skin structure and organization changes associated with cutaneous aging, dorsal skin of young (Y, mean age: 1.1 years) and middle-aged (MA, mean age: 11.3 years) animals were collected and evaluated by immunohistochemistry. HE and Sirius Red staining were used to measure epidermal and dermal thickness of NMR skin (Figure 1A, 1B). In human and mouse, skin thickness decreases with age. Surprisingly, NMR epidermis was significantly thicker in older animals (histogram in Figure 1A), and the global histomorphological structure was unchanged during aging. This increased epidermal thickness was confirmed in older 20-25 years NMR tissue samples (only available for histology analysis, Supplementary Figure 1A-1C). The number of epidermal cells (Supplementary Figure 1D) and epidermal layers (Supplementary Figure 1E) were significantly higher in the MA group. More precisely, cell counting per layer revealed an additional cellular layer ("layer 6") in older animals (Supplementary Figure 1F). No effect of age was observed on the dermal thickness between Y and MA animals (histogram in Figure 1B). Masson's trichrome staining for collagenous connective tissue fibers showed no difference between the 2 groups just like hyaluronic acid, a linear polysaccharide of the extracellular matrix contributing to contact inhibition and cell proliferation [21] or collagen I (found in all dermal layers; Figure 1C-1E). Epidermal layers immunostaining using anti-Keratin-14 (Krt14, basal
layer marker), Keratin-10 (Krt10, suprabasal layer marker), Loricrin (Lor, a major protein of the cell corneous envelope), Filaggrin (intermediate filament-associated protein of the outer granular layer of the epidermis), Laminin-5 (dermal-epidermal basement membrane protein), Integrin-beta4 (ITGB4, basal surface marker juxtaposed to the dermal-epidermal basement membrane) antibodies showed no difference between Y and MA animals (Figure 1F–1K). It has been discovered that epidermal proliferation declined with

age especially in humans [4]. Interestingly, in NMR, Ki67 index for cell proliferation revealed no significant difference between Y and MA NMR (Figure 1L). p16 (or CDKN2A), one of the most frequently used markers of senescence and an important hallmark of aging, was slightly detected in NMR skin but again, no differences were observed between the 2 age groups (Figure 1M) whereas this marker accumulates in an age-dependent manner in mammalian skin [22]. Thus, the results suggest that NMR skin aging is delayed.

<figure><img src="image_2.png" /><figcaption>**Figure 1. Naked mole-rat (NMR) skin histomorphology.** (A–C) HE, Sirius Red and Masson's Trichrome staining photomicrographs on sections of middle-aged and young NMR back skin. Histograms in A and B show the epidermal and dermal thickness (µm) of middle-aged and young NMR back skin. Scale bar = 50 or 100 µm. Bars: SEM. *represents differences between the age groups. ****p≤0.0001 in Student's t-test. n=4 animals per group. (D–M) Hyaluronic acid, Collagen I, Keratin-14 (Krt14), Keratin-10 (Krt10), Loricrin (Lor), Filaggrin, Laminin-5, ITGB4, Ki67 and p16 staining photomicrographs on sections of middle-aged and young NMR back skin. Cells were stained with either Alexa 546 or Alexa 488 fluorochromes. For hyaluronic acid, cells were stained with 3,3'-diaminobenzidine. Negative controls were performed in parallel with the samples substituting the primary antibody with the equivalent isotype. Histograms show the quantification of Krt14 and Krt10 expressions and Ki67+ cells per surface in both age group. Scale bar = 50 or 100 µm. Bars: SEM. A Mann–Whitney U test was used for Keratin-14 / Keratin-10 quantification and a Student's t-test was used for Ki67+ cells quantification. n=4 animals per group.</figcaption></figure>
scRNA-seq of NMR epidermis defines 3 cellular states for NMR keratinocytes

The skin epidermis is constantly rejuvenated by keratinocyte stem cells found in the basal layer that are capable of both proliferation and differentiation. But aging skin experiences a progressive decline in homeostatic and regenerative capacities, potentially attributed to degenerative changes in these stem cells and/or their niches [3-7]. Using single-cell RNAseq, we analyzed the cellular composition of the NMR epidermis to investigate the potential link with their slower rate of skin aging.

After a short enzymatic epidermis-dermis digestion and separation, epidermal cells were mechanically dissociated and directly loaded into the chromium controller (10x Genomics) without any FACS-sorting to prevent cells from stress (Figure 2A). All the samples were processed in parallel to minimize batch effects. Quality control (# UMI ≥ 2000, # genes ≥ 700, % MT ≤ 25%, number of cells) was made with FastQC, Cellranger and Seurat v3 pipeline (Supplementary Figure 2A-2D). Unsupervised clustering revealed that one sample from the MA group had a very distinctive pattern after Principal Component Analysis (PCA). UMAP (Uniform Manifold Approximation and Projection) visualization for all individual samples is presented in Supplementary Figure 2E. This sample (2V) was removed, keeping for the rest of our study 3 samples in the Y NMR group and 2 samples in the MA NMR group (see materials and methods). In total, we obtained 10,232 cells from 5 NMR epidermis for downstream analyses. Based on gene expression profiles and UMAP to visualize cell clusters, we found 11 cellular populations (numbered from 0 to 10, Figure 2B). To identify marker genes of each cluster, we used differential expression testing, with a Wilcoxon rank sum test. These marker genes were used to establish the cell identity of each cluster, together with additional markers identified through literature searching for cell types typically found in the human or mouse skin and mining existing single-cell transcriptomic skin data (Supplementary Data File 1). We found 2 different populations: keratinocytes representing 99,4% (10,171 cells) in clusters 0 to 9, and immune cells in cluster 10 representing 0.6% of epidermal cells (61 cells). This immune-isolated cluster 10 contained a majority of Cd207⁺/Ctss⁺/Mfge8⁺ LC [23, 24] (Supplementary Data File 1). We did not find cells expressing melanocytes genes, which is consistent with previous reports showing that melanocytes were only detected in the NMR dermis [17].

Cytokeratin expression can globally define the differentiation status of keratinocytes: *Krt5* and *Krt14*

are expressed in undifferentiated basal stem and proliferative cells, while *Krt1* and *Krt10* are found in differentiating suprabasal cells and *Lor* occupies 80% of the epidermal corneous envelope [25-27]. Thus, using these well-defined marker genes, we were able to confirm that in NMR, the keratinocyte population also separates into 3 main cellular states. Four clusters (clusters 1-3-6-7) expressed high levels of *Krt14* (*Krt14*<sup>high</sup>) and *Krt5* and were classified as basal cells. Clusters 0-2-4-5-9 expressed higher levels of *Krt10* (*Krt10*<sup>high</sup>) and *Krt1* and were categorized as suprabasal or intermediate cells. The last keratinocyte-cluster (cluster 8) was the only one to express very high levels of *Lor* (*Lor*<sup>high</sup>) and was considered as a cluster of terminally differentiated keratinized cells of the corneous layer (Figure 2C). The same results were confirmed with violin plots of these marker genes (Supplementary Figure 3A) and t-distributed Stochastic Neighbor Embedding (t-SNE) visualization (Supplementary Figure 3B, 3C). The most differentially expressed marker genes for each cluster were represented using heatmap visualization (Figure 2D) and 3 specific genes for each cluster were highlighted and color-coded to match cluster color. For each cluster, the most representative gene was projected onto the UMAP plot (Figure 2E).

In order to confirm these 3 cellular states in NMR keratinocytes, we identified other well-known skin layer genes in each cluster. In *Krt14*<sup>high</sup> basal cells clusters, cells from clusters 1-6-7 expressed specific basal layer genes such as *Tp73*, *Itga6*, *Itgb1*, *Cxcl14* and *Smoc2* [27-30]. Cells from cluster 1 expressed the Notch ligand *Dll1*, a marker of stem cells found in the basal layer in human and mouse skin [28, 31], *Igfbp3*, a gene exclusively expressed in the basal layer's proliferative keratinocytes in normal skin [4, 28, 32, 33], and *Socs3*, a gene known to promote maintenance of epidermal homeostasis [28]. Of note, *Dll1* and *Igfbp3* were also expressed in the two other basal clusters 6 and 7. Cluster 6 cells expressed *Bcam*, *Ascl2* and *Rarres2*, 3 genes found in the basal layer of the epidermis and known to stimulate keratinization [34, 35]. *Tsc22d1*, *Tnc*, *S100* family members (such as *S100a4* and *S100p*) and *Ccl2* expressed in Cluster 7 cells were reported to be restricted to basal/follicular stem cells and progenitors with high cell turnover notably during wound healing [28, 33, 36]. *Krt14*<sup>high</sup> cluster 3 cells specifically expressed classical epidermal cell cycle regulators such as *Cdc20*, *Ube2C* and *Cdkn3* [33, 37, 38] (Figure 2E). Concerning *Krt10*<sup>high</sup> suprabasal cells clusters 0-2-4-5-9, we found that *Notch* receptor gene was expressed in clusters 5 and 9 cells, while its target gene *Hes1* was found in cluster 0 and 4. It has been demonstrated that *Notch* activation induces terminal epidermal differentiation through *Hes1* which
A

(1) Skin biopsies of
young and aged NMR

<figure><img src="image_3.png" /></figure>

(2) Keratinocytes
viability ≥ 75%

<figure><img src="image_4.png" /></figure>

(3) Single cell partitioning;
Barcoding; Reverse transcription

<figure><img src="image_5.png" /></figure>

(4) Library
construction

<figure><img src="image_6.png" /></figure>

(5) Data analysis
& visualization

<figure><img src="image_7.png" /></figure>

B

Clusters - All cells

<figure><img src="image_8.png" /></figure>

C

<figure><img src="image_9.png" /></figure>

Krt10<sup>high</sup>
suprabasal cells

<figure><img src="image_10.png" /></figure>

Lor<sup>high</sup>
terminally
differentiated cells

<figure><img src="image_11.png" /></figure>

D

<figure><img src="image_12.png" /></figure>

E

BSC1 - Bcam (9.7%)

<figure><img src="image_13.png" /></figure>

BSC2 - Tsc22d1 (9.5%)

<figure><img src="image_14.png" /></figure>

BSC3 - Dll1 (15.8%)

<figure><img src="image_15.png" /></figure>

BPSC - Cdc20 (11.5%)

<figure><img src="image_16.png" /></figure>

Sp1 - F3 (15.8%)

<figure><img src="image_17.png" /></figure>

Sp2 - Endou (10.1%)

<figure><img src="image_18.png" /></figure>

Gr1 - Slurp1 (12.1%)

<figure><img src="image_19.png" /></figure>

Gr2 - Ccnf (9.8%)

<figure><img src="image_20.png" /></figure>

Gr3 - Cldn1 (2.3%)

<figure><img src="image_21.png" /></figure>

Cor - Krt23 (2.7%)

<figure><img src="image_22.png" /></figure>

Immune - Cd207 (0.6%)

<figure><img src="image_23.png" /></figure>

**Figure 2.** Cell populations characterization of NMR epidermis using single-cell RNA-sequencing (scRNA-seq). (A) Overview of the experimental workflow. (B) Global epidermal cell transcriptomes (n=10,232 cells) visualized with Uniform Manifold Approximation and Projection (UMAP), colored according to unsupervised clustering to better compare different cell subtypes. 11 different cell clusters were found, numbered from 0 to 10. (C) Expression of the 3 main marker genes during epidermal differentiation (Krt14, Krt10, Lor) projected onto UMAP. Cells were sub-divided into 3 main groups according to the differential epidermal cellular states: basal (Krt14<sup>high</sup>), suprabasal/intermediate (Krt10<sup>high</sup>) and corneous (Lor<sup>+</sup>) layers. Cells with the highest expression level are colored dark blue. (D) Heatmap of differentially expressed genes. For each cluster, the most differentially expressed genes and their relative expression levels in all sequenced epidermal cells are shown. Cells are represented in columns, and genes in rows. 3 selected genes for each cluster were color-coded and annotated on the right. The 10 clusters were identified as: basal and stem cells BSC1-2-3 for clusters 6-7-1; basal proliferating and stem cells
BPSC for cluster 3; spinous layer cells Sp1 and Sp2 for clusters 0-4; granular layer cells Gr1-2-3 for clusters 2-5-9 respectively; corneous layer cells Cor for cluster 8; immune for immune cells cluster 10. (E) Expression levels for each cell are color-coded and overlaid onto UMAP plot for one selected cluster-specific gene. Cells with the highest expression level are colored dark blue. % of cells in each cluster is annotated.

is expressed in the spinous layer. The expression of another Notch ligand, *Jag2* was restricted to basal layer's clusters 1-6-7 cells corroborating previous findings. *Jag1* on the other hand was only found in suprabasal cluster 9 cells [31]. In addition to *Hes1*, cluster 0 cells also expressed *F3* and *Ehf* reported to be located in suprabasal layers of the skin playing an essential role in keratinocyte differentiation [28, 39]. Cluster 4 analysis revealed the expression of *Endou*, *Ly6D* and *Mt4*, 3 genes previously associated with keratinocyte differentiation in the spinous layer of the skin [28, 30]. Of note, *Mt4*, reported to specify a transitory stage from basal to spinous layer [28, 36], was expressed in cells from cluster 0, 4 and 2. Cluster 2 cells strongly showed classical suprabasal gene expression such as *Krt1*, *Krt10* and *Slurp1*, known to be a marker of late differentiation state, predominantly expressed in the granular layer of the skin [40]. In cluster 5 cells, *Cnfn*, *Ly6g6c* and *Clic3* expression matched with a granular layer phenotype [28]. The last suprabasal cluster defined here, cluster 9, was also consistent with genes expressed by differentiating keratinocytes in the granular layer of human and mouse skins such as *Cldn1*, *SerpinB12* and *Lypd3* [28]. In addition, both clusters 5 and 9 cells revealed the expression of several strong makers such as *Hspb1* and *Grhl1* suggesting the possible granular keratinocytes phenotype of these two cluster cells [28, 38] (Figure 2E). Finally, in addition to *Lor*, cluster 8 cells expressed *Klk5/7*, *Krt23*, *Cdn* or *Lce5a* strongly suggesting these cells were terminally differentiated keratinocytes from the corneous layer [28] (Figure 2E, Supplementary Data File 1). We were then able to label these keratinocytes clusters. Undifferentiated cells from Clusters 1-6-7 were named basal and stem cells BSC3, BSC1 and BSC2 respectively. Cluster 3 was named basal proliferating and stem cells BPSC. Cluster 0 and 4 were defined as spinous layer cells Sp1 and Sp2, cluster 2-5-9 as granular layer cells Gr1-2-3 respectively. At last, cluster 8 was labeled as corneous layer cells Cor (Figure 2D, 2E).

The single cell transcriptomic analysis of NMR keratinocyte unraveled, for the first time in NMR, an epidermal differentiation process that follows 3 classical cellular states, from basal layer's proliferative stem cells and/or progenitors (BSC and PPSC clusters) to suprabasal (Sp1-2 and Gr1-2-3 clusters) and terminally differentiated cells (Cor cluster). Within each of these three cell population, each cluster could represent a specialized subpopulation.

## Cell populations and trajectories analysis of NMR keratinocytes define metaclusters but reveal no age-related differences

To study any transcriptional changes of the NMR epidermis during aging, we compared Y and MA NMR epidermis. Three Y NMR of 1.1 years old (n=6,014 cells) and two MA NMR of 11.3 years-old (n=4,218 cells) were analyzed. Y and MA cells seemed to be equally spread across clusters using UMAP visualization (Figure 3A, 3B). This observation was confirmed using t-SNE visualization (Supplementary Figure 4A). In addition, the clusters' cellular states, defined by the expression of *Krt14*, *Krt10* and *Lor* marker genes and visualized with violin plots, were identical between the two age groups (Figure 3C). The proportion of cell types in each cluster after normalization was not statistically different between the 2 conditions (Chi² statistic test pval ≤ 0,05; Figure 3D). To identify differentially expressed genes between Y and MA animals, we then compared the expression of genes between all cells from Y samples and all cells from MA samples. No gene was found differentially expressed using the following thresholds: p-value below 0.05 and average log fold-change (logFc) above 1. This result was surprising because it is now well known that aging induces the differential regulation of multiple growth-promoting genes. For example, 1237 genes were found differentially expressed in the skin of mice at 9 months, 15 months, and 24 months [41].

As keratinocytes differentiate, they undergo a process of transcriptional reconfiguration. Using Monocle 2, we reconstructed the single cell trajectory of Y and MA NMR keratinocytes to mimic their biological differentiation process. Of note, cluster 10 cells corresponding to epidermal immune cells was excluded from all Monocle analysis to focus only on the keratinocytes' differentiation program. Colored by cluster, we found the same linear trajectory for each age condition, with no branches, suggesting an absence of different cellular decision or alternative expression program between Y and MA keratinocytes (Figure 3E). For both Y and MA NMR, the distribution of keratinocytes along their process of cellular differentiation showed that, as expected, *Krt14*<sup>high</sup> clusters BSC1-2-3 were condensed at the root of the trajectory process, with cluster BPSC is spreading out more from this root showing their early engagement in the differentiation process. Then, intermediate *Krt10*<sup>high</sup>
A

Clusters - MA cells

<figure><img src="image_24.png" /></figure>

Immune

Clusters - Y cells

<figure><img src="image_25.png" /></figure>

Immune

B

Age - All cells

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

MA

Y

C

Krt14

<figure><img src="image_28.png" /></figure>

Krt10

<figure><img src="image_29.png" /></figure>

Lor

<figure><img src="image_30.png" /></figure>

Krt14

<figure><img src="image_31.png" /></figure>

Krt10

<figure><img src="image_32.png" /></figure>

Lor

<figure><img src="image_33.png" /></figure>

D

<figure><img src="image_34.png" /></figure>

E

<figure><img src="image_35.png" /></figure>

Clusters - MA cells

<figure><img src="image_36.png" /></figure>

Clusters - Y cells

<figure><img src="image_37.png" /></figure>

**Figure 3.** Single cell populations and trajectories analysis of young versus middle-aged NMR epidermal cells reveals no difference between the 2 age groups. (A) Uniform Manifold Approximation and Projection (UMAP), colored according to unsupervised clustering from middle-aged (MA; n=4,218 cells) and young (Y; n=6,014 cells) NMR epidermal cells. Cell populations from both MA and Y NMR were unchanged as the same 11 clusters were detected. (B) Cells from both MA (red) and Y (blue) NMR were jointly projected on the same UMAP plot, showing an overlap of the 2 age groups tested. (C) Violin plots of *Krt14*, *Krt10* and *Lor* marker genes expressed by MA (red) and Y (blue) epidermal cells (upper panel) and segregated by cluster (lower panel). No detectable changes in their expression were found. (D) Bar graph representing the relative proportion of epidermal cells in each cluster between MA (red) and Y (blue) animals. There was no significant difference in epidermal cell proportion between the 2 age groups using Chi² statistic test (pval ≤ 0,05). (E) Unsupervised differentiation trajectories constructed with M3Drop and Monocle v2.10.1 for MA and Y keratinocytes. Immune cells' cluster 10 has been excluded to focus on epidermal keratinocytes. The results showed a unique trajectory with no branches in both samples and the same repartition of cluster among the trajectory. BSC = basal and stem cells; BPSC = basal proliferating and stem cells; Sp = Spinous layer cells; Gr = granular layer cells; Cor = corneous layer cells.
clusters Sp1-2 and Gr1-2-3 cells moved towards the right of the trajectory, as they were more committed in the differentiation process (no cells were observable in the root). And finally, *Lor*<sup>high</sup> Cor cluster cells were distributed at the very end of the trajectory, for terminally differentiated keratinocytes [42] (Figure 3E). These results were also represented by Monocle states (state or “segment of the tree”, Supplementary Figure 4B) and for each cluster (on separate plots) for merged samples or Y and MA samples separately, to see precisely where each cluster was located (Supplementary Figure 4C).

Colored by pseudotime, the cells were also visualized according to a blue gradient which becomes lighter as cells move away from the root of the trajectory, reflecting their maturation state. No significant differences on pseudotime estimation were found between Y and MA NMR (Figure 4A). To validate our 3 cellular states in NMR keratinocytes, the 3 selected well-known gene markers of epidermal differentiation *Krt14*, *Krt10* and *Lor* [42] were plotted according to their mean expression levels (black line) along the pseudotime axis (Figure 4B). No changes between the 2 age groups were observed suggesting an absence of age-related effect on keratinocyte differentiation in NMR (Figure 4B). To assign a more precise cellular state of differentiation to the 10 clusters, we chose to individualize them on separate plots (Figure 4C). Thus, combined with data relative to specific gene expression obtained in Figure 2, pseudotime reconstruction allowed us to order the 10 clusters precisely and to partition them into 5 metaclusters (i.e group of clusters) according to their differentiation engagement: i) basal 1 containing BSC1-2-3 clusters; ii) basal 2 for BPSC cluster; iii) intermediate 1 containing Sp1-2 clusters; iv) intermediate 2 aggregating Gr1-2-3 clusters and finally v) terminally differentiated cells for the Cor cluster. The expression of *Krt14*, *Krt10* and *Lor* marker genes visualized with violin plots confirmed the cellular states of these newly defined metaclusters (Supplementary Figure 5A). The repartition of cell types between the Y and MA animals after normalization was not statistically significant in each metacluster (Chi² statistic test pval ≤ 0,05; Supplementary Figure 5B). Distributions of the frequency of Y and MA cells as a function of the pseudotime on its scale were also plotted. Wilcoxon test revealed that Y and MA NMR keratinocytes distributions were similar using the pseudotime estimation (p-value = 0.1673). According to pseudotime segments, no different gene markers were found between the segments (Supplementary Figure 5C). Then, we observed the most differentially expressed genes and their relative expression levels in all sequenced epidermal cells on a heatmap (Supplementary Figure 5D). Similar to the previous

analysis, five metaclusters of genes could be clearly identified, each with a different region of activity during the differentiation process. Unsupervised differentiation trajectories for Y and MA animals showed identical repartition of these metaclusters among the trajectories (Figure 4D).

To conclude, we identified 11 different clusters defining epidermal cells from Y and MA NMR and a single path cell trajectory. Using pseudotemporal ordering, we were able to identify 5 metaclusters. Basal metaclusters containing keratinocytes stem cells and proliferative progenitors seemed to contribute to all the other cellular states to promote epidermal differentiation, up to the upper corneous layer. Taken all the results together, comparing young and older NMR keratinocytes, we observed a well-maintained skin tissue and cellular identities throughout the years, contrary to what has been described in mouse or human [4, 27].

## Metaclusters characterization and enrichment analysis for specific functions of NMR keratinocytes

As basal keratinocytes transit through different cellular states until terminal differentiation, we examined the dynamics of the genetic profiles used to construct the trajectories using Monocle. We identified 6 gene modules following similar kinetic gene changes along pseudotime. The heatmap profiles were identical, suggesting that the same gene groups at the same expression level were orchestrating the differentiation process in both Y and MA NMR keratinocytes, independently of aging (Supplementary Figure 6A, left panel). We then investigated whether specific functions could be assigned to those 6 groups using Gene Ontology (GO). The most significant GO terms found in each activated gene group for both ages comforted their implication in their commitment : early activated cluster 5 for basal layer state, cluster 1 for cell cycle still in the early pseudo-timeline (note the sharp transition between these 2 states), mid-time activated cluster 6 for cell differentiation, cluster 3 at the end of differentiation program for an ultimate transition, cluster 2 strictly circumscribe at the end of pseudotime differentiation process and finally cluster 4 harbouring a biphasic expression pattern with functions involved at various times of the differentiation process supported by different genes (Supplementary Figure 6A, right panel).

Next, another GO analysis was performed on the different epidermal metaclusters defined above. We used the most representative markers of each metacluster and described a repartition of function for each of them (Supplementary Figure 6B; Supplementary Data File 2). Basal 1 metacluster (BSC1-2-3 clusters) was characterized by developmental and cell adhesion
<figure><img src="image_38.png" /><figcaption>**Figure 4.** Pseudotime analysis and metacluster identification in young versus middle-aged NMR keratinocytes. (A) Colored pseudotemporal ordering of keratinocytes for MA (middle-aged) and Y (young) NMR separately. No significant differences were observed between the 2 groups. (B) Relative expression of *Krt14*, *Krt10* and *Lor* markers plotted along the pseudotime axis for each condition. The same relative expression profile was found between MA and Y NMR, validating the 3 cellular states for NMR keratinocytes and the no age-related impact of aging in the differentiation process. (C) Pseudotime estimation showing epidermal progression along the differentiation process with individualized cluster visualization, for each condition. Five metaclusters were defined according to their differentiation engagement: basal 1, basal 2, intermediate 1, intermediate 2 and corneum. (D) Unsupervised differentiation trajectories for MA and Y animals. The same repartition of metacluster among the trajectory was observed. BSC = basal and stem cells; BPSC = basal proliferating and stem cells; Sp = Spinous layer cells; Gr = granular layer cells; Cor = corneous layer cells.</figcaption></figure>
processes relative to their basal stem cell/progenitor phenotypes. Indeed, those cells were likely prone to establish cell junction to the extracellular matrix (ECM) and to commit to a global developmental process of tissue formation. As expected, basal 2 metacluster (BPSC cluster) was well defined by a wide range of cell cycle functions relative to their mitotically active progenitor cell type. Intermediate 1 (Sp1-2 clusters) and intermediate 2 (Gr1-2-3 clusters) metaclusters were characterized by cell differentiation and developmental processes confirming these clusters as a population of differentiated keratinocytes of the spinous and granular layers respectively. Lastly, corneum metacluster (Cor cluster) was represented by active lipid metabolism as well as cell death regulation and reactive oxygen species (ROS) reaction terms (Supplementary Figure 6B; Supplementary Data File 2). The characteristics of metaclusters were well maintained in both age groups.

Furthermore, due to the link between ROS and aging, we performed an enrichment analysis for our metaclusters for aging and ROS-related genes. Surprisingly, none of the aging-related gene analyzed were found enriched in any of our metaclusters. However, 4 ROS-related genes were significantly enriched in corneum metacluster (Supplementary Figure 6C and Supplementary Data File 3). Indeed, *Prdx6*, *P4hb*, *Sod1* and *Prdx2* were enriched in the Corneum metacluster (Supplementary Figure 6C and Supplementary Data File 3; pval 1,01e-02). In NMR, the abundantly expressed antioxidant enzyme *Sod1* activity was 1.35-fold higher as compared to age-matched mice but age did not have an effect on its activity in both species. Peroxiredoxins, especially *Prdx2* and *Prdx5*, were shown to be expressed at lower levels in NMR livers, which could result in increased levels of ROS and suggest that the long-lived NMR can thrive despite elevated oxidative stress [43]. Of note, *Prdx 2* and *5* were also enriched in the immune cluster. Those immune cells expressed high levels of transcripts related to regulation of immune system process and response to stimulus, showing that their main function is to be involved in resistance to skin pathogens. Those terms were also found at a lower level in the corneum metacluster where immunity meets external antigens, and in the basal 1 metacluster where keratinocytes have key roles in immune defense (Supplementary Figure 6D and Supplementary Data File 3).

### Igfbp3 is overexpressed in older NMR skin

Previous studies have demonstrated that cutaneous gene expression levels change with age, but NMR cutaneous gene expression showed a remarkable stability with age (see above). However, as subtle differences might exist, we performed a more in-depth analysis between Y and

MA animals within each cluster independently. Using a p-value below 0.05 and an average logFc above 1, only 2 differentially expressed genes in BSC2 cluster (*ENSHGLG00000002542/Rps2* and *Igfbp3*) were overexpressed in MA NMR (Supplementary Data File 4). *Rps2* gene encodes a ribosomal protein that is a component of the 40S subunit, with unknown role in skin homeostasis. Interestingly, *Igf/Igfbp* signaling plays an important role in aging processes and mediates epidermal proliferation and skin peripheral immunity. *Igfbp3*, one of six IGF binding proteins, was found located in the basal/germinative layer of the epidermis and modulates epidermal homeostasis. Here, using UMAP visualization, we retrieved its expression in basal stem and proliferative clusters BSC1-2-3 and BPSC (Figure 5A). In mice and humans, *Igfbp3* transcript is significantly reduced with increasing age [4, 28, 32, 33]. In NMR, *Igfbp3* expression was almost 4 time higher in the BSC2 cluster of older animals as compared to the younger ones (Figure 5B, 5C and Supplementary Data File 4). Confirming our scRNA-seq data, IGFBP3 nuclear expression assessed by immunohistochemistry was increased in older animals (Figure 5D). Moreover, 5 differentially expressed genes in immune cluster (*S100a9*, *ENSHGLG00000004445/S100a8*, *Krt17*, *S100p*, and *Krt5*) were overexpressed in older animals. Interestingly, no differentially expressed genes typical of LC were found, suggesting the stability of these cells over time in the NMR (Supplementary Data File 4).

### Wound healing is not affected by age in NMR skin

The ability of the skin to heal quickly and completely is a major function that characterizes the skin status. Epidermal stem cells, and their crucial role in the cutaneous healing process, have been extensively studied in mice and humans [33, 44] but never in NMR. Stem cell dysfunction is a hallmark of aging and is associated with impaired wound healing in older rodents and humans [5, 7, 45]. The delayed skin aging in the NMR prompted us to investigate the wound healing process in young and older NMR. Punch biopsy devices were used to create surgical wounds on the backs of Y and MA NMR and the wound closure rate was measured (percentage of wound closure at day 0) for the first time in these animals. Unlike normal mice that possess a healing time course of approximately 5 days, NMR appeared to heal at a slower rate of 22 days but with less scarring (Figure 6A, 6B). In humans and mice, wound closure occurs more rapidly in younger than older individuals [5, 7, 45]. Our results showed that in NMR, despite an age difference of 11 years, the wound closure rate of MA animals followed the same kinetic as the one of Y animals (Figure 6B). Next, we investigated re-epithelialization at day 7 post-wounding
as an essential component of wound healing and used as
a defining parameter of its success (Supplementary
Figure 7A). A well-organized neo-epidermis of the
same thickness between Y and MA was measured
(Figure 6C and Supplementary Figure 7B). In all
individuals, the same levels of expression of basal

Krt14 and differentiated Krt10 cell markers were found
and the neo-epidermal tongue expanse was substantially
identical between the 2 groups (Figure 6D, 6E). The
epidermal tongue proliferation rate measured using
Ki67 immunostaining remained stable between Y and
MA animals (Figure 6F).

<figure><img src="image_39.png" /><figcaption>**Figure 5.** *Igfbp3* expression in young versus middle-aged NMR keratinocytes. (A) Expression levels for *Igfbp3* gene onto Uniform Manifold Approximation and Projection (UMAP). Note that its expression is restricted to basal stem and proliferative cell clusters BSC1-2-3 and BPSC. (B) Expression levels for *Igfbp3* gene onto UMAP for MA (middle-aged) and Y (young) NMR. The dotted circle delimits BSC2 cluster. (C) Violin plots of *Igfbp3* genes expressed by MA (red) and Y (blue) epidermal cells. (D) IGFBP3 staining photomicrographs on sections of MA and Y NMR back skin. Cells were stained with Alexa 546 fluorochromes. Negative controls were performed in parallel with the samples substituting the primary antibody with the equivalent isotype. Scale bar = 50 µm. Histogram showing the quantification of nuclear IGFBP3 expression in both age group. Bars: SEM. \*represents differences between the age groups. \*\*\*\*p≤0.0001 in Student's t-test. n=4 animals per group.</figcaption></figure>
<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /><figcaption>**Figure 6.** NMR wound healing. (A) Photomicrographs of back skin wounds for middle-aged and young animals from day 0 to day 22 post wounding. (B) Wound closure rate (% of wound surface at day 0) from day 4 to day 22 post wounding. No statistical difference was found in the 2 groups using Student's t-test. n=4 animals per group. (C) HE staining photomicrographs on sections of middle-aged and young NMR back skin 7 days post wounding. (D, E) Krt14 and Krt10 staining photomicrographs on sections of middle-aged and young NMR neoepidermal tongues in the wound at day 7. Histograms show the level of expression of Keratin-14 and Keratin-10 in the neoepidermis of wounded samples (mean grey value per field) and the percentage of re-epithelialization was calculated according to the formula: $[length of the extending epidermal tongues]/[length of the wound] \times 100\%$. (F) Ki67 staining photomicrographs on sections of middle-aged and young NMR neoepidermal tongues in the wound at day 7. Histogram show the percentage of Ki67+ cells in epidermal wound edges representing the growth fraction of an epidermal tongue cell population. Legends: epidermal wound edges (EpiD), dermis (D), fibrino-leukocyte exudate (FLE). In (D–F), cells were stained with either Alexa 546 or Alexa 488 fluorochromes. Negative controls were performed in parallel with the samples substituting the primary antibody with the equivalent isotype. Scale bar = 50 or 100 µm. Dash lines represent the limit between epidermis and dermis. Bars: SEM. Student's t-test were used. n=2 animals per group.</figcaption></figure>
Metaclusters enrichment analysis for wound healing related genes revealed that 3 genes were significantly enriched from a list of 74 genes in Basal 1 metacluster: *Cav1*, *Dst*, *Esp* (Supplementary Figure 6C and Supplementary Data File 3; pval 5,86e-03). Mutations in the mouse *dystonin* (*Dst*) and *desmoplakin* (*Esp*) genes have previously been shown to cause skin fragility [46, 47]. *Cav1* has been found in undifferentiated epidermal basal layer keratinocytes and considered by others as a potential epidermal stem cell marker [48]. Moreover, some studies have proven that disruption of *Cav1* expression can lead to age-related changes in skin and aberrant wound healing [49]. Interestingly, *Igfbp-3* can be internalized independent of IGFs into the cell through caveolae-mediated endocytosis.

DISCUSSION

Age-related epidermal changes can be caused by
multiple elements leading to: reduction of keratinocyte
stem cell compartment, alteration of cellular
proliferation, thinning of the epidermis, slowing of the
keratinocytes turnover, flattening of the dermal-
epidermal junction or the regulation of the stem cell
niche. The results of our investigation clearly showed
that NMR skin ages differently than rodent or human
skin [3, 4, 27, 38, 50, 51]. First, skin structure and
organization appear to be well maintained between
young and older NMR. Second, using scRNA-seq, we
found that their epidermal cells displayed 2 main
populations, immune cells (1 cluster) and
keratinocytes, subdivided into 10 clusters. The immune
cells cluster was easily distinguished from the other
cell clusters with a characteristic signature (*Cd207*,
*Ctss*, *Mfge8*) corresponding mainly to LC. By ordering
the 10 keratinocytes clusters, we found a classical
differentiation process including basal stem cells and
proliferative keratinocytes, followed by spinous then
granular layer cells and finally terminally differentiated
cornified cells. The results of the present study suggest
that NMR epidermal homeostasis remain unaltered for
the first 11 years of their life and for some parameters
up to 25 years old. Indeed, epidermal cell identities,
cellular states and differentiation trajectory were highly
similar in older animals. Recently, single cell analysis
was successfully used to detect differences between
young and adult mice epidermis, showing that
young epidermal cells progenitors appeared more
homogenous compared with adult cells progenitors.
Moreover, these authors found only one differentiation
trajectory in young mice but 2 in older mice, while
others have found up to 3 distinct trajectories in adult
mice interfollicular epidermal cells [28, 33, 36]. Using
also single-cell transcriptomic data, The Tabula Muris
Consortium assessed cell-type-specific age-related
changes in *Mus musculus* tissues. In skin, these authors

found that the relative abundance of keratinocyte stem
cells decreased significantly with age with at least 20
differentially modulated genes in skin computed by age
(3, 18, 21 and 24 months) [50]. Another recent article
confirmed that human epidermal thickness decreased
constantly over time, with thinner skin already
observed in MA individuals (35-48 years old). In
addition, unlike the NMR, thousands of differentially
expressed genes were identified in the human aging
skin. Further analysis of these differentially expressed
genes highlighted that MA individuals were closer to
old individuals at the transcriptomic level, indicating
the early onset of human skin aging [51]. These
data could indicate that NMR are better protected
against skin aging through the mainstay of gene
expression, contrary to what has been shown in mice
and humans [41, 51].

Performing a deeper analysis within each cluster
individually, we found 2 genes overexpressed in basal
stem cells of older animals and 5 genes overexpressed
in immune cells of older animals. In basal stem cells,
*Igfbp3* was four times higher in older animals, *while* in
mice and humans, there is a large decrease of this gene
during skin aging [4]. This analysis also revealed that
the level of LC appeared unaffected through aging in
the NMR. While, in mouse and human epidermal aging,
a decrease in the number of LC has been proven; their
numbers and maturation state are reduced in aged mice
starting at 12 months of age [6]. In addition, in NMR
cutaneous immune cells, there were overexpressed
genes belonging to the *S100* and *Keratin intermediate
filament* families, known to be regulators of
inflammation and immunity in the skin. To date, only
few studies using scRNA-seq profiling have focused
their attention on NMR immunity, but no data on skin
immune cells were available. First, authors found that
NMR immune system was characterized by a high
myeloid to lymphoid cell ratio as compared to mice but
only healthy young NMR were used [52]. More
recently, NMR aged 3 and 11 years were used as
compared to mice aged 3 and 12 months. They
discovered that NMR owned an additional pair of
cervical thymi and that they had no thymic involution
up to 11 years of age, suggesting that they displayed a
delayed immunosenescence [53]. Thus, we hypothesize
that the maintenance of cellular compartments in the
older NMR, especially the stem cell pool through high
*Igfbp3* expression, coupled with an increase skin
immunity, could explain their skin slower rate of aging.
In addition, a recent study suggests that NMR skin
could be immune to senescent cells build-up through a
unique mechanism: senescent cell death, an immune-
independent apoptosis of senescent cells. This
mechanism could reduce the chance of senescent cell
accumulation, a hallmark of skin aging [54].
Age-related skin changes have a deleterious impact on water loss, body temperature, and wound healing, making the skin more susceptible to long-standing wounds. Surgical excision wounds performed on the back of the NMR, showed an impressive maintenance of healing kinetics in older animals. A well-organized neo-epidermis with a clear neo-epidermal tongue expanse was found regardless of the age of the animals. Once again, these results are in contrast with what has been shown in mouse and human [5, 7, 45]. For example, wound closure has been proven to occur more rapidly in young mice (6 months) than in middle-aged (15 months) or aged mice (26 months; [45]). This important finding demonstrates that NMR skin remains protected against one of the worst consequences of skin aging.

At this point, several other elements could be discussed. First, the MA animals were the oldest animals currently available in our laboratory. However, 11 year-old NMR are considered old in the wild where only few animals live beyond 10 years and most die at age 1-2 years, but are considered middle-aged in captivity where mortality increases well beyond 25 years [8-10]. Second, it is important to note that while MA NMR skin structure was found to be similar to that of young animals, the thickness of the epidermis appeared paradoxically increased at an older age, as a consequence of increased number of epidermal layers. Marked thickness of the corneous layer in the adult NMR has previously been described [18] and lack of fur has also been described to be compensated by a thicker epidermis [17]. The reason for such an increase remains unknown. It could result from a reaction to the higher cumulative skin friction in their tunnel habitats for older animals. It could also be due to the overexpression of *Igfbp3* in the basal/ germinative layer of the epidermis in older animals, leading to a modulation of the early stages of epidermis differentiation. Third, it should also be noted that we focused our study on skin epidermal keratinocytes although it has been demonstrated that human dermal fibroblasts undergo a partial loss of their cellular identity during aging [27]. Finally, alterations in hair follicle stem cell niche -that occur in mice aging- was not possible to address in this hairless rodent [38].

To conclude, NMR live more than 37 years, while they should only live for 6 years, based on their size. Therefore, the maintenance of their epidermal stem cell and LC compartments -possibly in relation with increased *Igfbp3* expression and markers of immunity- appear to contribute to their remarkable longevity and resistance to skin aging. Our data also confirm the relevance of the NMR as a model organism for studying cutaneous biology and disease resistance and provide

valuable clues to better understand the cutaneous biology of this animal.

# MATERIALS AND METHODS

## Animals

Young (Y, mean age: 1.1 years) and middle-aged (MA, mean age: 11.3 years) males were obtained from four distinct colonies maintained at the Ecole Nationale Vétérinaire de Maisons-Alfort (EnVA). Animal experiments were performed according to experimental protocols following European Community Council guidelines and approved by the Institutional Animal Care and Use Committee of the Ecole Nationale Vétérinaire de Maisons-Alfort (n°016).

## Histology and immunohistochemistry

MA and Y NMR skin samples were collected under anesthesia by performing 5mm skin punch biopsies on their back skin. The skin was then immerged in formalin for 12 to 24 hours at 4° C, and embedded either in paraffin or OCT compound for sectioning.

Formalin-fixed paraffin-embedded skin sections (5 µm) were stained with Hematoxylin and Eosin, Sirius Red and Masson's Trichrome using standard procedures. Formalin-fixed paraffin-embedded skin sections or cryo-sections were also labeled with various antibodies (Supplementary Table 1). Briefly, for formalin-fixed paraffin-embedded sections, a prior antigen retrieval was performed by incubating the slides with citrate buffer (pH=6) or EDTA buffer (pH=9) at 98° C for 20 min. Sections were then washed with Tris buffer saline (TBS) (137mM NaCl and 20mM Tris) / Tween 20 (TBST) 0.1% and blocked for 2 hours with TBST/Normal Goat Serum (NGS) 5%. The sections were incubated with the different primary antibodies diluted in blocking buffer overnight at 4° C for intracellular staining, followed by secondary antibody goat anti-mouse Alexa 546 (dilution 1:1000, Life Technologies), or goat anti-rabbit Alexa 488 (dilution 1:1000, Life Technologies) for 45 min at room temperature. Slides were also incubated with 4',6-diamidino-2-phenylindole (DAPI, Sigma-Aldrich) for 5 min for nuclei staining and mounted using Fluoromount G (Southern Biotech, Birmingham, AL, USA). For

hyaluronic acid staining, the slides were incubated 2h with HA-binding protein (at 1,6/100 dilution) with 0.2% bovine serum albumin and 0.02% TritonX-100 in PBS. HA-binding protein staining was achieved with 1/100 strepavidin HRP (Perkin Elmer, Waltham MA, USA) solution for 30 min at room temperature and 3,3'-diaminobenzidine (Liquid DAB + substrate Chromagen
System; Dako, Santa Clara, CA, USA) for color
development. Finally, sections were counterstained with
hematoxylin (Sigma-Aldrich, St. Louis, MO, USA) and
mounted using DPX (HA staining in brown). Negative
controls were performed in parallel with the samples
substituting the primary antibody with the equivalent
isotype.

The sections were visualized under an Olympus BX63F microscope equipped with an Olympus DP73 camera (and pictured with Metamorph software; Olympus, Tokyo, Japan). Using Image J software, epidermal, dermal thickness and number of epidermal layers and cells were measured. Epidermal cells or layers were counted in 10-3mm2 randomized surface squares in the epidermis of each specimen. Quantification was done on various parts of the skin, dissected in a gridded fashion and 3 different sectors were analyzed; for each different sector of the skin, 3 different areas of 10-3mm2 were counted. Similarly, immunolabeled cells were counted and fluorescence intensity was measured with mean gray value per field. For IGFBP3 staining, cell nuclei were automatically segmented and counted using a deep learning-based segmentation algorithm called Cellpose (https://www.cellpose.org/).

Single-cell RNA-sequencing

Epidermal cells from 3 Y and 3 MA NMR were used. All the samples were collected and processed at the same time. NMR were anesthetized by the inhalation of 4.9% isoflurane delivered at a flow rate of 300 ml min⁻¹ in ambient air and ethanol was applied topically to the dorsal skin for 20 seconds. Punch biopsy devices were used to create 5 mm surgical wounds on the backs of the animals. To preserve cell viability, skin pieces were immediately placed in PBS/BSA 3% medium on ice until dissection. The wounds were rinsed with betadine solution. Using a scalpel, the subcutaneous fat was removed from the dermis. The cleaned tissues were placed dermis side down in a sterile petri dish in a trypsin 0.25%/EDTA 1mM solution at 37° C for 1h. The enzymatic digestion was stopped by addition of PBS/FBS 10%. This floating procedure allows the digestion and dermis-epidermis separation by lifting the epidermis up straight above the dermis. Epidermis was then mechanically dissociated and minced in small pieces with razor blades in PBS and then triturated in 25mL, 10mL then 5mL pipets. 3 rounds of washes were performed starting with 100µm filtration, 70µm and 40µm and interspersed with centrifugation cycles at 150g at 4° C for 5 minutes. Finally, the cell pellet was resuspended in PBS/BSA 3% and cell viability was measured using trypan's blue method on Kova slides. All single cell suspensions had a viability greater than 75%. Single cells then were captured in droplet

emulsion using the Chromium Controller (10x
Genomics, Pleasanton, CA, USA), and scRNA-seq
libraries were constructed according to the 10x
Genomics protocol using the Chromium Single-Cell 3'
Gel Bead and Library V2 Kit (10x Genomics,
Pleasanton, CA, USA). In brief, cell suspensions were
diluted in PBS/FBS 3% to a final concentration of 1 x
10⁶ cells/mL (1,000 cells per µL). Cells were loaded in
each channel with a target output of 5,000 cells per
sample. All reactions were performed in a C1000 Touch
Thermal Cycler (Bio-Rad Laboratories, Hercules, CA,
USA) with a 96 Deep Well Reaction Module. Twelve
cycles were used for cDNA amplification and sample
index PCR. Amplified cDNA and final libraries were
evaluated using a Bioanalyzer 2100 (Agilent
Technologies, Santa Clara, CA, USA) with a high
sensitivity chip. Samples were sequenced on an HiSeq
4000 instrument (Illumina, San Diego, CA, USA).

scRNA-seq pre-processing

scRNA-seq data analyses were performed by
GenoSplice technology (http://www.genosplice.com).
Sequencing data quality analysis was performed using
FastQC v0.11.2 on 6 NMR epidermal cells samples, 3
from a Y NMR group (samples 1J, 5J, 6J) and 3 from an
MA NMR group (samples 2V, 3V, 4V). For read
alignment and unique molecular identifiers (UMI)
quantification, CellRanger software v3.0.2 was used on
Heterocephalus glaber genome HetGla_female_1.0 with
default parameters and gene annotation from Ensembl
96. The 6 expression matrices containing the UMI
counts were merged, and only the genes with UMI ≥ 1
in at least one cell were kept. The following filters were
applied to generate a global matrix used in further
analysis: cells with UMI ≥ 2000, number of detected
genes ≥ 700, and cells with UMI in mitochondrial genes
≤ 25%. For UMIs normalization, Seurat v3.1.1 was used
[PMID:29608179], and global-scaling normalization
method was applied with a scale factor of 10,000 and
log-transformation of data. This was followed by a
scaling linear transformation step, to avoid highly
expressed genes having higher weight in downstream
analysis.

Clustering and marker genes

PCA was performed on the scaled data, with a
Jackstraw plot to choose how many PCs to retain as an
input for Seurat clustering step. All samples were
evenly distributed among clusters, except for sample 2V
that showed a very distinctive pattern and was grouped
in separated clusters. This sample was therefore
removed from the expression matrix in further analysis,
keeping 3 samples in the Y NMR group (1J, 5J, 6J) and
two samples in the MA NMR group (3V, 4V). A second
clustering step was performed with the five remaining samples using default parameters, Louvain algorithm as the clustering method and a resolution parameter defining the clusters granularity set to 0.4. Marker genes defining each cluster were found via differential expression testing, with a Wilcoxon rank sum test and a log fold change threshold (logFc) of 1. The found marker genes were first compared to a list of 98 genes of interest for cluster identification, extracted from previous single-cell transcriptomic skin data [21, 22]. Then, we also compared the marker genes to published scRNA-seq transcript or protein expression data documenting skin layer specificity (Supplementary Data File 1). Initial 11 clusters were then defined separately based on a specific set of genes and next grouped in the following metaclusters: basal 1, basal 2, intermediate 1, intermediate 2, corneum and immune cells.

## Trajectory

Single-cell trajectory was constructed using M3Drop [PMID:30590489] and Monocle v2.10.1 [PMID:28825705]. Input genes for Monocle trajectory construction were selected using an unsupervised approach via M3Drop result, which identifies differentially expressed genes based on a Michaelis-Menten function for the relationship between mean expression and dropout rate, the relevant genes being the ones shifting above a fitted curve. The default Monocle workflow was then performed to generate the trajectories, with focus on known markers expression profiles (*Krt14*, *Krt10*, *Lor*) to select the appropriate trajectory orientation. We utilized Monocle 2.6.425 to order cells in pseudotime based on their transcriptomic similarity. Pseudotime-dependent genes have been identified using differential GeneTest function from Monocle. Those genes were the most significant in this context and Monocle reconsidered this list of 820 genes and kept the majority of them for the analyses (791 pseudotime-dependent genes). The differential test was performed on genes used to construct Monocle trajectory on cells from Young samples and selected according to q-value (lower than 0.05). Clustering of these genes was done by Monocle functions. Heatmaps for Y and MA samples used this clustering of genes.

## Enrichment analysis

Using the Cytoscape v3.7.2 software and BiNGO plugin (binomial statistical test and Bonferroni correction, threshold with adjusted p<0.01), we determined which Gene Ontology (GO) categories are statistically overrepresented in a set of genes of a biological network in the newly defined metaclusters on Mus Musculus annotations. For specific analysis, a Fisher exact-test was performed on a list of genes compiled

from REACTOME (ROS) and GO terms (aging and wound healing) databases. Results with uncorrected p-value ≤ 0.05 were considered enriched.

## Wound healing procedure

NMR were anesthetized and punch biopsy devices were used to create four 5 mm surgical wounds on the backs of the animals. All tissues above the panniculus carnosus were excised. Wounds were left uncovered until analysis. Wound tissues were harvested at day 7 post injury and stored at -80° C after formalin fixation and sucrose washing/immersion. Standardized images of the wounds were obtained at various time points, with a Sony Cybershot 10.1-megapixel DSC-W180 digital camera (Sony, Tokyo, Japan). The images were used to measure the wound closure rate using ImageJ software every 2-3 days.

## Statistics

The data were analyzed using Chi² statistic test, Wilcoxon test, Unpaired Student's t-test, Mann-Whitney U test and One-way ANOVA. n indicates the number of independent experiments performed. In all histograms, asterisks correspond to: \*p < 0.05, \*\*p < 0.01, \*\*\*p < 0.001, \*\*\*\*p < 0.0001.

## Abbreviations

NMR: Naked mole-rats; scRNA-seq: RNA-sequencing; Krt: Keratin; Lor: Loricrin; Y: Young; MA: Middle-aged; UMAP: Uniform Manifold Approximation and Projection; t-SNE: t-distributed Stochastic Neighbor Embedding; BSC: basal and stem cells; BPSC: basal proliferating and stem cells; Sp: spinous layer cells; Gr: granular layer cells; Cor: corneous layer cells.

## AUTHOR CONTRIBUTIONS

Study conception and design: M.V, G.F, S.A, R.H.F. Acquisition of data: A.S, T.J, P.M, C.L, D.D.M, L.L, R.H.F. Analysis and interpretation of data: A.S, T.J, F.S, C.G.F, P.M, M.V, G.F, S.A, R.H.F. Drafting and/or critical revision of manuscript: A.S, T.J, F.S, C.G.F, P.C, M.V, G.F, S.A, R.H.F. All authors have read and approved the last version of the manuscript.

## ACKNOWLEDGMENTS

We thank Genosplice Technology (http://www.genosplice.com) for scRNA-seq data analyses and Daniel Vaiman, Dany Nassar, Bénédicte Oulès, Matthieu Tihy, Maria Sbeih and Mathieu Castela for helpful discussions. We also thank Aude Fert, Michèle Oster, Sandrine Lamblin, Romain Morichon
(Imaging Platform UMS30 LUMIC), IMAG'IC and HISTIM platforms of Institut Cochin, and Thomas Lilin (Ecole Nationale Vétérinaire d'Alfort) for their technical support and Tristan Denis for graphical design. M.V. is grateful to the Ecole Nationale Vétérinaire d'Alfort (EnvA) and the Association De Prévoyance Santé (ADPS) partners of Allianz, for their invaluable help in establishing and funding the animal facility.

## CONFLICTS OF INTEREST

The authors declare that they have no conflicts of interest.

## FUNDING

This study was supported by the IdEx Université de Paris, ANR-18-IDEX-0001, Société Française de Dermatologie and Fondation pour la Recherche en Physiologie.

## REFERENCES

1. Eyerich S, Eyerich K, Traidl-Hoffmann C, Biedermann T. Cutaneous Barriers and Skin Immunity: Differentiating A Connected Network. Trends Immunol. 2018; 39:315-27.
https://doi.org/10.1016/j.it.2018.02.004
PMID:29551468

2. Blanpain C. Stem cells: Skin regeneration and repair. Nature. 2010; 464:686-7.
https://doi.org/10.1038/464686a
PMID:20360726

3. Goodell MA, Rando TA. Stem cells and healthy aging. Science. 2015; 350:1199-204.
https://doi.org/10.1126/science.aab3388
PMID:26785478

4. Giangreco A, Qin M, Pintar JE, Watt FM. Epidermal stem cells are retained in vivo throughout skin aging. Aging Cell. 2008; 7:250-9.
https://doi.org/10.1111/j.1474-9726.2008.00372.x
PMID:18221414

5. Ashcroft GS, Mills SJ, Ashworth JJ. Ageing and wound healing. Biogerontology. 2002; 3:337-45.
https://doi.org/10.1023/a:1021399228395
PMID:12510172

6. Xu YP, Qi RQ, Chen W, Shi Y, Cui ZZ, Gao XH, Chen HD, Zhou L, Mi QS. Aging affects epidermal Langerhans cell development and function and alters their miRNA gene expression profile. Aging (Albany NY). 2012; 4:742-54.
https://doi.org/10.18632/aging.100501
PMID:23178507

7. Eming SA, Martin P, Tomic-Canic M. Wound repair and regeneration: mechanisms, signaling, and translation.

Sci Transl Med. 2014; 6:265sr6.
https://doi.org/10.1126/scitranslmed.3009337
PMID:25473038

8. Lee BP, Smith M, Buffenstein R, Harries LW. Negligible senescence in naked mole rats may be a consequence of well-maintained splicing regulation. Geroscience. 2020; 42:633-51.
https://doi.org/10.1007/s11357-019-00150-7
PMID:31927681

9. Lowe R, Danson AF, Rakyan VK, Yildizoglu S, Saldmann F, Viltard M, Friedlander G, Faulkes CG. DNA methylation clocks as a predictor for ageing and age estimation in naked mole-rats, Heterocephalus glaber. Aging (Albany NY). 2020; 12:4394-406.
https://doi.org/10.18632/aging.102892
PMID:32126024

10. Hochberg ME, Noble RJ, Braude S. A Hypothesis to Explain Cancers in Confined Colonies of Naked Mole Rats. bioRxiv. 2016.
https://doi.org/10.1101/079012

11. Buffenstein R, Ruby JG. Opportunities for new insight into aging from the naked mole-rat and other non-traditional models. Nature Aging. 2021; 1: 3-4.
https://doi.org/10.1038/s43587-020-00012-4

12. Ruby JG, Smith M, Buffenstein R. Naked Mole-Rat mortality rates defy gompertzian laws by not increasing with age. Elife. 2018; 7:e31157.
https://doi.org/10.7554/eLife.31157
PMID:29364116

13. Lagunas-Rangel FA, Chávez-Valencia V. Learning of nature: The curious case of the naked mole rat. Mech Ageing Dev. 2017; 164:76-81.
https://doi.org/10.1016/j.mad.2017.04.010
PMID:28472634

14. Edrey YH, Park TJ, Kang H, Biney A, Buffenstein R. Endocrine function and neurobiology of the longest-living rodent, the naked mole-rat. Exp Gerontol. 2011; 46:116-23.
https://doi.org/10.1016/j.exger.2010.09.005
PMID:20888895

15. Liang S, Mele J, Wu Y, Buffenstein R, Hornsby PJ. Resistance to experimental tumorigenesis in cells of a long-lived mammal, the naked mole-rat (Heterocephalus glaber). Aging Cell. 2010; 9:626-35.
https://doi.org/10.1111/j.1474-9726.2010.00588.x
PMID:20550519

16. Kim EB, Fang X, Fushan AA, Huang Z, Lobanov AV, Han L, Marino SM, Sun X, Turanov AA, Yang P, Yim SH, Zhao X, Kasaikina MV, et al. Genome sequencing reveals insights into physiology and longevity of the naked mole rat. Nature. 2011; 479:223-7.
https://doi.org/10.1038/nature10533
in normal and diseased human skin. J Cutan Pathol.
2000; 27:108-11.

https://doi.org/10.1034/j.1600-0560.2000.027003108.x

PMID: 10728810

35. Moriyama M, Durham AD, Moriyama H, Hasegawa K, Nishikawa S, Radtke F, Osawa M. Multiple roles of Notch signaling in the regulation of epidermal development. Dev Cell. 2008; 14:594-604.

https://doi.org/10.1016/j.devcel.2008.01.017

PMID: 18410734

36. Haensel D, Jin S, Sun P, Cinco R, Dragan M, Nguyen Q, Cang Z, Gong Y, Vu R, MacLean AL, Kessenbrock K, Gratton E, Nie Q, Dai X. Defining Epidermal Basal Cell States during Skin Homeostasis and Wound Healing Using Single-Cell Transcriptomics. Cell Rep. 2020; 30:3932-47.e6.

https://doi.org/10.1016/j.celrep.2020.02.091

PMID: 32187560

37. Quek LS, Grasset N, Jasmen JB, Robinson KS, Bellanger S. Dual Role of the Anaphase Promoting Complex/Cyclosome in Regulating Stemness and Differentiation in Human Primary Keratinocytes. J Invest Dermatol. 2018; 138:1851-61.

https://doi.org/10.1016/j.jid.2018.02.033

PMID: 29526760

38. Ge Y, Miao Y, Gur-Cohen S, Gomez N, Yang H, Nikolova M, Polak L, Hu Y, Verma A, Elemento O, Krueger JG, Fuchs E. The aging skin microenvironment dictates stem cell behavior. Proc Natl Acad Sci USA. 2020; 117:5339-50.

https://doi.org/10.1073/pnas.1901720117

PMID: 32094197

39. Rubin AJ, Barajas BC, Furlan-Magaril M, Lopez-Pajares V, Mumbach MR, Howard I, Kim DS, Boxer LD, Cairns J, Spivakov M, Wingett SW, Shi M, Zhao Z, et al. Lineage-specific dynamic and pre-established enhancer-promoter contacts cooperate in terminal differentiation. Nat Genet. 2017; 49:1522-8.

https://doi.org/10.1038/ng.3935

PMID: 28805829

40. Favre B, Plantard L, Aeschbach L, Brakch N, Christenzaech S, de Viragh PA, Sergeant A, Huber M, Hohl D. SLURP1 is a late marker of epidermal differentiation and is absent in Mal de Meleda. J Invest Dermatol. 2007; 127:301-8.

https://doi.org/10.1038/sj.jid.5700551

PMID: 17008884

41. Srivastava A, Barth E, Ermolaeva MA, Guenther M, Frahm C, Marz M, Witte OW. Tissue-specific Gene Expression Changes Are Associated with Aging in Mice. Genomics Proteomics Bioinformatics. 2020; 18:430-42.

https://doi.org/10.1016/j.gpb.2020.12.001

PMID: 33309863

42. Fuchs E. Epidermal differentiation. Curr Opin Cell Biol.
1990; 2:1028-35.

https://doi.org/10.1016/0955-0674(90)90152-5

PMID: 1712211

43. Fang X, Seim I, Huang Z, Gerashchenko MV, Xiong Z, Turanov AA, Zhu Y, Lobanov AV, Fan D, Yim SH, Yao X, Ma S, Yang L, et al. Adaptations to a subterranean environment and longevity revealed by the analysis of mole rat genomes. Cell Rep. 2014; 8:1354-64.

https://doi.org/10.1016/j.celrep.2014.07.030

PMID: 25176646

44. Yang R, Liu F, Wang J, Chen X, Xie J, Xiong K. Epidermal stem cells in wound healing and their clinical applications. Stem Cell Res Ther. 2019; 10:229.

https://doi.org/10.1186/s13287-019-1312-z

PMID: 31358069

45. Cohen BJ, Danon D, Roth GS. Wound repair in mice as influenced by age and antimacrophage serum. J Gerontol. 1987; 42:295-301.

https://doi.org/10.1093/geronj/42.3.295

PMID: 3571865

46. Groves RW, Liu L, Dopping-Hepenstal PJ, Markus HS, Lovell PA, Ozoemena L, Lai-Cheong JE, Gawler J, Owaribe K, Hashimoto T, Mellerio JE, Mee JB, McGrath JA. A homozygous nonsense mutation within the dystonin gene coding for the coiled-coil domain of the epithelial isoform of BPAG1 underlies a new subtype of autosomal recessive epidermolysis bullosa simplex. J Invest Dermatol. 2010; 130:1551-7.

https://doi.org/10.1038/jid.2010.19

PMID: 20164846

47. Favre B, Begré N, Borradori L. A recessive mutation in the DSP gene linked to cardiomyopathy, skin fragility and hair defects impairs the binding of desmoplakin to epidermal keratins and the muscle-specific intermediate filament desmin. Br J Dermatol. 2018; 179:797-9.

https://doi.org/10.1111/bjd.16832 PMID: 29878302

48. Tan DW, Jensen KB, Trotter MW, Connelly JT, Broad S, Watt FM. Single-cell gene expression profiling reveals functional heterogeneity of undifferentiated human epidermal cells. Development. 2013; 140:1433-44.

https://doi.org/10.1242/dev.087551

PMID: 23482486

49. Egger AN, Rajabiestarabadi A, Williams NM, Resnik SR, Fox JD, Wong LL, Jozic I. The importance of caveolins and caveolae to dermatology: Lessons from the caves and beyond. Exp Dermatol. 2020; 29:136-48.

https://doi.org/10.1111/exd.14068 PMID: 31845391

50. Tabula Muris Consortium. A single-cell transcriptomic atlas characterizes ageing tissues in the mouse. Nature. 2020; 583:590-5.
https://doi.org/10.1038/s41586-020-2496-1
PMID:32669714

51. Zou Z, Long X, Zhao Q, Zheng Y, Song M, Ma S, Jing Y, Wang S, He Y, Esteban CR, Yu N, Huang J, Chan P, et al. A Single-Cell Transcriptomic Atlas of Human Skin Aging. Dev Cell. 2021; 56:383-97.e8.
https://doi.org/10.1016/j.devcel.2020.11.002
PMID:33238152

52. Hilton HG, Rubinstein ND, Janki P, Ireland AT, Bernstein N, Fong NL, Wright KM, Smith M, Finkle D, Martin-McNulty B, Roy M, Imai DM, Jojic V, Buffenstein R. Single-cell transcriptomics of the naked mole-rat reveals unexpected features of mammalian immunity. PLoS Biol. 2019; 17:e3000528.
https://doi.org/10.1371/journal.pbio.3000528
PMID:31751331

53. Emmrich S, Tolibzoda Zakusilo F, Trapp A, Zhou X, Zhang Q, Irving EM, Drage MG, Zhang Z, Gladyshev VN, Seluanov A, Gorbunova V. Ectopic cervical thymi and no thymic involution until midlife in naked mole rats. Aging Cell. 2021; 20:e13477.
https://doi.org/10.1111/acel.13477
PMID:34596321

54. Kawamura Y, Oka K, Takamori M, Sugiura Y, Oiwa Y, Fujioka S, Homma S, Miyawaki S, Narita M, Fukuda T, Suematsu M, Bono H, Okano H, et al. Senescent cell death as an aging resistance mechanism in naked mole-rat. bioRxiv. 2020.
https://doi.org/10.1101/2020.07.02.155903
SUPPLEMENTARY MATERIALS

Supplementary Figures

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

**Supplementary Figure 1.** Epidermal morphology of young, middle-aged and aged naked mole-rat (NMR). (A) HE staining with a higher magnification of an aged (25 years-old) NMR skin. Scale bar = 100 µm. (B) Histogram shows the increased epidermal thickness of aged and middle-aged as compared to young NMR back skin. (C) Examples of two aged (20 years-old) NMR skin histology using HES staining. Scale bar = 50 µm. (D) Histogram showing the number of epidermal cells in middle-aged versus young NMR per surface. (E) Histogram showing the number of epidermal layers in middle-aged vs young NMR. (F) Histogram showing cell number per layer in middle-aged versus young NMR, per surface. An additional external layer was found in middle-aged samples. n=4 animals per group except for the aged animals (n=3). Bars: SEM. \*represents differences between the age groups. \*\*\*\*p < 0.0001 in ANOVA and Student's t-test. ND: non detectable. ns = non significant.
<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /><figcaption>**Supplementary Figure 2.** **Quality controls of scRNA-seq.** (A) Violin plot of the number of UMs per cell for each animal. A minimum threshold is applied at 2,000 (red line). (B) Violin plot of the number of genes detected per cell for each animal. A minimum threshold is applied at 700 (red line). (C) Violin plot of the percentage of mitochondrial reads per cell for each animal. A minimum threshold is applied at 25% (red line). (D) Bar plot of the number of cells for each animal. (E) Uniform Manifold Approximation and Projection (UMAP) representation for all samples grouped and individualized, colored by individual, including the outlier (2V) removed in the rest of the study.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>**Supplementary Figure 3.** scRNA-seq of NMR epidermal cells define 3 cellular states. (A) Violin plots of *Krt14*, *Krt10* and *Lor* marker genes expressed by NMR epidermal cells, leading us to confirm 3 cellular states of keratinization. (B) t-distributed Stochastic Neighbor Embedding (t-SNE) visualization of epidermal cells clustering. (C) t-SNE visualization of *Krt14*, *Krt10* and *Lor* marker genes expressed by NMR epidermal cells. BSC = basal and stem cells; BPSC = basal proliferating and stem cells; Sp = Spinous layer cells; Gr = granular layer cells; Cor = corneous layer cells.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /><figcaption>**Supplementary Figure 4. scRNA-seq of young versus middle-aged NMR epidermis.** (A) t-distributed Stochastic Neighbor Embedding (t-SNE) visualization of epidermal cells where both middle-aged (MA, red) and young (Y, blue) NMR were jointly projected on the same plot, showing an overlap of the 2 age groups tested. (B) Unsupervised differentiation trajectories for merged, MA and Y keratinocytes colored by Monocle state of differentiation. A single and similar trajectory in both samples was found. (C) Unsupervised differentiation trajectories for merged, MA and Y keratinocytes segregated by cluster. No difference between the 2 age groups was found. BSC = basal and stem cells; BPSC = basal proliferating and stem cells; Sp = Spinous layer cells; Gr = granular layer cells; Cor = corneous layer cells.</figcaption></figure>
<figure><img src="image_66.png" /><figcaption>A</figcaption></figure>

Krt14

Krt10

Lor

<figure><img src="image_67.png" /></figure>

Identity

<figure><img src="image_68.png" /><figcaption>B</figcaption></figure>

Metacluster

<figure><img src="image_69.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_70.png" /></figure>

**Supplementary Figure 5. Metaclusters identification and pseudotime estimation and of young and middle-aged NMR epidermal cells.** (A) The expression of *Krt14*, *Krt10* and *Lor* marker genes visualized with violin plots confirmed the cellular states of the newly defined metaclusters for middle-aged (MA, red) and young (Y, blue) NMR. (B) Bar graph representing the relative proportion of epidermal cells in each metacluster between MA (red) and Y (blue) animals. No significant difference between the 2 age groups was noted (Chi2 statistic test). (C) Distributions plotting the frequency of MA and Y cells as a function of the pseudotime on its scale. Distributions were similar between the 2 age groups. (D) Heatmap of most differentially expressed genes for each newly defined metacluster. Cells are represented in columns, and genes in rows.
<figure><img src="image_71.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_74.png" /></figure>

**Supplementary Figure 6. Metaclusters characterization of NMR keratinocytes.** (A) Heatmap of pseudotemporal dynamics of the pseudotime-dependent genes in MA (middle-aged) and Y (young) NMR (left panel). Each row (i.e., gene) is normalized to its peak value along the pseudotime. 6 distinct clusters or stages during pseudotime are represented by colored bars on the side. The most significant enriched GO terms in each gene cluster are listed (right panel). The gene ordering is strictly the same between MA and Y heatmaps. (B) Barplot representing a selection of the most significant enriched Gene Ontology (GO) terms (from BiNGO analysis) in each keratinocytes metaclusters (except immune cluster) sorted by p-value. Barplot represents the adjusted p-value. (C) Venn diagram depicting the genes overlap significantly overrepresented in basal1 (pval 5,86e-03) and corneum metaclusters (pval 1,01e-02) as compared respectively with wound healing and ROS gene lists. (D) Bubble plot representing the immune Gene Ontology (GO) terms released by the BiNGO analysis. The bubbles represent the adjusted p-value of these terms in the immune cluster and in the other clusters when they emerged.
<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /><figcaption>**Supplementary Figure 7.** Wound healing of young versus middle-aged NMR. (A) Photomicrographs of back skin wounds for middle-aged and young NMR at day 0 and day 7 post wounding. Skin wounds were harvested at day 7 to analyze re-epithelialization. (B) Histogram showing the neoepidermal thickness (μm) of middle-aged and young NMR wounded back skin at day 7. Bars: SEM. A Mann-Whitney U test was used. n=2 animals per group.</figcaption></figure>
Supplementary Table

<table><caption>Supplementary Table 1. Characteristics of antibodies used for immunohistochemistry.</caption><thead><tr><th>Reagent or resource</th><th>Dilution</th><th>Species</th><th>Source</th><th>Identifier</th><th>Section method</th></tr></thead><tbody><tr><td colspan="6"><strong>Antibodies</strong></td></tr><tr><td>Keratin-14</td><td>1:1000</td><td>rabbit</td><td>Biolegend</td><td>905301</td><td>paraffin</td></tr><tr><td>Keratin-10</td><td>1:1000</td><td>rabbit</td><td>Biolegend</td><td>905401</td><td>cryo</td></tr><tr><td>Loricrin</td><td>1:1000</td><td>rabbit</td><td>Biolegend</td><td>905101</td><td>paraffin</td></tr><tr><td>Filaggrin</td><td>1:500</td><td>rabbit</td><td>Biolegend</td><td>905801</td><td>paraffin</td></tr><tr><td>Laminin-5</td><td>1:200</td><td>rabbit</td><td>Abcam</td><td>ab14509</td><td>cryo</td></tr><tr><td>ITGB4</td><td>1:250</td><td>rabbit</td><td>Abcam</td><td>ab182120</td><td>paraffin</td></tr><tr><td>Ki67</td><td>1:400</td><td>rabbit</td><td>Abcam</td><td>ab15580</td><td>cryo</td></tr><tr><td>Collagen I</td><td>1:250</td><td>rabbit</td><td>Abcam</td><td>ab21286</td><td>paraffin</td></tr><tr><td>p16</td><td>1:100</td><td>rabbit</td><td>ThermoScientific</td><td>PA1-30670</td><td>cryo</td></tr><tr><td>Hyaluronic acid</td><td>1,6:100</td><td>bovine</td><td>Calbiochem</td><td>385911</td><td>paraffin</td></tr><tr><td>IGFBP3</td><td>1:500</td><td>rabbit</td><td>Abcam</td><td>272857</td><td>cryo</td></tr></tbody></table>
Supplementary Data Files

Please browse Full Text version to see the data of Supplementary Data Files 1 to 4.

**Supplementary Data File 1.** (A) List of 98 characteristic genes of epidermal populations. (B) Top 25 marker genes for each cluster of epidermal cells from young and middle-aged naked mole-rats. (C) Differentially expressed gene signatures for each cluster of epidermal cells from young and middle-aged naked mole-rats.

**Supplementary Data File 2.** (A) Differential expression analysis comparing metacluster basal1 subpopulations of naked mole-rat keratinocytes with all the other metaclusters : markers genes of metacluster basal 1. (B) BiNGO Enrichment analysis using markers genes of metacluster basal 1 of naked mole-rat keratinocytes with all the other metaclusters. (C) Differential expression analysis comparing metacluster basal2 subpopulations of naked mole-rat keratinocytes with all the other metaclusters : markers genes of metacluster basal 2. (D) BiNGO Enrichment analysis using markers genes of metacluster basal 2 of naked mole-rat keratinocytes with all the other metaclusters. (E) Differential expression analysis comparing metacluster intermediate1 subpopulations of naked mole-rat keratinocytes with all the other metaclusters : markers genes of metacluster intermediate 1. (F) BiNGO Enrichment analysis using markers genes of metacluster intermediate1 of naked mole-rat keratinocytes with all the other metaclusters. (G) Differential expression analysis comparing metacluster intermediate2 subpopulations of naked mole-rat keratinocytes with all the other metaclusters : markers genes of metacluster intermediate 2. (H) GO Enrichment analysis using markers genes of metacluster intermediate2 of naked mole-rat keratinocytes with all the other metaclusters. (I) Differential expression analysis comparing metacluster corneum subpopulations of naked mole-rat keratinocytes with all the other metaclusters : markers genes of metacluster corneum. (J) GO Enrichment analysis using markers genes of metacluster corneum of naked mole-rat keratinocytes with all the other metaclusters. (K) Differential expression analysis comparing metacluster immune subpopulations of naked mole-rat keratinocytes with all the other metaclusters : markers genes of metacluster immune. (L) GO Enrichment analysis using markers genes of metacluster immune of naked mole-rat keratinocytes with all the other metaclusters.

**Supplementary Data File 3.** (A) Specific enrichment analysis of naked mole-rats keratinocytes for ROS, aging and wound healing functions. (B) Databases used for specific enrichment analysis of naked mole rats keratinocytes for ROS, aging and wound healing functions.

**Supplementary Data File 4.** (A) Differentially expressed genes by cluster between young and middle-aged animals. (B) Differentially expressed genes by metacluster between young and middle-aged animals.