Multiscale spatial mapping of cell populations across anatomical sites in healthy human skin and basal cell carcinoma

Significance Single-cell RNA sequencing (scRNAseq) has revolutionised cell biology, enabling high-resolution analysis of cell types and states within human tissues. Here, we report a comprehensive spatial atlas of adult human skin across different anatomical sites and basal cell carcinoma (BCC)— the most common form of skin cancer—encompassing in vivo optical coherence tomography, scRNAseq, global spatial transcriptomic profiling, and in situ sequencing. In combination, these modalities have allowed us to assemble a comprehensive nuclear-resolution atlas of cellular identity in health and disease.

Our understanding of how human skin cells differ according to anatomical site and tumour formation is limited.To address this, we have created a multiscale spatial atlas of healthy skin and basal cell carcinoma (BCC), incorporating in vivo optical coherence tomography, single-cell RNA sequencing, spatial global transcriptional profiling, and in situ sequencing.Computational spatial deconvolution and projection revealed the localisation of distinct cell populations to specific tissue contexts.Although cell populations were conserved between healthy anatomical sites and in BCC, mesenchymal cell populations including fibroblasts and pericytes retained signatures of developmental origin.Spatial profiling and in silico lineage tracing support a hair follicle origin for BCC and demonstrate that cancer-associated fibroblasts are an expansion of a POSTN+ subpopulation associated with hair follicles in healthy skin.RGS5+ pericytes are also expanded in BCC suggesting a role in vascular remodelling.We propose that the identity of mesenchymal cell populations is regulated by signals emanating from adjacent structures and that these signals are repurposed to promote the expansion of skin cancer stroma.The resource we have created is publicly available in an interactive format for the research community.Human skin comprises two main layers, the epidermis and dermis.It is home to cells of the innate and adaptive immune systems and is extensively vascularised and innervated.Skin in different body sites performs specialised functions that are reflected in differences in the nature and density of adnexal structures-hair follicles (HFs), sebaceous glands, and sweat glands.As a result of its exposure to ultraviolet (UV) light and other carcinogens, the skin is susceptible to developing cancer, and basal cell carcinoma (BCC) is the most common type of cancer affecting humans (1).
In order to gain a fuller understanding of skin cell organisation and function in health and disease, a number of research teams, including our own, have used single-cell RNA sequencing (scRNAseq) to create a comprehensive catalogue of the cell types present in healthy adult skin, developing skin and common inflammatory skin diseases (2).However, we currently lack a full understanding of whether the cell types in skin from different anatomical sites differ, and we do not have high-resolution spatial maps of where specific cell types are located.
The issue of cellular identity is particularly interesting because during development the dermis of the trunk and limbs is derived from the somatic mesoderm, whereas the dermis of the anterior head and neck is derived from the cranial neural crest (3).The question of spatial cellular location is also important since lineage tracing studies in mice have shown that func tionally distinct subpopulations of dermal cells have different spatial organisation (4)(5)(6).
A further important question is how the cell populations present in healthy skin relate to those in skin cancer (7).BCC is a locally invasive skin cancer composed of islands of tumour cells embedded in a dense fibroblastic stroma along with vascular elements that preferentially affects the face (8).The cell of origin of BCC has long been a subject of debate, with some studies pointing to a HF origin (9,10).Growth of a BCC requires both expansion of stromal elements and the growth of new blood vessels (11).Altered proliferation of epithelial elements is driven primarily by cell-intrinsic genetic alterations (12), including PTCH1, TP53, NOTCH1, NOTCH2, and FAT1 mutations in the case of BCC (13).However, the origins and drivers of the associated stromal changes are not fully understood.Fibroblasts within tumour stroma are referred to as "cancer-associated fibroblasts" (CAFs) (14) and are con sidered a potential target for therapeutic strategies (15), but it is not clear whether they have a distinct cellular identity compared to healthy dermal fibroblasts.
The goals of the Human Cell Atlas ( 16) are to define the states of all cell populations within tissues in both health and disease and to understand their spatial localisation with respect to one another.Here, we report a comprehensive spatial atlas of cell populations in human skin across multiple scales incorporating Visium spatial transcriptomics (ST) at 55 μm resolution, in situ sequencing (ISS) of transcripts at subcellular resolution, <10-μmresolution optical coherence tomography (OCT) imaging of live human skin, and scRNAseq of skin from multiple body sites and BCC.Our analysis has allowed us to create maps of cellular populations at single-cell resolution for healthy skin at different ana tomical sites and for BCC.We believe that this will be a valuable resource for the community.

Cartography of a Multiscale Atlas of Healthy Human Skin and BCC.
To understand the landscape of skin cell populations, we studied healthy skin across multiple anatomical sites and BCC (Fig. 1A).The human skin exhibits substantial morphological variation between sites, such as differences in epidermal thickness, HF, sebaceous gland and eccrine coil density (SI Appendix, Fig. S1 A and B).
Previous scRNAseq studies of human skin have largely focused on the skin of the trunk and limbs (2,(17)(18)(19)(20).However, BCC is more common in facial skin (7,8).We therefore augmented existing publicly available datasets through scRNAseq of full thickness skin biopsies from 15 healthy donors and eight BCCs across five facial skin sites (ear, nose, cheek, forehead, and temple).Additionally, we used two complementary approaches to localize cells within tissue sections: global ST (Visium) and targeted ISS (21,22) (Fig. 1A).We also quantified skin vasculature structure with angiographic OCT across anatomical sites.OCT images of in vivo human skin obtained at 10 μm resolution through the principle of low-coherence interferometry (23) yielded high-resolution 3D images of cutaneous vasculature (Fig. 1A and SI Appendix, Fig. S1C).OCT showed dif ferences in the dermal density of vascular networks between face and body skin (SI Appendix, Fig. S1D).
scRNAseq in Healthy Skin and BCC Samples.We sequenced 233,379 cells, of which 134,839 cells passed quality control (QC).We have integrated our sequenced cells with two publicly available datasets from 11 healthy body sites (Inguinal, Arm) to give a combined dataset of 155,401 cells (Fig. 1B and SI Appendix, Fig. S2 A-C).Unsupervised clustering revealed a total of 30 clusters corresponding to 16 cell types (Fig. 1B and SI Appendix, Fig. S2D).Cell types and sub-cell types were annotated with previously described markers (2,17,19,(24)(25)(26) (Fig. 1C and SI Appendix, Fig. S2E and Table S1).
The mesenchymal cells that constitute the dermis of the face are derived from neural crest (3), whereas those of the body are derived from somatic mesoderm.Previous studies have revealed develop mental gene expression differences in fibroblasts at different ana tomical sites in both human and mouse (27,28).We confirmed that adult skin cells retained signatures of their embryological ori gins: we found differences in expression of Hox genes and mesen chymal neural crest markers between face, including facial BCC, and body skin in scRNAseq and spatial transcriptomic datasets (SI Appendix, Fig. S3 A and B), particularly for fibroblasts, pericytes, SMC (smooth muscle cells), vascular endothelial cells, and chon drocytes (SI Appendix, Fig. S3 C and D).These findings demonstrate that signatures associated with developmental origin can be detected in a cell type-specific manner in adult face and body skin.
Although we found skin cell populations to be conserved (Fig. 1B), we identified a mesenchymal population that was pres ent only in a full thickness ear sample.On the basis of marker expression (including high expression of ACAN and COL2A1), we identified these cells as extra-articular chondrocytes (Fig. 1 B-D and SI Appendix, Fig. S2 E and F).We additionally identified a skeletal muscle cell population that was present only in a fore head sample (Fig. 1 B -D and SI Appendix, Fig. S2 E and F).The face contains 20 distinct flat skeletal muscles that attach to the skull, including the forehead (29).The chondrocytes and skeletal muscle cells are not strictly skin cell types and were not represented across all the samples; therefore, we did not study them further.
Differences in the relative abundance of cell types were observed.T cells were overrepresented in facial skin and BCC compared to body skin.Natural Killer/Innate Lymphoid Cells (NK/ILC) and B cells were overrepresented in BCC compared to healthy skin (face and body sites).NK/ILC were also more abundant in face compared to body sites.Fibroblasts were more abundant in body skin compared to facial skin.A difference in the abundance of keratinocytes was also observed between face and body skin, which likely reflects in the efficiency of keratinocyte isolation (Fig. 1D and SI Appendix, Fig. S2D).

ISS.
Having catalogued cell populations across anatomical sites and in BCC, we localized them within the tissue by combining scRNAseq data with global ST and targeted ISS (21,22).
We prepared ST libraries for nine body, 13 face, and eight BCC sections (Fig. 1A and SI Appendix, Fig. S4A).The presence of high-quality RNA throughout sections was verified through house keeping gene expression via multiplex RNA FISH (SI Appendix, Fig. S4B).The percentage of mitochondrial genes and Unique Molecular Identifier counts in the Visium sections confirmed that the libraries were of high quality (SI Appendix, Fig. S4 C and D), and the few low-quality spots were removed.To deconvolve cell types present within each Visium 55 μm grid spot, we used a published algorithm, cell2location ( 22)-a Bayesian model that can resolve fine-grained cell types in ST data to create comprehensive cellular maps of diverse tissues, permitting the inference of cell types present on the basis of an annotated scRNAseq dataset (Fig. 1 A and B).
An advantage of ST is the ability to analyse the entire transcrip tome in whole tissue sections.However, resolution is limited to 55-μm grid spots, and highly expressed genes can dominate sequencing libraries in these spots.To compensate for this, we performed ISS.This approach utilises rolling-circle amplification in combination with barcoded padlock-probe circularization (30).
Based on published scRNAseq data (2,17,19,20,(24)(25)(26) and our sequencing of facial skin, we defined a panel of 165 marker genes representing all skin cell populations (SI Appendix, Tables S2-S4).The presence of high-quality RNA was verified through housekeep ing gene expression via the ISS workflow (SI Appendix, Fig. S4E).We then performed ISS of four body, four face, and four BCC sections (Fig. 1A and SI Appendix, Fig. S4 A and F).According to Qian et al. (21), 50 optimally chosen markers were sufficient for accurate classification of 28 cell populations.Since the number of cell populations present in our human skin dataset is comparable to their analysis, our panel of 165 marker genes should be adequate for cell typing.To do so, we used the Spot-based Spatial cell-type Analysis by Multidimensional mRNA density estimation (SSAM) pipeline (31), which allows inference of cell types from mRNA sig nals in a cell segmentation-free manner by correlating gene expres sion with an annotated scRNAseq dataset.To demonstrate the power of the two computational spatial methods, we predicted and pro jected the main skin cell types such as the keratinocyte clusters, fibroblast clusters, endothelial cell clusters, pericyte clusters, and SMC clusters onto the ST sections and the localisation of ISS reads (SI Appendix, Fig. S5 A-C).
We next applied the same methodologies to all the 30 clusters found in our dataset (Fig. 1B and SI Appendix, Fig. S2 D and E and Table S1) to create a spatial skin atlas.To facilitate the use of the atlas, we generated a public web interface (https://spatial-skinatlas.cellgeni.sanger.ac.uk/) (32), where researchers can explore spatial gene expression patterns and skin cell population predictions.This resource contains an explorer called CELLxGENE for scRNA-seq data and ST data and another explorer called Vitessce (33) for ISS data.

The Architecture of Cutaneous Blood Vessels Reflects Vessel Size
and Anatomical Site and Is Altered in BCC.In healthy skin, the cutaneous vasculature is organised as a superficial (in the papillary dermis) and a deep (in the reticular dermis) vascular plexus connected by perforating vessels.In BCC, tumour epithelial cells are embedded in a dense vascular and fibroblastic stroma.Growth of the tumour to a macroscopic size requires both expansion of stromal elements and the growth of new blood vessels (11).Since the cutaneous vasculature is a dynamic 3D network, histological sections do not accurately capture its structure in vivo.Therefore, we used angiographic (speckle contrast) OCT (34) to generate 3D images of vasculature density across facial and body sites from the same 16 healthy individuals (Fig. 2 A and B) and of matched healthy and BCC areas from the same 11 patients (Fig. 2 C and D).
In healthy skin, vascular density was significantly higher in facial skin compared to body skin in both papillary and reticular dermis (Fig. 2 A and B and SI Appendix, Fig. S1 C and D).However, we found no significant differences in segment length, number of segments, tube thickness, or number of branching points in face skin samples (SI Appendix, Fig. S1E).
In BCC patients, OCT images showed a disorganised vascular plexus with tortuous vessels (Fig. 2C and SI Appendix, Fig. S6A), as previously published (35,36).Quantification revealed a signif icant increase in blood vessel density in the tumour compared to adjacent skin from the same patient (Fig. 2D).
We next analysed the vascular cell populations present in our integrated scRNAseq dataset.In contrast to previous studies (20), we were able to identify only one population of VEC (vascular endothelial cells) and one of LEC (Fig. 1B).Regarding perivascular cell populations, we found three populations of pericytes (Fig. 2E).In keeping with previous studies (20), there were two major pericyte populations, which we denoted RGS5+ and TAGLN+ based on high expression of those genes (SI Appendix, Fig. S2E).We additionally identified a population expressing high levels of DES and ACTA2, which is a gene signature of SMC (SI Appendix, Fig. S2E).Gene ontology (GO) analysis revealed SMC to have signatures of muscle and contractile function, whereas RGS5+ pericytes exhibited signa tures of blood vessel development and morphogenesis.TAGLN+ pericytes had signatures of both blood vessel development and con tractile function (SI Appendix, Fig. S6B).
Despite the increased density of blood vessels in facial skin, we did not find a significant difference in the abundance of VEC, LEC (Fig. 1D) or pericyte subpopulations (Fig. 2F) between healthy face and body skin.However, there was a selective expan sion of the RGS5+ pericytes and a reduction in TAGLN+ pericytes in BCC compared to healthy skin (face and body) (Fig. 2F).
Immunostaining of RGS5 in papillary and reticular dermis showed RGS5+ cells to be in close contact with NCAM+ cells, a pan vascular marker, in healthy skin and BCC (Fig. 2G).TAGLN+ cells colocal ized and were in close contact with NCAM+ cells in the papillary and reticular dermis and in small and large blood vessels in healthy face and body skin.The colocalization of TAGLN+ pericytes with small blood vessels was lost in BCC (Fig. 2G).In healthy skin, DES and ACTA2 RNA hybridization was observed in the arrector pili muscle (APM) cells in proximity to HFs.In BCC, SMC (ACTA2+ DES+ cells) were diffusely distributed within the stroma (Fig. 2G).RGS5, TAGLN, and ACTA2 distributions were confirmed in the global ST data and in ISS (SI Appendix, Fig. S6 C and D).
In healthy skin, ST deconvolution analysis revealed VEC to be localized in superficial and deep vascular plexuses (Fig. 2H).Both RGS5+ and TAGLN+ pericytes colocalized with VEC and with one another throughout the dermis (Fig. 2H).In BCC, the distri bution of RGS5+ pericytes closely paralleled VEC, whereas the distribution of TAGLN+ pericytes was more focal and less abundant (Fig. 2H).This is in contrast to healthy skin.In both healthy skin and BCC, SMC were adjacent to HFs and in the lower dermis, suggesting an association with the deep vascular plexus (Fig. 2H).Pearson correlation coefficients (PCC) of spot normalized cell abun dances using the cell2location predictions in all the ST samples confirmed a strong colocalization of RGS5+ and TAGLN+ pericytes with VEC and with one another.Although fewer in number, TAGLN+ pericytes in BCC still colocalized with VEC (Fig. 2I).
We also applied CellPhoneDB (37), a receptor-ligand interaction analytical tool; this revealed a high number of predicted interactions (>100) between RGS5+ and TAGLN+ pericytes and VEC (SI Appendix, Fig. S6E).Significant biological ligand-receptor inter actions are predicted (SI Appendix, Fig. S6F) for RGS5+ pericytes with VEC and TAGLN+ pericytes with VEC, as those predicted cell clusters are part of the same spatial microenvironment in ST data.
At the spatial single-cell level, RGS5+ and TAGLN+ pericytes were also localized to both small and large blood vessels.SMC were restricted to larger vessels and additionally colocalized with the APM (Fig. 2J and SI Appendix, Fig. S6 C and D).
Fibroblast Subpopulations Are Spatially Restricted within the Dermis of Healthy Skin and BCC Stroma.Previous studies in both mouse (4, 5) and human (2,6,17,19,20,24,25) have revealed the presence of multiple fibroblast subpopulations in the dermis.The number of subpopulations found in human skin varies between studies and spatial information at the tissue level and level of anatomical sites is lacking.
Unbiased clustering of our integrated scRNAseq dataset iden tified four main populations of fibroblasts which we designated APOD+, SFRP2+, PTGDS+, and POSTN+ based on specific gene expression (Fig. 3A).The same number of fibroblast clusters was present in face, body, and BCC skin.However, there were significant differences in their abundance and distribution.
APOD+ fibroblasts were overrepresented in healthy facial skin in comparison with body skin and BCC (Fig. 3B).They highly coexpressed APOD and APOE, also expressed in immune cells (SI Appendix, Fig. S2E).GO analysis revealed potential functions in cell migration and mobilisation, characteristic of immune cells and implicated in humoral immune responses (SI Appendix, Fig. S7A).APOD in situ hybridization coupled with vimentin (VIM) immunostaining showed an association of APOD+ fibro blasts and blood vessels (VIM+) (Fig. 3C).The APOD distribution is shown in global ST data and ISS data (SI Appendix, Fig. S7  B and C).ST deconvolution analysis and PCC confirmed APOD+ fibroblasts to be colocalized with VEC in facial and body areas and in BCC (Fig. 3 D and E).Receptor-ligand analysis revealed a high number of predicted interactions (>100) between APOD+ fibroblasts and VEC (SI Appendix, Fig. S7 D and E).The increased prevalence of APOD+ fibroblasts in facial skin (Fig. 3B) is in keep ing with the increased vascular density compared to body skin (Fig. 2B and SI Appendix, Fig. S1 C and D).
SFRP2+ fibroblasts were overrepresented in body sites compared to facial sites and BCC (Fig. 3B).They coexpressed high levels of SFRP2 and WISP2 (SI Appendix, Fig. S2E).GO terms included extracellular matrix (ECM) structure and organisation, but also negative regulation of BMP signalling, which is known to control HF growth (38); these are characteristics of reticular fibroblasts (SI Appendix, Fig. S7A).Moreover, WISP2 is a Wnt inhibitor (39).Epidermal Wnt signalling is important in the regulation of HF development and cycling (40).SFRP2 RNAscope spots were present  throughout the dermis (Fig. 3C).The distribution of SFRP2 was confirmed in global ST data and ISS data (SI Appendix, Fig. S7 B  and C).ST deconvolution analysis confirmed SFRP2+ fibroblasts to be present throughout the dermis (Fig. 3D).The greater abun dance of SFRP2+ fibroblasts in body skin (Fig. 3B) could be related to denser collagen and the lower quantity of epithelial structures in the body compared to face (SI Appendix, Fig. S1 A and B).
The abundance of the PTGDS+ fibroblasts was similar between face and body skin and BCC (Fig. 3B).GO terms included ECM organisation and remodelling but also collagen fibril organisation and response to TGF-ß, which are characteristics of papillary fibro blasts (SI Appendix, Fig. S7A).PTGDS RNA spots were located in the papillary dermis in healthy face and body areas and around the tumour islands in BCC (Fig. 3C).The PTGDS distribution is shown in global ST data and ISS data (SI Appendix, Fig. S7 B and C).ST deconvolution confirmed PTGDS+ fibroblasts to be located in papillary dermis and in proximity to follicular structures including eccrine coils (Fig. 3D).Furthermore, computational prediction in ISS showed PTGDS+ fibroblasts to be frequently localized next to the upper HF areas such as bulge and infundibulum (Fig. 3F).
POSTN+ fibroblasts showed a trend of higher abundance in face skin compared to body skin and were significantly increased in BCC compared to healthy skin (Fig. 3B).Like PTGDS+ fibroblasts, GO terms showed function in ECM deposition and organisation and collagen fibril organisation but also characteristics of endoderm formation and development (SI Appendix, Fig. S7A), which is inter esting given that cancer progression shares some gene expression profiles with developmental pathways (41).POSTN RNAscope spots were mainly around HFs in body and face skin and were considerably increased around the tumour islands in BCC (Fig. 3C).Expression of POSTN in ST and ISS confirmed these results (SI Appendix, Fig. S7 B and C).ST deconvolution and PCC analysis showed POSTN+ fibroblasts to colocalize with PTGDS+ fibroblasts in both healthy body and face skin (Fig. 3 D and E).Computational prediction in ISS revealed a localization next to the lower shaft and bulb of HF in healthy skin (Fig. 3F).Receptor-ligand analysis revealed a high number of predicted interactions (>150) between POSTN+ and PTGDS+ fibroblasts (SI Appendix, Fig. S7D); mul tiple collagen subtypes and integrins are implicated.This suggests reciprocal cell-cell communication involving different fibroblast subpopulations within the HF niche (SI Appendix, Fig. S7D).
ST deconvolution and ISS projection confirmed that the expanded POSTN+ fibroblasts were localized in proximity to tumour islands (Fig. 3 D and F).Other fibroblast subpopulations were also present in the stroma but distant from the tumour islands.

Lack of Distinct Epithelial Cell Populations in BCC.
The intimate association of POSTN+ fibroblasts with epithelial tumour islands could potentially be mediated by an inductive signal arising from the epithelial cells.We hypothesised that since the POSTN+ fibroblasts were associated with HFs in healthy skin and expanded in BCC, the epithelial compartment of BCC would have characteristics of abnormal HFs.To test this, we explored the epithelial cell populations in BCC compared to healthy epidermis.
To be able to analyse keratinocytes in the interfollicular epider mis (IFE) and in follicular structures, we enriched our datasets with scRNAseq of 6217 cells from IFE and pilosebaceous units (PSU) microdissected from healthy scalp, an area containing many follicular structures, of which 4,565 cells passed the QC steps (SI Appendix, Fig. S8A).In order to perform high-resolution clus tering and refine cell annotations, we integrated our sequenced cells with one publicly available scRNAseq dataset derived from healthy scalp skin (20,561 epithelial cells; SI Appendix, Fig. S8B) (18).We used those cells as a template for the initial steps of clustering (Fig. 4A) and then analysed only our dataset comprising 11,966 epithelial cells from healthy skin (face and body) and BCC.
Unbiased clustering identified a total of 10 cell populations comprising undifferentiated (basal) and differentiated (spinous and granular) keratinocytes (K) specific to the IFE, dividing and transitional K, and K specific to adnexal structures such as upper HF K, inner and outer bulb K, sebocytes, and sweat gland secre tory luminal cells (Fig. 4 B and C and SI Appendix, Fig. S8C and Table S5).We plotted key markers in ISS (Fig. 4D).
We did not find any distinct epithelial cluster characteristic of BCC.However, markers in the hedgehog signalling pathway such as PTCH1/2, HHIP, EPCAM, a diagnostic marker for BCC (42), and IGKC, a immunologic marker of solid cancer (43), were highly ex pressed in BCC basal keratinocytes compared to healthy (SI Appendix, Fig. S8D).Of known epidermal stem cell markers (SI Appendix, Fig. S9A), TP63 and LGR5 were overexpressed in BCC (SI Appendix, Figs.S8D and S9B).We found an increase of outer bulb K in BCC compared to healthy face skin.We also found an increase in transi tional K in BCC and in healthy face compared to body skin.Surprisingly, there was also an increase of IFE basal and spinous K in healthy body compared to healthy face and BCC (Fig. 4E).ST deconvolution revealed epithelial clusters to localize to distinct regions within the epidermis and cutaneous appendages (SI Appendix, Fig. S8D).Tumour regions exhibited signatures of the inner bulb K and, predominantly, of the outer bulb K (Fig. 4F).Additionally, BCC cells in the inner and outer bulb K clusters showed high expression of PTCH1/2, HHIP (hedgehog signalling pathway), EPCAM and IGKC compared to healthy cells in these same clusters (Fig. 4G).
The most striking difference in BCC versus healthy epidermis was high expression of KRT17 in the nodules and invasive tumours (Fig. 4D).This was confirmed in the scRNAseq data and independently in the ST data (Fig. 4H and SI Appendix, Fig. S8 F and G).Immunostaining of KRT17 showed high expression in the inner bulb and inner root sheath and low expression in the outer bulb and outer root sheath in healthy skin (Fig. 4I).These results lead us to hypothesise that BCC tumour cells retain the transcriptional signatures of their cells of origin.
We then performed in silico single-cell trajectory and pseudotime gene analysis using Monocle 3 (44) in BCC keratinocytes (Fig. 4 J and K) and independently in all epithelial cells (SI Appendix, Fig. S8H).This trajectory inference method revealed numerous "root nodes" in the outer bulb cluster in BCC (Fig. 4J).Pseudotime showed the given lineage for each cell based on the distance from the cell of origin, such as outer bulb K in the BCC (Fig. 4K).
Collectively, our data suggest that malignant epithelial cells of BCC could arise from cells of the inner and outer HF bulb.

Discussion
Through a combination of OCT, scRNAseq, global ST, and tar geted ISS, we have created a cellular-resolution atlas of cell popu lations across multiple anatomical sites in healthy human skin and BCC.Using computational spatial mapping, we probabilistically inferred the localisation of cell populations in tissue sections.Mesenchymal cell subpopulations exhibited characteristic spatial distributions within the dermis, and these same cell populations were repurposed in the stroma of BCC.
We found that the same skin cell populations were present in body and facial skin despite differences in morphology and embryological origin.However, in keeping with a previous study (27), we observed differences in Hox gene expression in fibro blasts derived from the body versus the face.We additionally showed differences in expression of neural crest markers, not only in fibroblasts but also in other cell types including pericytes, Schwann cells, VEC, and melanocytes.These findings indicate that fibroblast, pericyte, and VEC differentiation programmes can be concurrently activated in dermal lineages derived from cranial neural crest and somatic mesoderm.This is of particular interest in regard to oncogenic specificity.We know that onco genic alterations to DNA are not transforming in all cellular contexts (7).It has been shown that anatomic position determines oncogenic responses due to preexisting transcriptional pro grammes of the cell of origin (45).
In line with previous studies, we identified two main pericyte subpopulations (2,(17)(18)(19)(20) and a SMC cluster that also had a perivascular location.RGS5+ and TAGLN+ pericytes corresponded to the populations reported by Reynolds et al. (20).GO and the location of TAGLN+ pericytes suggest a role in blood vessel devel opment and contractility, similar to mesh/thin-strand pericytes previously described in other organs (46).In BCC, there was expansion of the RGS5+ pericytes and a reduction of APOD+ fibroblasts and TAGLN+ pericytes.The distribution of the RGS5+ pericyte subpopulation closely paralleled with VEC, whereas the distribution of the APOD+ fibroblasts and TAGLN+ pericytes was more focal.This contrasts with healthy skin where APOD+ fibro blasts, RGS5+ and TAGLN+ pericytes populations largely colo calize, implying selective expansion of the RGS5+ pericytes during tumour angiogenesis.Genetic ablation of RSG5 in a mouse cancer model reduces tumour angiogenesis (47).These observations sug gest that RGS5+ pericytes may have a role in vascular remodelling during cancer neovascularization and therefore in enhancing tumour progression.
We identified four main fibroblast subpopulations localized to the following spatial contexts: in association with blood vessels (APOD+ fibroblasts); in the interstitial dermis with a concentra tion in the reticular dermis (SFRP2+ fibroblasts); concentrated adjacent to HFs, near the upper HF and in the papillary dermis (PTGDS+ fibroblasts); and near the HF bulb (POSTN+ fibro blasts).While the number of fibroblast subpopulations identified in previous studies has varied, our demonstration of distinct spatial contexts for each of the populations sets a minimum bound on this number.It is probable that additional fibroblast populations will be identified in the future by enriching scRNAseq datasets in mesenchymal cells through increasing the depth of reads and improvement of spatial transcriptomic technologies.
In nodular and infiltrative BCCs, the CAFs corresponded to fibro blast subpopulations in healthy skin.However, there was selective expansion of the POSTN+ fibroblasts which were present adjacent to the tumour islands.This suggests that malignant epithelial cells that retain transcriptional signatures of the HF may induce POSTN+ fibroblasts or promote their selective expansion.POSTN has been previously found in invasive BCCs (48,49).
The spatial localization of mesenchymal populations relative to epithelial appendages or epithelial tumours suggests that reciprocal signalling is important in the establishment and main tenance of cellular identity.We identified multiple subpopula tions of keratinocytes within the skin that correlate with analogous subpopulations in adult mouse skin (50), although we did not find subpopulations specific to the hair bulge, but instead found two subpopulations which localize to the hair bulb.This could reflect differences in human skin architecture and the fact that 90% of human HFs are in anagen compared to mouse HFs (51).
The cell of origin of human BCC has been a subject of debate.Lineage tracing in the mouse supports a follicular origin (10,52) and keratin expression in human BCC is in keeping with this model (9).KRT17 was up-regulated in BCC.In mice, KRT17 is implicated in hair cycle regulation and tissue repair; genetic abla tion inhibits tumour development and growth (53,54).KRT17 is also overexpressed in other cancers suggesting a more general role in tumorigenesis (55).In addition, Lgr5+ stem cells in the lower bulge have been shown to play a role in BCC formation (52,56).We observed that LGR5 was highly expressed in outer and inner bulb K and dividing K in BCC compared to healthy epidermis suggesting a role of those specific stem cells in human BCC.However, murine studies have shown that activation of oncogenic hedgehog signalling in bulge stem cells is not able to induce BCC (57), whereas activation in the IFE can (58).Spatial prediction of keratinocyte subpopulations and in silico lineage tracing support a follicular origin of BCC, specifically the bulb area.
An understanding of the location and potential function of cell populations in healthy skin and how they are subverted in skin cancer is likely to inform therapeutic approaches.This could include inhibition, ablation, or selective expansion of specific fibroblast or pericyte subpopulations.The rational design of such therapeutics will require a mechanistic understanding of the intrinsic transcriptional regulators of these cell populations and the mechanisms through which reciprocal signaling is coordinated with neighbouring cells.
In summary, our cellular-resolution map of human healthy skin and cutaneous malignancy should be a valuable resource for the skin research and dermatology communities and forms a corner stone of the Human Cell Atlas.

Materials and Methods
All the materials and methods used in this study, including histology, immunohistochemistry, RNAscope, OCT imaging, tissue dissociation, scRNAseq, ST (10X Visium), ISS, and data analysis are described in detail in SI Appendix, Supplementary Materials and Methods.The key resources are summarised in a key resource table in SI Appendix, Supplementary Materials and Methods.
Ethics.The study was sponsored by Guy's and St Thomas' NHS Foundation Trust and King's College London and was subject to both institutional and external research ethics council (REC) review (REC reference 19/NE/0063).Following informed consent, excess skin samples were obtained from patients undergoing skin surgery.
Data, Materials, and Software Availability.The raw and processed datasets for scRNAseq and Visium ST generated in the current study are available on ArrayExpress under "E-MTAB-13085" (59) and "E-MTAB-13084" (60) accessions and an interactive tool to display and download scRNAseq, ST and ISS data is available at https://spatial-skin-atlas.cellgeni.sanger.ac.uk/ (32).Previously published data were analysed in this work (17,19).All other data are included in the manuscript and/or SI Appendix.
human cell atlas | skin | basal cell carcinoma | single cell RNA sequencing | fibroblasts

Fig. 1 .
Fig. 1.Cartography of a multiscale human skin atlas across anatomical sites and cancer.(A) ScRNAseq was performed for skin samples derived from multiple facial sites and compared with previously published scRNAseq of body sites.We also profiled skin samples from these conditions with global ST which permitted us to examine the global transcriptional profile at 55 μm resolution; ISS which mapped skin at subcellular resolution with a custom panel of 165 genes (SI Appendix, Tables S2-S4) designed to represent all the cellular populations present in human skin; and OCT imaging which permitted us to assess the morphology of the skin at 10 μm resolution in vivo.The number of samples for each modality is indicated.The location and identity of cells were inferred probabilistically in the global ST data and in the ISS data according to our annotated scRNA-seq dataset through computational analyses.(B) Uniform Manifold Approximation and Projection (UMAP) and clustering of 155,401 cells from 33 donors (11 healthy body sites, 14 healthy face sites, and 8 BCC patients) representing 30 skin cell populations and 16 cell types.Each cell is represented by a single point with the two-dimensional location of this point on the UMAP plot corresponding to the global transcriptional state of the cell.Cells with similar transcriptional profiles form clusters. (C) Dot plot showing well-known marker genes specific to each skin cell type.For each cell type, the percentage of cells expressing the marker (diameter) and the average log2 normalized expression (color) are shown.(D) Abundance of skin cell types in healthy face and body areas and in BCC from the face.The relative proportion of sequenced cells assigned to each cellular population is illustrated for body, face, and BCC face.Significant differences (one-way ANOVA test) are indicated (*<0.05,**<0.01, and ****<0.0001).

Fig. 2 .
Fig. 2. Healthy morphology of cutaneous vasculature and localisation of vascular cell populations and their alterations in the stroma of BCC.(A) 3D morphology of vascular networks is compared for representative OCT scans of facial and body skin from the same individuals.OCT images (500 frames) illustrate reflectivity (grayscale) overlaid with the blood flow (red) (Upper).The skin vasculature network is also shown as a 3D computational reconstruction (Lower).Image size 6 mm × 6 mm; (scale bar: 500 μm.) (B) Quantification of skin vascular density across facial and body sites from the same individuals.Vascular density was quantified for a total of six sites (3 facial sites and 3 body sites) across 16 individuals.Statistically significant differences (one-way ANOVA test) are indicated (****<0.0001).(C) Comparison of vascular architecture in healthy skin areas and BCC sites from the same patients.Image size 6 mm × 6 mm (scale bar: 500 μm).(D) Quantification of skin vascular density in healthy and BCC sites from the same patients.Vascular density was quantified for a total of three facial sites across 11 individuals.Statistically significant differences (one-way ANOVA test) are indicated (***<0.001).(E) UMAP plot of pericytes (RGS5+ pericytes, TAGLN+ pericytes), SMC, and VEC subpopulations present in healthy skin (face and body) and BCC.(F) Comparison of the abundance of different pericyte subpopulations (as a proportion of total pericytes) in healthy facial and body skin and BCC.Statistically significant differences (one-way ANOVA test) are indicated (ns, not significant, **<0.01,***<0.001,and ****<0.0001).(G) Immunostaining of NCAM, RGS5, and TAGLN in the papillary and reticular dermis in healthy body and face skin and in BCC.RNAscope imaging of ACTA2 and DES in healthy body and face skin and in BCC (DAPI, blue).(H) VEC and pericyte cell populations annotated in scRNAseq data were computationally predicted on global ST sections using cell2location.Predicted cell abundances shown by color gradients per spot in tissue architecture (H&E) images (representative samples are shown).(I) Cell cluster colocalization analysis for VEC with RGS5+, TAGLN+ pericytes and SMC.The barplot shows PCC of cell2location predictions per spot normalized cell abundances across all spots of Visium samples (each individual bar represents a Visium sample).(J) Computational spatial mapping of VEC and pericyte subpopulations in ISS sections using SSAM (representative samples are shown).Localization of basal and suprabasal keratinocytes is shown in order to indicate skin structures.

Fig. 3 .
Fig. 3. Localisation of fibroblast populations in healthy human skin and expansion of the POSTN+ subpopulation in the stroma of BCC.(A) UMAP plot of four main fibroblast subpopulations in human healthy skin and BCC (APOD+ fibroblasts; SFRP2+ fibroblasts; POSTN+ fibroblasts; PTGDS+ fibroblasts).(B) Comparison of the abundance of fibroblast subpopulations (proportion of total fibroblasts) between body and facial skin and BCC face.Statistically significant differences (one-way ANOVA test) are indicated (ns, not significant, **<0.01,***<0.001,and ****<0.0001).(C) RNAscope of specific markers of each fibroblast subpopulation: APOD, SFRP2, PTGDS, POSTN (red), DAPI (blue) and immunostaining of VIM in healthy skin and BCC tissue sections.(D) Fibroblast populations annotated in scRNAseq data were computationally predicted on global ST sections using cell2location.Predicted cell abundances shown by color gradients per spot in tissue architecture (H&E) images (representative samples are shown).(E) Cell cluster colocalization analysis for VEC with APOD+, SFRP2+, PTGDS+ and POSTN+ fibroblasts and POSTN+ fibroblasts with APOD+, SFRP2+, and PTGDS+ fibroblasts.The barplot shows PCC of cell2location predictions per spot normalized cell abundances across all spots of Visium samples (each individual bar represents a Visium sample).(F) Computational spatial mapping of fibroblast subpopulations in ISS sections using SSAM (representative samples are shown).Localization of basal and suprabasal keratinocytes is shown.

Fig. 4 .
Fig. 4. Conservation of epithelial cell populations between normal skin and BCC.(A) UMAP of body, face skin (including IFE and PSU cells), and BCC integrated with a publicly available healthy scalp skin dataset.(B) Unsupervised clustering identified a total of 10 clusters.(C) Dot plot showing average log2 normalized expression of specific marker genes for each epithelial subpopulation.(D) Spatial localisation of ISS reads for the different epithelial subpopulations in ISS images from representative healthy body, face, and BCC face samples (DAPI, gray).KRT14, POSTN, DLL1 represented the IFE basal K; Dividing K coexpressed CDK1 and PCNA; KRT17 was a marker of the inner bulb cluster.SFRP1; WIF1 marked the outer HF bulb.KRT1, KRT10, and LY6D marked spinous IFE keratinocytes; KRT19 and DCD marked sweat glands; MGST1 marked sebaceous glands.Markers such as FOS, JUN or KRT6A and KRT6B, highly expressed by transitional K and upper HF K, were not included in our ISS gene panel.(E) Comparison of the abundance of the epithelial subpopulations in body, facial skin, and BCC.Statistically significant differences (one-way ANOVA test) are indicated (**<0.01 and ****<0.0001).(F) Cell2location prediction of keratinocyte populations corresponding to inner and outer bulb clusters.Predicted cell abundances shown by color gradients per spot.Three representative samples are shown.(G) Volcano plot highlighting differential expression of PTCH1/2, HHIP, EPCAM, and IGKC in BCC compared to healthy cells in the inner and outer bulb K clusters.(H) Differential expression of KRT17 in scRNAseq data (Left) and global ST data (Right) between BCC and healthy skin (face and body).Statistically significant differences (one-way ANOVA test) are indicated (***<0.001).(I) Immunostaining of KRT17 in healthy face skin and BCC (DAPI, blue).(J) Single-cell trajectory gene analysis using Monocle 3 showing root nodes in BCC epithelial cells.(K) Monocle pseudotime showing given lineage for each cell based on the distance from the cell of origin in BCC epithelial cells.