Decoding the Human Epidermal Complexity at Single-Cell Resolution

The epidermis is one of the largest tissues in the human body, serving as a protective barrier. The basal layer of the epidermis, which consists of epithelial stem cells and transient amplifying progenitors, represents its proliferative compartment. As keratinocytes migrate from the basal layer to the skin surface, they exit the cell cycle and initiate terminal differentiation, ultimately generating the suprabasal epidermal layers. A deeper understanding of the molecular mechanisms and pathways driving keratinocytes’ organization and regeneration is essential for successful therapeutic approaches. Single-cell techniques are valuable tools for studying molecular heterogeneity. The high-resolution characterization obtained with these technologies has identified disease-specific drivers and new therapeutic targets, further promoting the advancement of personalized therapies. This review summarizes the latest findings on the transcriptomic and epigenetic profiling of human epidermal cells, analyzed from human biopsy or after in vitro cultivation, focusing on physiological, wound healing, and inflammatory skin conditions.


Introduction
The human epidermis, the outermost layer of the skin, is composed of stratified squamous epithelium that acts as a protective barrier against external insults; hence, any alteration can drastically impair patients' life quality [1][2][3].
The basal layer represents the proliferative compartment of the epidermis, in which both epithelial stem cells (SCs) and transient amplifying progenitors (TACs) reside. As keratinocytes leave the basal layer and move towards the skin surface, they exit the cell cycle and undergo a terminal differentiation program, generating the suprabasal epidermal layers (spinosum, granulosum, lucidum, and corneum) [4,5]. This process, known as cornification or keratinization, consists of morphologic and metabolic changes whose endpoint is to form the cornified layer or the stratum corneum, the outermost layer of the epidermis [6].
In the last four decades, the ex vivo restoration of functional epithelia allowed the treatment of different medical needs [7][8][9][10]. An astonishing example has been the lifesaving treatment of third-degree burn and junctional epidermolysis bullosa (JEB)-affected patients [10] thanks to the knowledge acquired during the studies conducted by H. Green in the 1980s [7]. The successful outcome of therapeutic approaches relies on a fine understanding of the molecular mechanisms and pathways driving human keratinocytes' organization and regeneration. In this context, single-cell technologies have already and will further increase our knowledge on in vivo and in vitro cellular complexity.
Single-cell RNA sequencing (scRNA-seq), unlike bulk analysis, allows researchers to decipher the complexity of biological systems at the single-cell level. In 2006, a pioneering Given the considerable skin-extension and site-specific functions, experimental settings (in vivo or in vitro) ( Table 1), biopsy withdrawal and processing methods (Table 1), the anatomical areas of collection (Figure 1), and the age and pathology of the donors ( Figure 2) must be considered to properly analyze single cells' published data.
The less invasive and more common biopsy techniques for research purposes are shaved and punch biopsies, which enable the collection of samples from both the epidermis and the underlying dermis. Depending on the cell populations of interest, a skin biopsy can be processed in different ways: heat, chemical reagents, enzymes, or mechanical digestion. However, heat treatment may cause thermal damage, and chemical reagents may alter cellular electrolyte equilibrium. Hence, the most used methods rely on enzymatic activity or mechanical separation to divide the epidermis from the dermis [21]. Several proteases (dispases, trypsins, pancreatin, pronase, and thermolysin) have been tested to separate different layers of the dermal-epidermal junction [21,22].
A clinically validated method specifically used for the in vitro culture of human primary keratinocytes consists of a biopsy cleaning step to remove the adipose tissue and partially the dermis. Then, the tissue is minced and incubated in a trypsinizing flask with a trypsin/EDTA solution (mixture of trypsin, chymotrypsin, and elastase) for 30 min at 37 • C. Cells are recovered after each round of trypsinization and plated onto a lethally irradiated feeder layer [23]. This procedure allows the collection of both interfollicular and follicular keratinocytes, including SCs used for cell and gene therapy applications (unpublished data and [10,24,25]).
Another popular method takes advantage of dispase I, one of the most used enzymes to gently separate the dermis from the epidermal layer. It acts on extracellular matrix (ECM) proteins, including fibronectin, collagen IV, and to a lesser extent collagen I, allowing the collection of the entire epidermis [22,26]. Then, trypsin/EDTA must be used to gather single-fibroblast or single-keratinocyte suspensions [27]. Notably, most cell types display a distinctive stress signature related to the dissociation step, which must therefore be taken into account for further downstream analysis [28].
Suction blistering is another approach used to collect skin cells while avoiding extensive enzymatic digestion. High-vacuum mechanical forces are applied to the patient skin in this in vivo procedure, allowing the collection of epidermal and upper dermal cells [29,30]. Recently, a scRNA-seq comparison between suction blistering and normal punch biopsies has been reported. Some cell types were under-represented in suction blistering biopsies; however, both sampling techniques shared most of the identified pathways [30]. In addition, the authors claimed that suction blistering led to a better transcriptomic resolution of skin cells, also presenting the possibility to combine interstitial fluid analysis at the protein level [29,30].
types display a distinctive stress signature related to the dissociation step, which must therefore be taken into account for further downstream analysis [28].

Single-Cell Molecular Profiling of In Vivo Epidermis
Skin diseases can be ascribed to alterations in both basal or suprabasal keratinocytes, to other epidermal cells, or to cell-cell communication networks established in vivo between the different skin layers. In light of this, the single-cell analysis of the entire skin biopsy could be relevant to better understand both healthy and pathological conditions.

Single-Cell Molecular Profiling of In Vivo Epidermis
Skin diseases can be ascribed to alterations in both basal or suprabasal keratinocytes, to other epidermal cells, or to cell-cell communication networks established in vivo between the different skin layers. In light of this, the single-cell analysis of the entire skin biopsy could be relevant to better understand both healthy and pathological conditions.
The first scRNA-seq report of human keratinocyte heterogeneity was published by Cheng et al., who studied epithelial complexity in different anatomical sites by collecting nine biopsies among foreskin, scalp, and truncal skin (Table 1, Figure 1A-C) [31]. They found 11 main clusters, 8 of which were identified as keratinocytes, 2 as melanocytes, and 1 composed of immune cells. Among the eight keratinocyte clusters, they distinguished two basal clusters expressing KRT14 and KRT5, one formed by spinous keratinocytes expressing KRT1 and KRT10 and one granular cluster marked by FLG and LOR. In addition, they identified four other clusters defined as WNT1 (SFRP1 + and FRZB + ); follicular (MGST1 + and S100A2 + ); mitotic (CCND1 + and PCNA + ); and channel (ATP1B3 + and GJB2 + ) ( Table 2).
WNT1 and follicular clusters were mainly present in scalp-derived skin, confirming their role in follicular hair biogenesis, while channel cluster was highly under-represented in all samples. Basal keratinocytes were characterized by the expression of different markers: The basal1 cluster was marked by CXCL14 and DMNK, while the basal2 cluster was marked by CCL2 and IL1R2. Only foreskin biopsies showed a third basal cluster, referred to as basal3, which was enriched for AREG expression, an EGFR ligand that promotes keratinocyte proliferation [31]. The relevance of AREG-EGFR signaling in the proliferative context has also been suggested in a single-cell dataset generated from long-term expanded skin [43].
A comparison between adult trunk-( Figure 1C), fetal-(7-10 weeks post-conception), and atopic dermatitis-or psoriasis-derived skin (Figure 2A) was performed by Reynolds et al., who characterized epidermal non immune, dermal non-immune, antigen-presenting, and lymphoid and mast cells groups for a total of 528,253 cells. Among the adult keratinocytes, the authors identified a cluster of undifferentiated (KRT14 + and KRT5 + ), one of differentiated (KRT1 + and KRT10 + ), and one of proliferating cells(CDK1 + and PCNA + ) ( Table 2), as previously reported [25,31]. Fetal-derived keratinocytes differed from the adult ones due to the expression of KRT8, KRT18, and KRT19. Trajectory analysis suggested that the undifferentiated cluster contains SCs, which give rise to proliferating progenitors and differentiated cells [25].
In addition, both basal clusters are enriched for cell-cycle genes (similar to the basal3 cluster from Cheng et al. [31] and the mitotic cluster from Zou et al. [44], later discussed in this review) even though cell-cycle regression was used in the analysis, meaning that differences in their transcriptome are not only strictly related to the cell cycle itself. BAS-III (ASS + , COL17A1 + , and POSTN + ) and BAS-IV (GJB2 + and KRT19 + ) ( Table 2) clusters expressed genes located at the rete ridges level, as confirmed by immunofluorescence analysis [33].
This model has been confirmed through pseudotime trajectory analysis, in which basal cells give rise to spinous and granular cells, but no clear consensus exists on the trajectory position of the proliferative cluster compared with the quiescent one [33,35,44].
Wiedemann et al. recently studied the differences in thebehavior of cells derived from the sole, the hip, and the palm (Table 1, Figure 1E-G). Particularly, they investigated the common features and differences among palmoplantar and non-palmoplantar areas [38]. The authors analyzed the transcriptomic profile of 15,423 cells, 9471 of which were keratinocytes. Three clusters were identified as basal (basal I-III) (KRT14 + , KRT15, TP63 + , ITGB1 + , and ITGB4 + ), two as spinous (KRT10 + , KRT1 + , SBSN + , and KRTDAP + ) and one as granular (FLG + , LOR + , and SPINK5 + ) ( Table 2). Basal-cell clustering is mostly coherent with the clusters identified in the foreskin by Wang et al. Nonetheless, Wiedemann et al. identified an additional spinous cluster in adult-derived skin, highlighting peculiarities present in the palm, the sole, and the hip. Monocle trajectory analysis confirmed that basal I and II clusters are at the starting point of the pseudotime trajectory, followed by basal III, spinous, and granular clusters [38].

The Wound-Healing Process Analyzed at Single-Cell Level
Martinengo et al. estimated a global pooled prevalence of 2.21 per 1000 people for chronic wounds of mixed etiologies [53]. In chronic wounds, inflammation, fibrosis, and necrosis coexist, leading to a hampered skin regeneration process. Understanding the complex human wound-healing mechanisms would greatly impact the success of diagnosis and therapy of both acute and chronic wounds. ScRNA-seq data generated from chronic wounds represent a powerful tool, offering new perspectives on the healing process [37,41,54].
Comparisons of chronic non-healing pressure ulcer skin (PU), normal acute wounds (AWs), and skin from matched healthy donors have been performed (Table 1, Figure 2B).
A total of 1170 single cells were collected and analyzed using Smart-seq2 technology. In agreement with Cheng et al., keratinocytes were classified into four clusters (KC1-4). KC2 and KC4 were classified as proliferating basal cell clusters: the first was characterized by the expression of genes linked to adhesion and SC signature; the last was characterized by the expression of immunomodulatory genes. The remaining KC1 and KC3 clusters represented the spinous and granular cells, respectively. All the identified clusters were detected in both PU-and AW-derived samples, but changes in their relative abundance were reported. An increase in the percentage of spinous (KC1) and granular (KC3) cells was evident in AW, compared with healthy skin. In contrast, in PU-derived skin, fewer spinous (KC1) and basal (KC2 and KC4) keratinocytes were detected, while a major persistent granular (KC3) cluster was present. Moreover, an increase in the proportion of immune cells was shown in chronic wounds compared with the other samples. Over-representation analysis (ORA) showed that PU keratinocytes overexpressed neutrophil-mediated immunity and apoptosis-related genes [41]. These findings were also confirmed in diabetic foot ulcers, another type of chronic wound, where inflammation and programmed cell death pathways were both upregulated in dermal fibroblasts [50].
In order stratify the patients with PU, the most variable genes on all PU keratinocytes were used to perform PCA. Two groups of patients were identified: t the PU_G2 group was characterized by the enrichment in keratinocytes expressing proliferation-related markers, while keratinocytes from PU_G1 group displayed the upregulation of genes involved in MHC-II-mediated antigen presentation and IFN-γ response. The stimulation of human primary keratinocytes with cell-free wound fluid from PU_G1 patients (but not PU_G2) induced keratinocytes expression of MHC-II-related genes. This effect was blocked by neutralizing IFN-γ in the wound fluids, suggesting that IFN-γ may account for MHC II expression in PU keratinocytes [41].
These different mechanisms in PU highlight the importance of molecular wound diagnosis. Indeed, distinct molecular hallmarks highly correlate with different clinical outcomes that could be exploited in novel personalized medicine strategies concerning chronic wounds. For instance, a targeted molecular diagnosis could highlight the pathogenic mechanisms occurring in a patient's wound. This could enable the repurposing of anti-IFN-γ antibodies previously developed for other disease treatments, to improve wound healing by blocking IFN-γ pathways [41,55,56].
Another therapeutic opportunity was investigated by Singh et al. Their effort focused on the possibility of a pharmacological reversal of DNA hypermethylation, since mouse model studies showed that this might be a feasible solution to rescue tissue regeneration [37,57]. Epithelial-mesenchymal transition (EMT) is responsible for the initiation of re-epithelialization required for wound closure. The loss of this epithelial plasticity leads to chronic wound persistence. The DNA methylation profile of chronic wound edges (WEs) compared with unwounded samples (UWs) revealed that EMT-related genes and their upstream regulator TP53 were hypermethylated in WEs. ScRNA-seq data analysis from 25,168 cells from chronic WEs and 25,561 cells from UW skin allowed researchers to distinguish two keratinocyte clusters, identified as Kera1 (KRT14 + and KRT1 + ) and Kera2 (KRT19 + and KRT7 + ) (Table 1, Figure 2B) [37]. WE samples were marked by the absence of Kera2 cells, which expressed genes relevant to the cell plasticity and the metabolism switch required during the EMT and the preneoplastic progression [58,59]. Differences also arose in the Kera1 cluster, which displayed a lower expression of TP53 target genes in WE-derived cells, probably due to TP53-promoter hypermethylation, as suggested by methylome data. This finding has been further validated in in vitro studies and mouse models [37]. Thus, the presence of TP53-demethylated locus in a diseased hypermethylated genome paves the way for further investigations. New therapeutic strategies might act on the peculiar epigenetic landscape of the wound-compartmentalized microenvironment.

Inflammatory Skin Diseases: Towards Precision Medicine Using Single-Cell Data
Non-communicable inflammatory skin diseases (ncISDs), such as psoriasis or atopic dermatitis, are a major cause of global disease burden due to their frequency, heterogeneity and complexity. Driven by an intricate interplay of genetics and environmental factors, they are a crucial challenge of modern medicine [60,61]. Due to the lack of models able to predict therapeutic responses, many individuals do not benefit from available therapies [46,60]. Thus, new omic approaches could provide a deep phenotyping of key cell types in ncISDs, offering new information to tackle these diseases [46,51]. Currently, most of the studies are focused on psoriasis ( Figure 2C) [25,46,47,51] and atopic dermatitis (AD or eczema) (Figure 2A) [25,40,46,47], while some others also investigate lichen planus (LP) ( Figure 2D) [46], erythrokeratodermia variabilis ( Figure 2E) [47], lupus erythematosus (LE) ( Figure 2F) [42,48], and vitiligo ( Figure 2G) [29].
TNFα and IL-17 cytokines contribute to the dysregulation of the immune response in keratinocyte-driven rashes through NF-kB activation. The systemic blockade of these cytokine pathways is beneficial, but a specific targeting in the accessible skin tissue would reduce systemic side effects. A20 is a promising targetable NF-kB-inhibiting partner protein. The scRNA-seq analysis of 42,105 cells showed that A20 overexpression inhibited the expression of an inflammatory-genesignature upregulated in psoriasis, AD, and erythrokeratodermia variabilis samples. Based on this finding, the A20 in vivo upregulation could represent a therapeutic path to dampen skin inflammation in a variety of ncISDs [47].
The scRNA-seq analysis of psoriasis, atopic dermatitis, and healthy control biopsies enabled the identification of four clusters of undifferentiated, differentiated, proliferating, and inflammatory differentiated cells. The inflammatory-differentiated cluster, expressing inflammation markers (ICAM1, TNF, and CCL20), as well as low levels of undifferentiationrelated (TP63 and ITGA6) and differentiation-related (KRT1 and KRT10) markers, was expanded in psoriasis skin [25]. Previous reports have shown that psoriasis lesional skin was enriched in inflammatory-differentiated keratinocytes and mitotic cells, suggesting an increase in the cell plasticity in disease states [31].
He et al. provided a further single-cell characterization of AD biopsies (Table 1). AD keratinocytes were characterized by the upregulation of epidermal proliferation-associated genes, including S100 and protease inhibitors SERPINB4, consistent with epidermal hyperplasia. The lesional proliferating and suprabasal keratinocytes from their dataset overexpressed KRT6, KRT6A, and KRT16, which are usually enriched in hyperproliferative and wound-healing states [40]. Similar results were confirmed via a scRNA-seq analysis of suction blistering and punch biopsy of both healthy and AD skin [30] (Table 1).
Notably, ncISDs (particularly LP, AD, and psoriasis) have also been studied through spatial transcriptomics using Visium technology by 10× Genomics. From each group of patients, a lesional and non-lesional section was withdrawn [46] (Table 1). This assay is not formally at single-cell resolution, as it generally covers between 1 and 10 cells per spot. Therefore, the authors did not directly analyze single keratinocytes but implemented a density-based clustering method. This allowed them to correlate cytokine expression to responder gene signatures, according to spatial features, paving the way for curative treatment strategies, such as antigen-specific immunotherapies [46], as already suggested for pemphigus vulgaris [62,63].
Although all these strategies might prove successful, further investigations of cell-cell interaction in lesional and non-lesional samples are required, as they might shed light on the subsequent processes that trigger inflammation and lesion formation [42,48]. The identification of the specific disease-causing antigens, and the subsequent targeting of cytokine-producing immune cells in the inflammatory microenvironment, will speed up the development of ncISD-tailored therapies [46].

Single-Cell Molecular Profiling of In Vitro Cultured Human Primary Keratinocytes
Regenerative therapies to restore damaged or diseased epithelium are routinely applied in clinics. Several patients have been successfully treated, completely saving or changing the lives of those affected by extensive burns, skin genetic diseases, and burned-cornea blindness [7,9,24]. The success of all these different approaches relies on the possibility of in vitro culture epithelial SCs. In 1975, H. Green was able to overcome the limitations previously observed in the cultivation of epidermal cells in surface cultures by coculturing them with mitotically inactivated murine fibroblasts [64]. This protocol has been used since 1984 to treat burned patients using autologous cultivated epidermal graft [7].
Given their potency and plasticity, in vitro cultured keratinocytes have been characterized over time. Already in 1987, H. Green was able to describe three clonal types: holoclones, meroclones, and paraclones [65]. Holoclone-forming cells displayed the highest proliferative potential and unique self-renewal capacity; therefore, they are considered epithelial SCs. Conversely, meroclone-and paraclone-forming cells showed less proliferative potential and are normally referred to as TACs. As progenitors, they lack self-renewal capacity and progressively lose their proliferative potential, giving rise to terminally differentiated cells. These assumptions were initially based on clinical observations, in which permanent epithelial regeneration was only possible when an adequate number of holoclone-forming cells were present in the grafted culture [66]. Nonetheless, the formal proof that holocloneforming cells are SCs only became possible in 2017, when a patient affected by JEB was successfully treated using corrected autologous keratinocytes. Clonal tracing analysis has been performed using the provirus integration site as a unique clonal marker, allowing researchers to track cell progeny in the transduced grafted skin [10].
To better understand the biology on which these functional differences rely, the molecular characterization of the three clonal types has been carried out. TP63, YAP, and FOXM1 were found to be the crucial transcription factors (TFs) that sustain epithelial SC selfrenewal [36,[67][68][69]. The microarray analysis of holoclone-, meroclone-, and paraclonederived progenies identified a list of 526 genes differentially expressed in holoclones, as compared to meroclones and paraclones, hence defined as holoclone signature [36].
To gain insight into clone-founding cells, scRNA-seq was applied to the clinical-grade culture of human primary keratinocytes extracted from two healthy donors' truncal skin biopsies (Table 1, Figure 1C) [23] The transcriptomic profiles of 7354 cells were analyzed. Three clusters expressed high levels of clonogenic/basal markers (KRT14 + , TP63 + , ITGA6 + , and ITGB1 + ), whereas differentiated cells, identified as clusters TD1 and TD2, expressed SERPINB3, SFN, KRT10, IVL, and SPINK5. Among the three clonogenic clusters, the one expressing the holoclone signature to the highest extent was named H, containing holoclone-forming cells. M and P clusters expressed lower levels of holoclone signature and contained meroclones and paraclone-forming cells. The H cluster expressed high levels of genes linked to cell cycle, DNA repair, microtubule organization, and YAP and FOXM1 signaling pathways. Monocle analysis confirmed that the H cluster is the starting point of a unique and mainly linear trajectory that proceeds through M and P clusters, giving rise to TD1 and TD2 cells. Singlecell data were further validated using gain-and loss-of-function experiments, confirming the crucial role of FOXM1 as a downstream target of YAP to sustain SC self-renewal potential during in vitro culture [36,67].
The presence of an adequate number of SCs in the graft is one of the main concerns during translational phases [70]. Jayarajan et al. performed a single-cell analysis to study the anoikis-preventing effect of Rho-associated kinase inhibitor (ROCKi) in maintaining keratinocyte SC self-renewal ( Figure 1K, Table 1) [45,71].
ROCKi enhanced cell clonogenic potential after a 6-day treatment. In line with Enzo et al., they identified clusters of SCs, TACs, and differentiated cells in untreated conditions. Surprisingly, they showed a ROCKi-driven reduction in the percentage of holocloneforming cells, reversible upon withdrawal. The single-cell profiling of long-term-treated keratinocytes could shed light on molecular mechanisms responsible for SC reduction, possibly confirming previously published data [72,73].
The epigenetic profiling of in vitro cultured keratinocytes at the single-cell level was published by Khavary's lab. To study the gene network controlling cell fate, they developed perturb-ATAC (assay for transposase-accessible chromatin) (Table 1, Figure 1A), a method able to measure the impact of CRISPR modification on chromatin accessibility in each cell [32]. To this aim, undifferentiated and calcium-induced differentiated keratinocytes were subjected to single-cell ATAC sequencing (scATAC-seq). This allowed the alignment of single cells onto a unique and mainly linearly pseudotime differentiation trajectory. The differentiation was driven by 67 TFs that clustered in 3 modules, which were differentially activated along the differentiation route and were able to recapitulate (i) the proliferation state and progenitors' mitosis control; (ii) the mid-differentiation and the cell-cell adhesion; and (iii) the late keratinization in terminally differentiated cells. The authors identified TFs required for the differentiation program and examined epigenetic changes upon single or combined gene silencing of those TFs. This method could ideally identify co-regulated regulatory elements, providing insights for further biological validations [32].

Network Identification Using Single-Cell Dataset
Studies involving cell-cell interactions might play a pivotal role in understanding physiological and pathological tissue conditions. Starting from single-cell data, different approaches have been used to investigate such synergy in the human epidermis.
Network studies have emerged from the re-analysis of Cheng's dataset, through which a network of transcription factors modules controlling self-renewal and differentiation has been established. However, the biological validation of these results was only partial [74].
Notably, scRNA-seq data provide novel opportunities to study receptor-ligand interactions, identifying networks of communicating cells in tissues. In a previous study by Wang et al., cell-cell interactions were determined using the SoptSC algorithm. They predicted the interactions occurring via WNT, JAK-STAT, NOTCH, and TGF-β signaling pathways. The WNT pathway was active in most of the basal and spinous cells, while TGF-β signaling was restricted to the basal layer. The upper layers of the epidermis were enriched in NOTCH4 and JAK/STAT signaling [33].
Other research groups studied cell-cell interactions in healthy and inflamed skin. Starting from scRNA-seq data, the differential expression of ligand-receptor pairs in AD versus healthy controls was analyzed [40]. CCL2 was significantly upregulated in lesional AD versus control basal keratinocytes. The gene was also abundantly expressed by a unique population of AD COL6A5 + and COL18A1 + inflammatory fibroblasts, unrecognized by previous single-cell studies on healthy dermal fibroblasts [75]. CCL2 receptors (CCR1 and CCR2) were expressed on macrophages and dendritic cells. Another chemokine, CCL27, was upregulated in lesional keratinocytes, not only in basal but also in suprabasal clusters. CCL27 upregulation suggests the establishment of signaling via CCR10 in T cells. Taken together, these findings highlight the interactions between immune and other cell types in skin [40].

Bioinformatic Tools for Cell Network Analysis
Starting from sc-RNA seq data, several user-friendly tools have been recently developed (Table 3). Here, we report those employed to study cell-cell interactions in healthy, wounded, and cancerous skin.
NicheNet, developed in 2019, enables the assessment of not only ligand-receptor interactions but also its predicted impact on intracellular signaling. As a result, it can predict which ligands influence gene expression in another cell, which target genes are affected by each ligand, and which signaling mediators may be involved [76]. Ji et al. applied NicheNet to study keratinocyte-predicted ligands that modulate the transcriptomic profile of squamous cell carcinoma (SCC) microenvironment-specific cells at the leading edge ( Figure 2H). The tumor-specific keratinocyte signaling with nearby cancer-associated fibroblasts was mediated by several receptor-ligand pairs, such as MMP9-LRP1 and TNC-SDC1. The tool also confirmed that tumor keratinocytes at the leading edge resembled an EMT-like population, as they expressed TGFB1 and integrins such as ITGA3 and ITGB1 [34].  Figure 1H). The analysis showed that the number of interactions between fibroblasts and undifferentiated keratinocytes decreased during aging [39].
Released in 2021, the CellChat database considers the additional effects that soluble and membrane-bound stimulatory and inhibitory cofactors exert on these interactions. This tool was employed by Sun et al. to identify intercellular communications driving skin regeneration for long-term expansion therapy [43]. Ligand-receptor pairs AREG-EGFR, CD96-NECTIN, and LAMIN-CD44 were identified as the most significantly upregulated signaling. These pathways might be essential for mechanical stretching and likely contribute to the maintenance of long-term skin regeneration. As expected, the EGF pathway resulted in upregulation in the expanded skin as well [43]. In the work of Guerrero-Juarez et al., CellChat allowed the description of the interactions between fibroblasts and the basal cell carcinoma (BCC) stroma in affecting tumor growth. Fibroblast-secreted WNT5A was identified as the ligand interacting with receptors FZD6 and FZD7 expressed by the basal layer's keratinocytes [49]. Yakupu et al. took advantage of CellChat to provide new insights into PU cellular connections, using data from Li et al. [41,54]. This analysis revealed that intercellular communication is enhanced in both number and strength, mostly between spinous keratinocytes and other clusters in PU, compared with AW samples. Nonetheless, a comparison with healthy skin is missing. The authors suggested that, in PU, spinous and mitotic cells mainly receive signals from melanocytes. In particular, the protease-activated receptor (PAR) signaling pathway mediated by the CTSG-F2RL1 ligand-receptor pair is one of the most activated in PU [54]. CTSG is a serine protease mainly expressed by melanocytes and involved in inflammation, while F2RL1 is widely expressed by fibroblasts and keratinocytes, where it plays an important role as a driver of inflammatory response. These data suggest that PARs might represent an important therapeutic target [54,78].
In 2023, Cang et al. developed a new tool called COMMOT that added spatial information to the study of cell-cell interactions. To validate it, they applied the algorithm to a wide range of datasets, demonstrating that it can consistently capture cell-cell communications (CCCs) already described in the literature [79]. COMMOT was used to study the role of CCCs in human epidermal development. Starting from receptor-ligand pairs annotated in CellChat DB, COMMOT predicted that the molecular interactions between GAS6 and PROS1 with TYRO3 were significant in granular cells and moderately present in basal cells. This notion was confirmed using IF and in situ hybridization (ISH) analysis. Furthermore, the authors analyzed four signaling pathways with well-established roles in epidermal homeostasis, namely WNT, TGF-β, NOTCH, and JAK/STAT. All the signaling cascades were mainly upward-directed, from the basal to suprabasal layers. Conversely, some signals were downward-directed toward the basal layers at the bottom of the ridges [79,80].

Conclusions
Single-cell techniques have revolutionized the field of biology, offering unprecedented insights into the complexity and heterogeneity of biological systems. Undoubtedly, the ability to profile individual cells has greatly expanded the understating of tissue organization and dynamics, unraveling the mechanisms that underlie normal development and disease progression. Furthermore, the development of multiomic approaches facilitates the description of the genomic and gene expression profile within each cell, providing a multifaceted and comprehensive understanding of cellular behavior.
However, it is important to be cautious of potential pitfalls. Considering the enormous amount of information derived from these techniques, standardizing bioinformatic analysis is crucial for the generation of comparable data. Misleading interpretations could also arise from sample-handling procedures, as they directly influence cell collection and stressrelated cellular responses. Therefore, in silico findings need careful experimental validation to avoid computational and sample-handling-derived artifacts.
Skin is a complex tissue, and its cellular heterogeneity is now being tackled by singlecell analysis, in both physiological and pathological conditions. The ability to characterize the epidermis architecture at the single-cell level, combined with cell-cell communication models, offers powerful insights into its many layers of biological complexity.
Despite the potentialities of these emerging tools, some biological questions remain unanswered. Indeed, further analyses are required to elucidate which pathways are involved in in vivo self-renewal and symmetric/asymmetric cell division. A starting point in deciphering epidermal stem-cell molecular profile involves in vitro experiments, in which they are activated in a condition mimicking wound healing.
The importance of unraveling skin complexity is linked to the possibility of (i) monitoring clinically relevant cell populations for advanced therapy (e.g., epidermal SCs in cell and gene therapy applications); (ii) defining a patient-specific diagnosis; and (iii) providing targets to improve precision medicine.
Excitingly, new frontiers are emerging in the field of single-cell analysis. Among them is the possibility to deconvolute cell complexity at the protein level, which promises an even more accurate characterization of biological processes.