Skip to main content

Single-cell analysis of human prepuce reveals dynamic changes in gene regulation and cellular communications



The cellular and molecular dynamics of human prepuce are crucial for understanding its biological and physiological functions, as well as the prevention of related genital diseases. However, the cellular compositions and heterogeneity of human prepuce at single-cell resolution are still largely unknown. Here we systematically dissected the prepuce of children and adults based on the single-cell RNA-seq data of 90,770 qualified cells.


We identified 15 prepuce cell subtypes, including fibroblast, smooth muscle cells, T/natural killer cells, macrophages, vascular endothelial cells, and dendritic cells. The proportions of these cell types varied among different individuals as well as between children and adults. Moreover, we detected cell-type-specific gene regulatory networks (GRNs), which could contribute to the unique functions of related cell types. The GRNs were also highly dynamic between the prepuce cells of children and adults. Our cell–cell communication network analysis among different cell types revealed a set of child-specific (e.g., CD96, EPO, IFN-1, and WNT signaling pathways) and adult-specific (e.g., BMP10, NEGR, ncWNT, and NPR1 signaling pathways) signaling pathways. The variations of GRNs and cellular communications could be closely associated with prepuce development in children and prepuce maintenance in adults.


Collectively, we systematically analyzed the cellular variations and molecular changes of the human prepuce at single-cell resolution. Our results gained insights into the heterogeneity of prepuce cells and shed light on the underlying molecular mechanisms of prepuce development and maintenance.

Peer Review reports


The prepuce, commonly referred to as the foreskin, comprises a specialized tissue structure that encapsulates the glans penis in males [1,2,3]. Despite its evident biological and physiological significance, our understanding of the prepuce's cellular heterogeneity and functional diversity remains limited. With the advent of single-cell RNA-sequencing (scRNA-seq) technologies, gene/transcript expression profile is feasible to be investigated at a single-cell resolution [4, 5]. Moreover, scRNA-seq also enables the investigation of cellular expression heterogeneity, inference of gene regulatory networks, and construction of cell-to-cell interactions [6,7,8,9]. Although several studies have investigated the transcriptional profile of human epidermis (e.g., truncal skin, and scalp) at the single-cell level [10,11,12], the cellular composition and heterogeneity of human prepuce were still largely unknown. Single-cell analysis of human prepuce is essential for unraveling its complex tissue composition and physiological functions.

Furthermore, with the development of prepuce, its cell composition and gene expression profile could be dynamically changed. The expression heterogeneity of prepuce cells is largely correlated with the composition of cell subtypes, while expression profiles of cell subpopulations are directly influenced by the gene regulatory networks (GRNs) formed by transcription factors (TFs) and downstream target genes [13,14,15]. Additionally, different cell types usually interact with each other to exert corresponding functions, which are mediated by related ligands and receptors [16,17,18,19,20]. Thus, cellular differences in human prepuce between children and adults could be closely correlated to the variations of GRNs and cell–cell interaction networks. Systematic investigation of the prepuce's cellular dynamics and molecular changes could provide a deeper understanding of the role of different cell types in maintaining prepuce health and preventing infections. It could also add valuable knowledge to the fields of human anatomy, biology, and reproductive health.

Here, we first explored the expression heterogeneity and cell type compositions of human prepuce of children and adults based on related scRNA-seq data. Then, we constructed the gene regulatory networks (GRNs) for the prepuce cells and investigated the gene regulatory dynamics among transcription factors and downstream target genes. The GRN differences between children and adults were also examined. Additionally, we built the cellular interaction networks among different cell subtypes and interrogated cell–cell communication variations between children and adults. The signaling pathways enriched by the ligand-receptor pairs of cell–cell interactions were explored as well.


Cell type identification for human prepuce based on scRNA-seq data

To explore the cell composition and expression heterogeneity of human prepuce at single-cell resolution, we first applied scRNA-seq to prepuce samples from 5 children and 4 adults. After removing the low-quality samples and cells (two unqualified samples of one child and one adult were excluded), a total of 90,770 qualified single cells were kept for 4 children and 3 adults (Fig. 1A). Based on the top 2000 variable genes, those prepuce cells from children and adults were mixed together in the UMAP plot (Fig. 1A). All these 7 samples showed similar cell distributions, suggesting that no batch effect exists among different samples. Using a graph-based clustering approach, the 90,770 cells of these 7 samples could be grouped into 27 different clusters with Seurat [21] (Fig. 1B). A series of marker genes with enriched expression were identified for different clusters (Fig. 1C, Supplementary Table S1). According to the known marker genes of corresponding cell types in CellMarker database [22], those 27 clusters could be further classified into 15 cell subtypes, including dendritic cells (DC), smooth muscle cells (SMC), fibroblast (Fibro), vascular endothelial cell (VEC), vascular endothelial-to-mesenchymal transition cell (VEndMT), lymphatic endothelial cell (LEC), lymphatic endothelial-to-mesenchymal transition cell (LEndMT), melanocyte (MC), keratinocyte/epithelium (KC/Epi), schwann cell (SC), Mast, Macrophages (Macro), neutrophils (Neutro), T/natural killer (NK) and B cells (Fig. 1D). Interestingly, fibroblast cells (55.52%) occupied the largest portion of cells, followed by SMC cells (10.49%), T/NK cells (9.50%), VEC (5.15%), Macro (4.93%), and LEC (4.83%).

Fig. 1
figure 1

Cell clustering and cell-type identification of human prepuce based on scRNA-seq data. A UMAP plot showing 90,770 qualified prepuce cells from 4 children and 3 adults based on the top 2000 variable genes in expression. B Clustering results of the prepuce cells of children and adults. C Expression heatmap for the markers of 27 different cell clusters. D UMAP plot displaying the 15 distinct cell subtypes of prepuce cells

Characterization of marker expression and cell composition in the human prepuce

For each of those 15 cell types, we identified a set of marker genes with enriched expression compared to other cell types (Fig. 2A, adjusted p-value < 0.05). For instance, DCN, LUM, and COL1A2 were highly expressed in fibroblast cells, while CD3D, CD3E, and CD3G showed enriched expression in T/NK cells (Fig. 2A). PECAM1, VWF, and ACKR1 were detected as the marker genes for VEC, while ACTA2, MYL9, and MYH11 exhibited high expression in SMC. CPA3, TPSB2, and CTSG were the marker genes of MAST, while SFN, KRT1, and KRT14 were with enriched expression in KC. Specifically, as shown in Fig. 2B-E, some marker genes like CD3D, CDH19, CD1C, and CPA3 were almost only expressed in T/NK, SC, DC, and Mast cells, respectively.

Fig. 2
figure 2

Maker expression profile and cell subtype composition of prepuce cells. A Expression distribution of selected markers for each type of prepuce cells. B-E Specific makers for T/NK, SC, DC, and Mast cell types. F Cell type composition of the prepuce in each individual of children and adults

We further examined the cellular compositions of different prepuce cell subtypes in children and adults. For each individual of children and adults, fibroblast cells accounted for the largest portion of cells (children: 46.46%-66.18%; adults: 52.6%-60.78%). However, the ranking of other types of cells varied greatly among different individuals (Fig. 2F). For example, smooth muscle cells were the second most abundant cells for child-2 (11.85%), child-3 (17.51%), child-4 (10.37%), and adult-1 (13.11%), while T/NK cells occupied the second largest fraction for child-1 (19.98%) and adult-2 (12.68%). Interestingly, we also found that the proportions of all the 15 cell types were highly variable among the 4 children in general, whereas several cell types showed relatively constant small fractions (< 1%) among the 3 adults (e.g., B cells, KC, Mast, neutrophils, and VEndMT). These results suggest that human prepuce of both children and adults have prominent differences in composition ratios of cell types.

Gene regulatory network inference for prepuce cells

To investigate the transcriptional regulatory profile of human prepuce, we inferred the gene regulatory networks (GRNs) in prepuce cells by employing SCENIC [23]. According to the expression associations between transcription factors (TFs) and downstream target genes, a total of 59 significant regulons were detected based on all 90,770 prepuce cells (each regulon contains one TF and its downstream target genes). These regulons involved 59 TFs and 2985 downstream target genes in total. As shown in Fig. 3A, different cell types can be distinguished in the t-SNE plot based on those 59 regulons, suggesting that the gene regulatory profiles of these cell types could be significantly different.

Fig. 3
figure 3

Gene regulatory network profile of human prepuce cells. A t-SNE plot showing the distribution of 15 cell types of prepuce based on 59 significant regulons. The cells in the plot were colored by corresponding cell types. B-F Profiles of cell-type specific regulons for KC, SC, and B cells, respectively. Venn plots displaying the intersection between downstream target genes of corresponding TF and the marker genes of related cell types. Bubble diagrams showing the enriched biological processes or KEGG pathways for the downstream target genes of the corresponding TF

Intriguingly, we detected a set of cell-type-specific regulons that mainly activated in certain cell types. For example, the regulon formed by TF KLF5 was primarily detected in KC, which contained 64 downstream target genes (Fig. 3B). We observed that 40 out of those 64 (62.5%) downstream target genes of KLF5 were the marker genes of KC (e.g., EHF, KRT14, CALML5, SERPINB5, PKP1, and LGALS7), indicating the coregulation of these genes and the significant constribution of KLF5 to the expression profile of KC. Those 64 target genes of KLF5 were significantly enriched in the biological processes of epidermis development, skin development, keratinocyte differentiation, and epidermal cell differentiation (Fig. 3B), further demonstrating its important regulatory role in the development of human prepuce. The regulons formed by TF SOX2 (66 downstream target genes) and SOX10 (113 downstream target genes) were enriched in the cell subpopulations of SC (Figs. 3C and D). Of note, 35 out of those 66 (53%) SOX2 target genes (e.g., L1CAM, COL9A3, CADM4, AATK, MPZ, and ATP1A2) and 55 out of those 113 (48.7%) SOX10 target genes (e.g., GPR155, L1CAM, CNPY2, GFRA3, COL9A3, and CADM4) were the marker genes of SC. The enriched KEGG pathways for those 66 SOX2 downstream target genes were cell adhesion molecules, PI3K − Akt signaling pathway, and ECM − receptor interaction, while those 113 SOX10 downstream target genes were mainly enriched in cell adhesion molecules and Wnt signaling pathway (Fig. 3C and D). The regulons of TF PRDM1 (18 downstream target genes) and XBP1 (99 downstream target genes) were mainly detected in B cells (Fig. 3E and F). Compared with the B cell marker genes, 6 out of those 18 (33.3%) PRDM1 target genes (e.g., RAB30, TMEM156, PDK1, MIR155HG, FBXW7, and PRDM1) and 29 out of those 99 (29.3%) XBP1 downstream target genes (e.g., C16orf74, IRF4, TPD52, RAB30, CD79B, and DERL1) overlapped with those markers. Gene functional enrichment analysis showed that those 18 PRDM1 target genes were mainly involved in the biological processes of positive regulation of proteolysis, negative regulation of triglyceride metabolic process, regulation of fibroblast apoptotic process, and fibroblast apoptotic process, while those 99 XBP1 target genes were mainly enriched the KEGG pathways of protein processing in the endoplasmic reticulum, biosynthesis of nucleotide sugars, amino sugar and nucleotide sugar metabolism, and fructose and mannose metabolism (Fig. 3E and F). Therefore, those cell-type-specific regulons could play crucial roles in regulating the gene expression of corresponding cell subpopulations. The results would benefit the understanding of cellular expression heterogeneity and the function of different cell subtypes of the human prepuce.

Dynamics of gene regulatory networks in human prepuce cells

To gain insights into the transcriptional regulatory differences between children and adults, we further separately constructed the GRNs for these two groups. As shown in Figs. 4A and B, each cell type could be distinguished from other cell types based on the detected regulons for children and adults. Interestingly, we found that a few TFs formed potential age-specific regulons, which were also differentially expressed between children and adults. For example, TF androgen receptor (AR) formed a specific regulon (involved 19 downstream target genes, such as BNC2, BOC, and CRNDE) in children, which showed a more enriched expression in children compared to adults (Fig. 4C). AR is a steroid-hormone-activated TF and has the potential to affect cellular differentiation and proliferation in corresponding tissues [24]. TF KLF10 also showed higher expression in children than adults and formed a children-specific regulon (involved 19 downstream target genes, such as CNOT8, NFYA, and PHF21A) (Fig. 4D). Previous studies have shown that KLF10 might play a critical role in regulating the circadian clock [25, 26]. Accordingly, we identified potential children-specific regulons that were not detected in adults, these regulons could be important in regulating the expression of related gene sets in children's prepuce cells.

Fig. 4
figure 4

Gene regulatory differences between the prepuce cells of children and adults. A t-SNE plot of prepuce cells for children based on detected regulons. B t-SNE plot of prepuce cells for adults based on inferred regulons. C Expression profile of children-specific regulon AR. D Expression distribution of children-specific regulon KLF10. E–F Regulons with the same TF but having different sets of downstream target genes in children and adults. t-SNE plot showing the activation state of a regulon in cells (blue: activated; gray: un-activated). Venn plots displaying the intersection between the downstream target genes of the same TF-regulon in children and adults. The significantly enriched biological processes were shown in the bubble diagram

We also observed that several TFs formed corresponding regulons in both children and adults, but they activated in different cell types and the downstream target genes varied greatly. For instance, TF JUND formed a specific regulon in B cells of children (containing 81 downstream target genes), whereas it generated a different regulon (including 51 downstream target genes) specific to the T/NK cells of adults (Fig. 4E). Most of the downstream target genes regulated by JUND in children were different from those regulated in adults (Fig. 4E), only sharing 8 genes (CCND2, TNFRSF12A, ANXA6, JUND, BTG1, ODC1, DUSP8, and PRR7). Gene functional enrichment analysis showed that those downstream target genes regulated by JUND in children were mainly involved in the biological processes of negative regulation of phosphorylation, regulation of protein catabolic process, negative regulation of protein phosphorylation, and regulation of protein ubiquitination, while those JUND targeting genes in adults were enriched in peptidyl-threonine dephosphorylation and inactivation of MAPK activity. The regulon of TF MITF was specific to MAST cells of children (having 34 downstream target genes), whereas it formed a specific regulon (including 227 downstream target genes) in the adult MC (Fig. 4F). We found that the downstream target genes regulated by MITF had 27 common genes (such as CLCN7, SLCO4A1, GNPTAB, GNAL, and MSC) between children and adults (Fig. 4F). Those 34 MITF downstream target genes in children were mainly enriched in the biological processes of pigmentation, developmental pigmentation, regulation of cell shape, pigment cell differentiation, and cellular pigmentation, while those 227 downstream target genes regulated by MITF in adults were involved in similar biological processes as well (adjusted p-value < 0.05). Consequently, we revealed the dynamics of gene regulatory networks formed by the same TFs in prepuce cells between children and adults, which could contribute to the expression dynamics of human prepuce development and maintenance.

Cell–cell communication network construction for human prepuce cells

To investigate the cellular communications among different cell types in children and adults, we constructed the cell–cell interaction networks based on the expression profiles of ligand-receptor pairs using CellChat [27]. By comparing the cell–cell interaction networks between children and adults, we found that fibroblast, VEndMT, and LEndMT cells exhibited relatively strong cellular communications with other cell types, whereas the cell types of T/NK, B, MC, Mast, KC, and Neutro showed relatively weak cell–cell interactions with other types of cells for both children and adults (Fig. 5A and B). However, cell types of VEC, LEC, Macro, SMC, DC, and SC generally had stronger cellular interactions with other types of cells in children compared to adults. For the outgoing cellular interaction strength, the largest was fibroblast in children followed by LEndMT, VEndMT, and LEC, whereas the strongest was LEndMT in adults followed by fibroblast, VEndMT, and DC (Supplementary Figure S1A). For the intensity of incoming cell–cell communications, the highest was VEC in children followed by VEndMT, Macro, and DC, while the largest was VEndMT in adults followed by LEndMT, VEC, and DC (Supplementary Figure S1B). For children, the lowest outgoing and incoming cellular interaction strength were KC and B cells respectively. In comparison, the weakest outgoing and incoming cell–cell interaction intensity were both Mast cells in adults.

Fig. 5
figure 5

Cell–cell communication networks among different prepuce cell types in children and adults. A Cell–cell interaction heatmap between any two cell types of children prepuce. B Cellular communication heatmap between cell types for the prepuce of adults. C Signaling pathways that have stronger interaction intensity in children compared to adults. D Signaling pathways that show stronger interaction intensity in adults compared to children

We also found that the cellular interaction intensities for some enriched signaling pathways formed by corresponding ligand-receptor pairs were significantly different between children and adults (Fig. 5C and D). For example, stronger cell–cell interactions were observed between B cells and DC or Macro cells in children compared to adults for BAFF signaling pathway (Fig. 5C). Our result was in line with previous reports that the cytokine of BAFF was critical for supporting the survival of mature naïve B cells, which was also required for the survival of autoimmune B cells and memory B cells [28,29,30]. More cellular communications were also detected among VEndMT, LEC, and VEC in children than that in adults for the EPHB signaling pathway (Fig. 5C). Previous studies have shown that the EPHB signaling pathway plays an important role in cell adhesion and migration, indicating that it could promote the development of children's prepuce [31,32,33,34].

In contrast, some signaling pathways showed opposite trends between children and adults (Fig. 5D). For instance, the cell–cell interactions enriched in the ADGRE5 signaling pathway were stronger in adults compared to children. ADGRE5 encodes the proteins of the EGF-TM7 subfamily of adhesion G protein-coupled receptors, which may be functionally important for cell adhesion as well as the recruitment, activation, and migration of leukocytes [35,36,37]. A similar phenomenon was observed for the ANGPT signaling pathway, where higher cell–cell interaction intensity was detected among VEndMT, LEndMT, SMC, Macro, LEC, and VEC in adults compared to children. These signaling pathways could be important for prepuce maintenance in adults, reflecting the functional differences of cellular communications between children and adults.

Cell–cell communication variations among different cell types of the human prepuce

We also detected the signaling pathways that were mainly enriched in children and adults, respectively. For example, the signaling pathways of CD96, EPO, IFN-1, CHEMERIN, and WNT were only identified as significant enrichment in children (Fig. 6A), suggesting that they could be mainly activated in children rather than adults. Cell–cell interactions were primarily detected between KC and T/NK cells for the CD96 signaling pathway. It has been suggested that CD96-mediated signaling had the potential to modulate the differentiation of effector T cells, which could have a co-stimulatory role in the activation and effector function of CD8 + T cells [38]. EPO signaling pathway was mainly enriched by the cell–cell interactions between fibroblast and other cell subtypes. Previous studies have revealed that EPO could enhance the differentiation of myofibroblasts, and also had the potential to accelerate skin wound closure [39,40,41]. CHEMERIN signaling pathway was mainly enriched by the interactions between Macro and other cell types, while IFN-I and WNT signaling pathways were mediated by the cell–cell communications among various cell types in the prepuce of children. Thus, these signaling pathways could be functionally important in the prepuce development of children.

Fig. 6
figure 6

Enriched cell-type specific signaling pathways of children and adults. A Signaling pathways mediated by ligand-receptor pairs significantly enriched in children. B Signaling pathways significantly enriched in adults mediated by ligand-receptor pairs

The signaling pathways mainly detected in adults included BMP10, NEGR, ncWNT, PRL, and NPR1 (Fig. 6B). BMP10 signaling pathway was enriched by the cellular communications among the cell types of VEndMT, LEndMT, fibroblast, LEC, VEC, and T/NK, while the NEGR signaling pathway was mainly activated by the interactions with SC. The activation of the NPR1 signaling pathway in adults’ prepuce was primarily mediated by the cell–cell interactions between VEndMT and other cell types. By contrast, ncWNT and PRL signaling pathways were enriched by the cellular interactions among different cell types in adults. The cell–cell communications mediated by corresponding ligand-receptor pairs for these pathways could play an important role in the prepuce maintenance of adults. Our results indicate that the cell–cell interaction networks among different cell types varied greatly between children and adults.


To the best of our knowledge, we are the first to systematically explore the cellular heterogeneity of human prepuce between children and adults from different aspects (including expression, gene regulation, and cell–cell communication) at the single-cell level. A total of 15 different cell types were identified in human prepuce cells and each cell type had a set of marker genes with significantly enriched expression. Fibroblast cells accounted for the largest fraction in both children and adults. We observed that the proportion of those identified cell types in different individuals varied greatly in general. A number of differentially expressed genes (DEGs) between children and adults for each cell type were detected. The quantity of DEGs for those 15 cell types ranged from 3 to 981 (Supplementary Figure S2), indicating the large expression variation between the prepuce cells of children and adults.

We also revealed the dynamics of gene regulation among different cell types as well as between the prepuce of children and adults. In total, 59 significant regulons were detected based on the 90,770 prepuce cells, which were involved in 59 TFs and 2985 downstream target genes. Interestingly, we detected a number of cell-type-specific regulons with enriched activation for certain cell types, such as the KC-specific regulon formed by TF KLF5, the regulons of TFs SOX2 and SOX10 specific to SC, as well as the regulons of PRDM1 and XBP1 enriched in B cells (Fig. 3). Surprisingly, a large fraction of those downstream target genes regulated by these TFs were the marker genes detected for corresponding cell types, suggesting that cell-type-specific TFs played critical gene regulatory roles and could significantly affect the expression profile of related cell types. We further identified two regulons formed by TFs AR and KLF10 that were mainly activated in children prepuce compared to adults. Moreover, we also uncovered several regulons shared between children and adults but activated in different cell types and regulated distinct sets of downstream target genes. For instance, TF JUND formed a specific regulon in B cells of children but its regulon was enriched in T/NK cells of adults, while MITF regulon was mainly activated in MAST cells of children but it formed a specific regulon of MC in adults (Fig. 4E). Consequently, the prepuce cells showed high gene regulatory dynamics between children and adults. Our findings could facilitate a better understanding of expression heterogeneity of prepuce cells during development.

We constructed the cell–cell communication networks among different cell types of the human prepuce and discovered large-scale variations between children and adults. Strong cell–cell interactions were observed among Fibroblast, VEndMT, and LEndMT cells, while relatively weak cellular communications were detected among T/NK, B, MC, Mast, KC, and Neutro. Furthermore, we also observed that some enriched signaling pathways exhibited highly different cellular interaction intensities between children and adults, such as BAFF, ADGRE5, EPHB, and ANGPT signaling pathways. Additionally, we identified a number of signaling pathways specific to children (e.g., CD96, EPO, IFN-1, CHEMERIN, and WNT signaling pathways) and adults (e.g., BMP10, NEGR, ncWNT, PRL, and NPR1 signaling pathways). Our results not only showed the cellular communication variations among different cell types between children and adults, but also revealed the signaling pathways that could contribute to the prepuce development in children and the prepuce maintenance in adults.

Given the results obtained from our study, several avenues for future investigations could be conducted. One promising direction is the incorporation of spatial transcriptomics techniques to complement our understanding of the prepuce's cellular organization and interactions within its microenvironment. Spatial transcriptomics can provide insights into the spatial distribution of cell types and their gene expression patterns [42,43,44], providing additional information for cell–cell communication networks and tissue organization within the prepuce. On the other hand, the ages of the four children (4, 4, 5, and 5) are very similar in this study, the age differences of adults (19, 32, and 37) may introduce certain biological variations. Further investigation of the longitudinal changes in cell heterogeneity across more groups with similar ages could facilitate better understanding of the development and maintenance of human prepuces. Additionally, exploring the prepuce's cellular dynamics in disease states may provide valuable insights into its tissue homeostasis and susceptibility to pathological conditions.

In conclusion, we dissected the prepuce cells in terms of cellular compositions, gene expression changes, gene regulatory dynamics, and cell–cell communication variations, which shed light on the cellular heterogeneity and underlying molecular mechanisms of human prepuce development and maintenance.

Materials and methods

Quality control of scRNA-seq data and cell type identification

The scRNA-seq approach of 10X Genomics was applied to the prepuce samples from 5 children and 4 adults. Two samples from one child and one adult were excluded due to the low quality of scRNA-seq data. The remaining qualified scRNA-seq data of 4 children (ages of 4, 4, 5, and 5) and 3 adults (ages of 19, 32, and 37) were used in the downstream analysis. We further employed Seurat (version 4.1.1) [21] to conduct cell quality control, dimensionality reduction, and clustering. Those cells with > 5,000 or < 200 expressed genes, > 10% mitochondrial genes, or > 0.01% red blood cell genes were removed. Then we used the function of FindVariableFeatures in Seurat to select the top 2000 variable genes in expression for dimensionality reduction. The function of FindClusters in Seurat was applied to group cells into different clusters with Uniform Manifold Approximation and Projection (UMAP). We determined the cell types of prepuce with corresponding cluster markers using the CellMarker database [22]. Specifically, the top 10 marker genes detected in each cell cluster were compared with the known makers of related cell types in CellMarker database. According to the overlap of marker genes, the cell clusters were grouped into corresponding cell types.

Gene functional enrichment analysis

We employed the R package of clusterProfiler [45] (version 3.18) to carry out Gene Ontology (GO) and KEGG pathway enrichment analysis. Adjusted p-value < 0.05 was used to define the significantly enriched biological processes and KEGG pathways. Gene set enrichment analysis (GSEA) was based on the Molecular Signatures Database (

Inference of single-cell gene regulatory networks

Single-cell gene regulatory networks of cells formed by transcription factors (TFs) and their downstream target genes were inferred using SCENIC (version 1.1.2) [23]. Specifically, the normalized gene expression matrix of cells was first used as the input for the R package of GENIE3 (version 1.16.0) to construct co-expression networks of genes. Next, RcisTarget (version 1.14.0) was employed to deduce the regulatory networks between TFs and downstream target genes. After the construction of gene regulatory networks, AUCell (version 1.16.0) was utilized to analyze the activity of predicted regulons in each cell. Each regulon was formed by the TF and its downstream target genes. Adjusted p-value < 0.05 was applied to define the significance of regulons.

Cell–cell communication network construction and signaling pathway analysis

We used CellChat (version 1.4.0) [27] to construct the cell–cell interaction network among different cell types of human prepuce. Cellular communications between two different cell types or within the same cell types were predicted based on the expression profile of known ligand-receptor pairs. CellChat computed an enrichment score for each potential ligand-receptor interaction by comparing the joint expression of the ligand in the sending cell type and the receptor in the receiving cell type against their individual background expression distributions. Statistically significant interactions between cell types were identified based on empirical p-values using a permutation-based test. By summarizing the probabilities of ligand-receptor pairs, the communication probability of a signaling pathway is computed by CellChat. Using manifold learning and quantitative comparisons, CellChat classified signaling pathways and identified conserved and context-specific pathways for children and adults.

Availability of data and materials

All raw and processed sequencing data generated in this study have been submitted to the National Omics Data Encyclopedia (NODE, database under accession number OEX020875.


  1. Haseebuddin M, Brandes SB. The prepuce: preservation and reconstruction. Curr Opin Urol. 2008;18:575–82.

    PubMed  Google Scholar 

  2. Romeus L, Liu C, Stanton R, Ellsworth P. Routine circumcision? The role of prepuce in syndactyly repair. J Pediatr Urol. 2020;16:497–9.

    PubMed  Google Scholar 

  3. Benson M, Hanna MK. Prepuce sparing: Use of Z-plasty for treatment of phimosis and scarred foreskin. J Pediatr Urol. 2018;14:545.e541-545.e544.

    Google Scholar 

  4. Haque A, Engel J, Teichmann SA, Lonnberg T. A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications. Genome Med. 2017;9:75.

    PubMed  PubMed Central  Google Scholar 

  5. Chen G, Ning B, Shi T. Single-Cell RNA-Seq technologies and related computational data analysis. Front Genet. 2019;10:317.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. Chen G, Schell JP, Benitez JA, Petropoulos S, Yilmaz M, Reinius B, Alekseenko Z, Shi L, Hedlund E, Lanner F, et al. Single-cell analyses of X Chromosome inactivation dynamics and pluripotency during differentiation. Genome Res. 2016;26:1342–54.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. Cao J, Packer JS, Ramani V, Cusanovich DA, Huynh C, Daza R, Qiu X, Lee C, Furlan SN, Steemers FJ, et al. Comprehensive single-cell transcriptional profiling of a multicellular organism. Science. 2017;357:661–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  8. Potter SS. Single-cell RNA sequencing for the study of development, physiology and disease. Nat Rev Nephrol. 2018;14:479–92.

    CAS  PubMed  PubMed Central  Google Scholar 

  9. Yuan D, Tao Y, Chen G, Shi T. Systematic expression analysis of ligand-receptor pairs reveals important cell-to-cell interactions inside glioma. Cell Commun Signal. 2019;17:48.

    PubMed  PubMed Central  Google Scholar 

  10. Cheng JB, Sedgewick AJ, Finnegan AI, Harirchian P, Lee J, Kwon S, Fassett MS, Golovato J, Gray M, Ghadially R, et al. Transcriptional programming of normal and inflamed human epidermis at single-cell resolution. Cell Rep. 2018;25:871–83.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. Nakamizo S, Dutertre CA, Khalilnezhad A, Zhang XM, Lim S, Lum J, Koh G, Foong C, Yong PJA, Tan KJ, et al. Single-cell analysis of human skin identifies CD14+ type 3 dendritic cells co-producing IL1B and IL23A in psoriasis. J Exp Med. 2021;218(9):e20202345.

    CAS  PubMed  PubMed Central  Google Scholar 

  12. Kim D, Chung KB, Kim TG. Application of single-cell RNA sequencing on human skin: technical evolution and challenges. J Dermatol Sci. 2020;99:74–81.

    CAS  PubMed  Google Scholar 

  13. Chan TE, Stumpf MPH, Babtie AC. Gene regulatory network inference from single-cell data using multivariate information measures. Cell Syst. 2017;5:251-267.e253.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Jackson CA, Castro DM, Saldi GA, Bonneau R, Gresham D. Gene regulatory network reconstruction using single-cell RNA sequencing of barcoded genotypes in diverse environments. Elife. 2020;9:e51254.

    CAS  PubMed  PubMed Central  Google Scholar 

  15. Iacono G, Massoni-Badosa R, Heyn H. Single-cell transcriptomics unveils gene regulatory network plasticity. Genome Biol. 2019;20:110.

    PubMed  PubMed Central  Google Scholar 

  16. Armingol E, Officer A, Harismendy O, Lewis NE. Deciphering cell-cell interactions and communication from gene expression. Nat Rev Genet. 2021;22:71–88.

    CAS  PubMed  Google Scholar 

  17. Almet AA, Cang Z, Jin S, Nie Q. The landscape of cell-cell communication through single-cell transcriptomics. Curr Opin Syst Biol. 2021;26:12–23.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Li Y, Hu X, Lin R, Zhou G, Zhao L, Zhao D, Zhang Y, Li W, Zhang Y, Ma P, et al. Single-cell landscape reveals active cell subtypes and their interaction in the tumor microenvironment of gastric cancer. Theranostics. 2022;12:3818–33.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. Kumar MP, Du J, Lagoudas G, Jiao Y, Sawyer A, Drummond DC, Lauffenburger DA, Raue A. Analysis of single-cell RNA-seq identifies cell-cell communication associated with tumor characteristics. Cell Rep. 2018;25:1458-1468.e1454.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Li R, Yang X. De novo reconstruction of cell interaction landscapes from single-cell spatial transcriptome data with DeepLinc. Genome Biol. 2022;23:124.

    PubMed  PubMed Central  Google Scholar 

  21. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, Lee MJ, Wilk AJ, Darby C, Zager M, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184:3573-3587.e3529.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, Luo T, Xu L, Liao G, Yan M, et al. Cell Marker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res. 2019;47:D721-d728.

    CAS  PubMed  Google Scholar 

  23. Van de Sande B, Flerin C, Davie K, De Waegeneer M, Hulselmans G, Aibar S, Seurinck R, Saelens W, Cannoodt R, Rouchon Q, et al. A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nat Protoc. 2020;15:2247–76.

    PubMed  Google Scholar 

  24. Yang J, Zhao YL, Wu ZQ, Si YL, Meng YG, Fu XB, Mu YM, Han WD. The single-macro domain protein LRP16 is an essential cofactor of androgen receptor. Endocr Relat Cancer. 2009;16:139–53.

    CAS  PubMed  Google Scholar 

  25. Ruberto AA, Gréchez-Cassiau A, Guérin S, Martin L, Revel JS, Mehiri M, Subramaniam M, Delaunay F, Teboul M. KLF10 integrates circadian timing and sugar signaling to coordinate hepatic metabolism. Elife. 2021;10:e65574.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Guillaumond F, Gréchez-Cassiau A, Subramaniam M, Brangolo S, Peteri-Brünback B, Staels B, Fiévet C, Spelsberg TC, Delaunay F, Teboul M. Kruppel-like factor KLF10 is a link between the circadian clock and metabolism in liver. Mol Cell Biol. 2010;30:3059–70.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12:1088.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. Schweighoffer E, Tybulewicz VL. BAFF signaling in health and disease. Curr Opin Immunol. 2021;71:124–31.

    CAS  PubMed  Google Scholar 

  29. Gross JA, Dillon SR, Mudri S, Johnston J, Littau A, Roque R, Rixon M, Schou O, Foley KP, Haugen H, et al. TACI-Ig neutralizes molecules critical for B cell development and autoimmune disease. impaired B cell maturation in mice lacking BLyS. Immunity. 2001;15:289–302.

    CAS  PubMed  Google Scholar 

  30. Rauch M, Tussiwand R, Bosco N, Rolink AG. Crucial role for BAFF-BAFF-R signaling in the survival and maintenance of mature B cells. PLoS ONE. 2009;4:e5456.

    PubMed  PubMed Central  Google Scholar 

  31. Batlle E, Henderson JT, Beghtel H, van den Born MM, Sancho E, Huls G, Meeldijk J, Robertson J, van de Wetering M, Pawson T, Clevers H. Beta-catenin and TCF mediate cell positioning in the intestinal epithelium by controlling the expression of EphB/ephrinB. Cell. 2002;111:251–63.

    CAS  PubMed  Google Scholar 

  32. Park I, Lee HS. EphB/ephrinB signaling in cell adhesion and migration. Mol Cells. 2015;38:14–9.

    PubMed  Google Scholar 

  33. Chumley MJ, Catchpole T, Silvany RE, Kernie SG, Henkemeyer M. EphB receptors regulate stem/progenitor cell proliferation, migration, and polarity during hippocampal neurogenesis. J Neurosci. 2007;27:13481–90.

    CAS  PubMed  PubMed Central  Google Scholar 

  34. Genander M, Halford MM, Xu NJ, Eriksson M, Yu Z, Qiu Z, Martling A, Greicius G, Thakar S, Catchpole T, et al. Dissociation of EphB2 signaling pathways mediating progenitor cell proliferation and tumor suppression. Cell. 2009;139:679–92.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. Tjong WY, Lin HH. The RGD motif is involved in CD97/ADGRE5-promoted cell adhesion and viability of HT1080 cells. Sci Rep. 2019;9:1517.

    PubMed  PubMed Central  Google Scholar 

  36. Hsiao CC, Keysselt K, Chen HY, Sittig D, Hamann J, Lin HH, Aust G. The Adhesion GPCR CD97/ADGRE5 inhibits apoptosis. Int J Biochem Cell Biol. 2015;65:197–208.

    CAS  PubMed  Google Scholar 

  37. Rosa M, Noel T, Harris M, Ladds G. Emerging roles of adhesion G protein-coupled receptors. Biochem Soc Trans. 2021;49:1695–709.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. Chiang EY, de Almeida PE, de Almeida Nagata DE, Bowles KH, Du X, Chitre AS, Banta KL, Kwon Y, McKenzie B, Mittman S, et al. CD96 functions as a co-stimulatory receptor to enhance CD8(+) T cell activation and effector responses. Eur J Immunol. 2020;50:891–902.

    CAS  PubMed  Google Scholar 

  39. Siebert N, Xu W, Grambow E, Zechner D, Vollmar B. Erythropoietin improves skin wound healing and activates the TGF-β signaling pathway. Lab Invest. 2011;91:1753–65.

    CAS  PubMed  Google Scholar 

  40. Sorg H, Harder Y, Krueger C, Reimers K, Vogt PM. The nonhematopoietic effects of erythropoietin in skin regeneration and repair: from basic research to clinical use. Med Res Rev. 2013;33:637–64.

    CAS  PubMed  Google Scholar 

  41. Brines M, Cerami A. Erythropoietin-mediated tissue protection: reducing collateral damage from the primary injury response. J Intern Med. 2008;264:405–32.

    CAS  PubMed  Google Scholar 

  42. Moffitt JR, Lundberg E, Heyn H. The emerging landscape of spatial profiling technologies. Nat Rev Genet. 2022;23:741–59.

    CAS  PubMed  Google Scholar 

  43. Longo SK, Guo MG, Ji AL, Khavari PA. Integrating single-cell and spatial transcriptomics to elucidate intercellular tissue dynamics. Nat Rev Genet. 2021;22:627–44.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. Williams CG, Lee HJ, Asatsuma T, Vento-Tormo R, Haque A. An introduction to spatial transcriptomics for biomedical research. Genome Med. 2022;14:68.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. Yu G, Wang LG, Yan GR, He QY. DOSE: an R/Bioconductor package for disease ontology semantic and enrichment analysis. Bioinformatics. 2015;31:608–9.

    CAS  PubMed  Google Scholar 

Download references


We thank Dr. Juan Wang for the helpful discussion.

Code availability

The source code used for analyzing the data of this study can be accessed in GitHub:


This study was supported by grants from Program of the Science and Technology Commission of Shanghai Municipality "Science and Technology Innovation Action Plan" Medical Innovation Research Special Project (Grant/Award Number: 22Y11905700), program of Shanghai Shen Kang Clinical Research Incubation Plan(Grant/Award Number: SHDC12018X01) and program of Shanghai Municipal Health Commission (Grant/Award Number: 201940293).

Author information

Authors and Affiliations



F.T. and G.C. conceived and supervised the study. Y.X., L.L., Y.Y., C.Z., P.L., Y.W., M.C., and J.W. analyzed the data. F.T., G.C., and J.W. wrote the manuscript. G.C. and F.T. revised and finalized the manuscript. All authors have read and agreed to the published version of the manuscript.

Corresponding authors

Correspondence to Fei Tan, Jiling Wen or Geng Chen.

Ethics declarations

Ethics approval and consent to participate

The Ethics Committee of Shanghai Skin Disease Hospital has approved the study on March 1st, 2021 with IRB approval number of SSDH-IEC-SG-029–3.0. All methods were carried out in accordance with relevant guidelines and regulation. The parental assent of the children and informed consent were obtained from all individual participants or legally authorized representatives before entering into the study.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Supplementary Table S1. Top 10 markers in each cell cluster used for determining corresponding cell types.

Additional file 2.

Additional file 3.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Tan, F., Xuan, Y., Long, L. et al. Single-cell analysis of human prepuce reveals dynamic changes in gene regulation and cellular communications. BMC Genomics 24, 514 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: