Two mouse models of RDEB were utilized in this study. C57BL6/J Col7a1^+/− mice, kindly provided by Jouni Uitto, MD, at Jefferson Medical College, were developed by targeted ablation of the COL7A1 gene through out-of-frame deletion of exons 14–18 (31). Breeding of this mouse strain gives rise to Col7a1^−/− (KO) mice. KO mice have a median life span of 2 days without treatment and rarely survive past 3 weeks of age (32). Col7a1^fNeo mice, kindly provided by Alexander Nyström at the University of Heidelberg, were developed on amixed C57BL/6 129sv background by replacing an 11-kb genomic fragment spanning exon 2 of Col7a1 with a targeting construct containing a phosphoglycerate kinase promoter-driven neomycin phosphotransferase (PGK-Neo) expression cassette (7). Breeding of this mouse strain gives rise to C7^hypo mice. C7^hypo mice have an improved median life span (around 12 days) than the KO mice and develop mitten deformities in long-term survived mice. Both RDEB mouse models developed blisters at birth and were validated for their genotype by PCR. All animal studies were conducted using protocols approved by the New York Medical College Institutional Animal Care and Use Committee (IACUC).

To derive mouse fibroblasts, paw skin was dissected from WT and KO mice and digested in 1% dispase overnight at 4˚C. The dermis was then separated from the epidermis, cut into small pieces, and cultured in a fibroblast medium for the outgrowth of fibroblasts. Human RDEB-specific fibroblasts were derived from the nonlesional, noninflamed skin of six individuals with RDEB after informed, written consent. The demographics of RDEB patient donors are described in Table 1 (29, 33, 34). The normal control fibroblasts included foreskin fibroblasts (NC1) (33) and two adult dermal fibroblast lines (NC2, cat # 106-05A, Sigma-Aldrich and NC3, cat # PCS-201-012, ATCC). Fibroblasts were maintained in low-glucose Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% FBS. For cytokine stimulation, each cell line was seeded in duplicate wells of six- or 12-well culture plates (Corning, Lowell, MA, USA) and cultured for 24 h. The medium was then replaced with a low-glucose DMEM medium (basal medium) for 24 h. For the human fibroblasts, recombinant human IL-1α (Biolegend, San Diego, CA, USA), IL-6 (R&D, Minneapolis, MN, USA), TNF (Biolegend, San Diego, CA, USA), IFN-γ (Biolegend), or TGF-β1 (R&D, Minneapolis, MN, USA) was added to a final concentration of 4 ng/mL, 2 ng/mL, 100 ng/mL, 100 ng/mL, and 5 ng/mL to the basal coculture medium, respectively. For the mouse fibroblasts, recombinant mouse IL-1α (Biolegend), TNF Biolegend, San Diego, CA, USA), IFN-γ (Biolegend, San Diego, CA, USA), or TGF-β1 (Biolegend, San Diego, CA, USA) was added to a final concentration of 4 ng/mL, 100 ng/mL, 100 ng/mL, and 5 ng/mL to the basal coculture medium, respectively. Cytokine treatment concentrations were chosen either based on previous literature (29) or from the lowest concentration needed to elicit the same maximal response of PDPN expression among each fibroblast cell line independently. At the end of the 24-h coculture, the media was collected for ELISA quantification of sST2, and the cells were trypsinized with 200 µL of TrypLE and washed with Magnetic-activated cell sorting (MACS) rinsing buffer. Collected cells were immediately stained for flow cytometry analysis.

Each fibroblast sample was blocked in MACS staining buffer and subsequently stained with PDPN (PE) and ST2 (647)-conjugated antibodies (Biolegend) for 30 min. Samples were washed and resuspended with 300 µL MACS rinsing buffer and immediately run on BD FACS Canto II Flow Cytometer. Using FCS Express 7 Flow, fibroblasts were gated from an FSC and SSC plot and measured for PE and 647 MFI to determine PDPN and membrane-bound ST2 expression, respectively.

To quantitate the cytokines in the skin, paw skin lysate was prepared following lysis of stripped skin tissue in the presence of protease inhibitors (Cell Signaling Technology, Danvers, MA, USA) and homogenization in gentleMACS™ M tubes using the gentleMACS Dissociator (Miltenyi Biotec, Charlestown, MA, USA). Protein concentrations were determined using the Bio-Rad Protein Assay (Bio-Rad Laboratory, Hercules, CA, USA). Paw skin lysate samples were diluted in assay diluent from the appropriate assay kit to a consistent concentration across all samples, either to 10 µg or 100 µg per 50 µL or 100 µL, depending on cytokine abundance and assay type, using the total concentration measured by the Bio-Rad Protein Assay. Mouse IL-1α, IL-1β, and Th1/Th2/Th17 cytokines were quantitated using the CBA Flex set (BD Bioscience, Franklin Lakes, NJ, USA). Mouse IL-1ra and TGF-β1 were quantitated using an ELISA kit (R&D). Mouse IL-33, mouse ST2, and human ST2 were quantitated using an ELISA kit (Thermo Fisher Scientific, Wilmington, DE, USA). ELISA and CBA assays were performed according to the protocols provided by each manufacturer. Optical densities (OD) from ELISA assays were measured with the FilterMax F5 Multi-Mode Microplate Reader (Molecular Devices, San Jose, CA, USA). Cytokine concentration was determined by plotting the OD or median fluorescence intensity (MFI) of each sample on each cytokine concentration curve and multiplying the result by a dilution factor to give a sample concentration in picograms per milligram. Plasma was collected from blood in EDTA-containing tubes and diluted 1:2 in the assay diluent from the appropriate assay kit. Plasma cytokine concentration was measured the same as the lysate but expressed in picograms per milliliter of concentrations. Supernatant samples from cell culture were also diluted 1:2 in assay diluent and had cytokines quantitated by ELISA. Boxplots, dot plots, and trendlines were generated with the ggplot2 package in R.

The collected tissues were embedded in Tissue-Tec OCT Compound and snap-frozen on dry ice. For each specimen, 6 µm serial sections were cut and stored at −80°C. Sections were fixed in freshly diluted 4% paraformaldehyde (Electron Microscopy Sciences, Hatfield, PA, USA) and blocked with M.O.M. blocking reagent (for antibodies raised in mice) or CAS-Block (for antibodies raised in other species of animals) with 0.1% Triton (Sigma, St. Louis, MO, USA). The slides were then incubated with primary and corresponding secondary antibodies and mounted in a Vectashield mounting medium containing DAPI (Vector Laboratories, Newark, CA, USA). Images were acquired using the EVOS M5000 imaging system (Thermo Fisher Scientific) or ZEISS LSM 980 plus confocal microscope using the same settings between the different groups in each set of experiments.

The skin was dissected from the paws of D11 WT and KO mice (N = 2), cut into small pieces, and dissociated into a single-cell suspension using a combination of enzymatic digestion and mechanical dissociation. Briefly, tissues were digested using the components in Multi Tissue Dissociation Kit 1 (Miltenyi Biotec), following the manufacturer’s recommendation for 5 h at 37°C. Samples were then transferred into gentleMACS™ C tubes and loaded onto the gentleMACS™ Octo Dissociator with Heaters (Miltenyi Biotec) with the program selected for skin dissociation. Cells were then passed through 70-μm and 40-μm filters (Falcon, Corning, NY, USA) sequentially. Dead cells were removed using Ficoll-Plaque PLUS (GE Healthcare, Chicago, IL, USA) separation. The viability of the cells was evaluated by Trypan blue staining. Library preparation and sequencing were done by Singulomics Corporation (https://singulomics.com/, Bronx, NY, USA). Viable cell suspensions were loaded into the Chromium Controller (10x Genomics, Pleasanton, CA, USA) to generate gel beads-in-emulsion (GEM) with each GEM containing a single cell as well as barcoded oligonucleotides. Next, the GEMs were placed in the SimpliAmp 96-well Thermal Cycler (Thermo Fisher Scientific, Wilmington, DE, USA), and reverse transcription was performed in each GEM (GEM-RT). After the reaction, the complementary cDNA was amplified and cleaned using Silane DynaBeads (10X Genomics, Pleasanton, CA, USA) and the SPRIselect Reagent kit (Beckman Coulter, Indianapolis, IN, USA). Amplified full-length cDNAs from poly-adenylated mRNA were then used to generate 5′ Gene Expression library (Chromium Next GEM Single Cell 5′ Reagent Kits v2 (Dual Index)), following the manufacturer’s instructions (10x Genomics, Pleasanton, CA, USA). Amplified cDNAs and the libraries were measured by the Qubit dsDNA HS assay (Thermo Fisher Scientific, Wilmington, DE, USA), and quality was assessed by the BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA). Libraries were sequenced on a NovaSeq 6000 instrument (Illumina, San Diego, CA, USA), and reads were subsequently processed using the 10x Genomics Cell Ranger analytical pipeline (v6.0.1) and mouse reference genome mm10. Dataset aggregation was performed using the cellranger aggr function, normalizing for the total number of confidently mapped reads across libraries. A total of 14,856 cells were sequenced, with 24,774 postnormalization mean reads per cell. Methods for bioinformatic analysis, including basic processing, cell assignment algorithms, analysis of subtype cells, GO term scoring, and trajectory inferences are described in supplementary data.

Globally differentially expressed genes with a log2-fold change greater than 1.2 and an adjusted p-value of less than 0.05 were uploaded into Qiagen’s ingenuity pathway analysis (IPA) system for core analysis.

The Chi-square test was used to determine the statistical significance of the quantities of each cell type between the WT and KO samples of the scRNAseq data. Statistical significance between groupings of WT vs. RDEB phenotypes in protein concentration and MFI was calculated by one-tailed Student’s t-tests, and significance between age categories and treatments of RDEB samples was calculated by ANOVA with Tukey’s correction.

To focus on an early phase of RDEB disease progression, we isolated single cells from the front paw skin of 11-day-old KO mice (n = 2) and WT controls (n = 2) for 10x Genomics scRNAseq. We confirmed that both KO mice exhibited typical histological features of RDEB skin, including separation of epidermis and dermis, dermal thickness, infiltration of immune cells, and hyperproliferation of epidermal cells, by sectioning and imaging their hind paws (Supplementary Figure S1). Following dimensionality reduction and uniform manifold approximation and projection (UMAP), cells of similar transcriptional profiles were clustered into 11 cell types identified as keratinocytes, fibroblasts, chondrogenic fibroblasts, vascular endothelial cells, lymphatic endothelial cells, perivascular cells, Schwann cells, mast cells, lymphocytes, neutrophils, and antigen-presenting cells (Figure 1A). The signature genes defining individual cell types are presented in Figure 1B and Data File S1. The KO samples contained significantly higher numbers of fibroblasts, perivascular cells, and immune cells, with lower numbers of keratinocytes and Schwann cells than the WT samples (Figure 1C).

Differential gene expression between KO and WT samples revealed global as well as cell-type-specific changes (Data Files S2, S3). Genes classified as hallmark signatures of cancer (35), such as those related to glycolysis, hypoxia-induced autophagy, stress response, fatty acid binding and lipid droplet accumulation, amino acid transport and metabolism, angiogenesis, and proinflammatory cytokines such as Il1a, Il1b, C3, and Cxcl1, were upregulated in RDEB samples across all cell types (Figure 1D).

Furthermore, several genes whose expression predicts the poor prognosis of many cancer types and represents potential anticancer targets, i.e., Lgals3, Smox, Ndrg1, Neat1, and Rbm39, were significantly upregulated in most of the KO cell types (Data File S2). Overall, these data confirmed the existence of an early proinflammatory environment in RDEB mouse skin. We included a table to highlight notable genes differentially expressed in different cell types (Table 2).

Unsupervised clustering within the fibroblast population identified nine subclusters represented by their most differentially expressed gene signatures (Supplementary Figure S2, Data File S4). Figure 2A highlights clusters that either were enriched (subclusters 0, 1, and 2), exhibited different trajectories (subcluster 3), or were significantly less abundant in the KO sample (subclusters 4 and 5). Using gene ontology (GO) enrichment analysis, we identified that subclusters 0, 1, 2, and 3 upregulated genes involved in stress response. More specifically, subcluster 0 and parts of subclusters 1 and 2 expressed genes involved in hypoxia-induced signaling and glycolysis, whereas subclusters 0 and 3 expressed genes related to cytokine and chemokine signaling such as Ccl2, Cxcl1, and Csf3 (Figure 2B).

While our samples are all from a single time point, the distribution of the subclusters suggests differentiation patterns that could be further elucidated through trajectory analysis. The resulting analysis suggested that subcluster 0 along with subclusters 1 and 2 represent two dynamic processes of inflammatory fibroblast activation and myofibroblast differentiation, respectively (Figures 2C, D). Fibroblast subcluster 0, and to a lesser extent subcluster 3, are enriched in gene encoding pattern recognition receptors that sense danger signals and trigger innate immune response, including Tlr2 and Tlr4 (Toll-like receptors), Nod2, Ifih1 (MDA5), and Ddx58 (RIG-1) (Figure 2E). These two subclusters were also enriched with the expression of receptors to inflammatory cytokines, i.e., IL-1 receptor (Il1r1), TNF receptor 2 (Tnfrsf1b), and IL-6 receptor (Il6ra) (Figure 2F). Meanwhile, the IFN-γ receptor (Ifngr) and TGF-β receptor 3 (Tgfbr3) were more ubiquitously expressed by all the subclusters of WT and KO fibroblasts (Figure 2F). Interestingly, upon transitioning from subcluster 3 into subcluster 0, the expression of Mmp2, Mmp3, Mmp13, Mmp14, and Ptgs2 was significantly upregulated, and the expression of Timp3 was downregulated (Figure 2D), all of which is consistent with the induction of IL-1 signaling (36, 37). Moreover, besides their Pdpn and Il1rl1 signature gene expression, subcluster 0 expressed high levels of C3, Saa3, Il11, and Lif (Figure 2G), which is reminiscent of the phenotype of IL-1-activated iCAFs. While Tgfb1 was expressed, the myofibroblast markers, like Acta2 and Tagln, were largely absent in subcluster 0 (Figure 2H). In addition, this subcluster is enriched with the angiogenesis and lymphangiogenesis markers Vegfa, Vegfb, and Pgf (Figure 2I). We postulate that subcluster trajectory 4–3–0 could be interpreted as a path through which RDEB fibroblasts acquire immune effector functions through activation by inflammatory cytokines such as IL-1 (Figure 2C); therefore, we categorized subcluster 0 as inflammatory fibroblasts. In contrast, enhanced expression of Tagln and Serpine2 and downregulation of Mmp2 in subclusters 1 and 2 suggest that subcluster trajectory 4–3–1–2 may represent a transition toward a myofibroblast phenotype (Figures 2C, D) (38–40). Specifically, upregulation of the stress response genes (Figure 2B) and low expression of Ptgs2 (COX2) and Acta2 (αSMA) (Figure 2D) were consistent with the phenotype of subclusters 2 and 1 being proto-myofibroblasts, an intermediate stage of fibroblasts that have responded to mechanical tension and maintain the capacity to differentiate into myofibroblasts upon further stimulation (41). Supportively, a reverse correlation of Tnc and Ccn2 expression was observed between two subcluster trajectories: subcluster trajectory 4–3–0 simultaneously upregulated Tnc and downregulated Ccn2, while subcluster trajectory 4–3–1–2 appeared to keep the two genes expressed more consistently (Figure 2D). This observation is consistent with the reported opposite effects of IL-1α and TGF-β on Ccn2 and Tnc expression (42–44).

scRNAseq analysis also demonstrated that the cells within the KO inflammatory fibroblast subset were mostly positive for Prdm1 (encoding B-lymphocyte-induced maturation protein 1 (BLIMP1)) and Lrig1, and negative for Dpp4 (CD26) (Figure 2J), which suggests that their lineage stems from papillary dermal fibroblasts based on published lineage tracing studies (45). Interestingly, Sca1(ly6a)⁺Dlk1^+/− cells, which characterize the fibroblasts in the lower dermis (hypodermal fibroblasts) (45), were also identified as part of the KO inflammatory fibroblast subset (Figure 2J). The recruitment of Sca1⁺ fibroblasts to skin wounds has been demonstrated to be among the initial waves of dermal repair (45). Moreover, Scar1⁺ oral fibroblasts exhibited immunomodulatory functions through the secretion of CCL2 (46). Therefore, the detection of Ly6a expression within the KO inflammatory fibroblast subset could represent the activation and recruitment of fibroblasts from the lower to upper dermis in response to dermal–epidermal separations in RDEB skin. Supportively, immunohistochemistry analysis demonstrated that the inflammatory fibroblasts in KO skin, represented by the expression of IL-1R1 and vimentin (VIM), were mainly identified at the dermal–epidermal junction where the blisters occur (Figure 2K).

Unsupervised clustering among immune cells identified 12 immune cell types/subtypes (Figure 3A; Data File S5). Clusters 0 and 1 were identified as αβ and γδ T cells, based on their expression of Trac and Trbc and Tcrg-c1 and Trdv4, respectively. Clusters 2 and 3 expressed neutrophil markers S100a8, S100a9, and Ngp and thus were identified as neutrophils, of which cluster 2 represented an activated neutrophil subpopulation unique in the KO sample. Similarly, clusters 4 and 5 together were identified as Langerhans cells based on their expression of Cd74 and Cd207, with cluster 4 representing an activated Langerhans cell subpopulation. The number of cells within clusters 6, 7, and 8 was significantly higher in the KO samples than in the WT samples, and these clusters were largely positive for the expression of H2-Ab1, Adgre1, Cd68, Fcrg1, and Itgam. For the scope of this manuscript, we grouped these clusters as macrophages/dendritic cells (Mφ/DC).

The αβ T cells contained conventional Cd4⁺Foxp3⁻ T helper cells, Cd8⁺ cytotoxic T cells, and Foxp3⁺ T regulatory (Treg) cells (Figure 3B). While Cd8 transcript levels were low in our data analysis, immunohistochemical (IHC) analysis demonstrated an abundant amount of CD8⁺ cells in the dermis of KO paw skin (Figure 3C). Meanwhile, CD8⁺ cells were absent in the WT paw skin (Figure 3C), which was consistent with the reported immunological features of murine paw skin (47). Therefore, we postulated that the low Cd8 detection in our KO dataset could be due to the limited sequencing depth of the scRNAseq technique not fully representing the transcriptional complexity of T cells. It may also be due to an incoherence between mRNA and protein expression (48). Nevertheless, our data demonstrated that the KO αβ T cells exhibited upregulated expression of both inhibitory (e.g., Pdcd1/PD-1, Ctla4) and costimulatory (e.g., Tnfrsf9/4-1-BB, Tnfrsf18/GITR, and Tnfrsf4/OX40) genes (Figure 3B), a phenomenon that has been documented in T cells within the tumor microenvironment (12, 49). Here, IHC analysis demonstrated that CD8⁺ T cells were PD-1⁺ when recruited to the KO mouse skin a week after birth (Figure 3C), suggesting that these cells were effector cells of adaptive immunity. The number of T cells and their expression of PD-1 remained elevated in the KO skin for the following 2 weeks (Figure 3C). Persistent PD-1 expression was also observed in T cells in the skin of C7^hypo mice that survived long-term and exhibited mitten deformity (Supplementary Figure S3A) consistent with the previous study (30) and suggestive of a process of T-cell exhaustion in the chronic inflammatory environment of RDEB skin. Interestingly, scRNAseq analysis revealed that Pdcd1⁺ T cells in KO samples can be separated into Tcf7-positive and negative subpopulations (Figure 3B). Tcf7 encodes T-cell factor-1 (TCF-1) that marks either naïve or central memory stem-like precursor CD8⁺ T cells with self-renewal capacity (50). Importantly, Pdcd1^highTcf7⁺ T cells with low expression of coinhibitory receptors such as Havcr2 (TIM-3), Tigit, and Lag3 in the KO sample (Figure 3B) were consistent with the phenotypes of progenitor-exhausted T cells that are long-lived, have stem cell-like properties, and respond to immune checkpoint blockade therapy (51). In contrast, the expression of Tox, which plays a crucial role in the generation and maintenance of an exhausted T-cell phenotype (52), and other coinhibitory receptors Havcr2, Tigit, and granzyme B (Gzmb), were more associated with the Pdcd1⁺Tcf7⁻ T cells in the KO samples (Figure 3B).

Correlated with increased expression of Pdcd1 in KO T cells, elevated expression of Cd274 (PD-L1) and/or Pdcd1lg2 (PD-L2) was observed in Langerhans cells, neutrophils, and macrophages/dendritic cells (Figures 3D, E). Interestingly, although CD68⁺ cells were identified in the dermis of newborn KO paw skin, they had no expression of PD-L1 (Figure 3D). In accordance with the timing of PD-1⁺ T-cell recruitment, PD-L1 was identified in CD68⁺ cells in the paw skin of 7-day-old and older KO mice, suggestive of the activation of the PD-1/PD-L1 axis in RDEB immune cells (Figure 3D). PD-L1 expression was also observed in the CD68⁺ cells in C7^hypo mice with mitten deformity (Supplementary Figure S3B) but was not identified in the WT CD68⁺ cells at any time points (represented by the d11-timepoint in Figure 3D). In addition, although the expression of costimulatory molecules such as Tnfsf9 (4-1BBL), Tnfsf4 (OX40L), and CD86 was increased in some KO immune cell types, the extent of their upregulation was overall less than the elevation of immune suppressive molecules such as Lgals9, Lgals3, Mif, CD74, Adora2a, Adora2b, Arg1, Arg2, and Il4i1 (Figure 3E).

Significant upregulation of Il1a (IL-1α), Il1b (IL-1β), and their antagonist Il1rn (IL-1ra) was identified in multiple immune cell types in the KO sample (Figure 3E), further supporting the subcluster trajectory of inflammatory fibroblast activation. Moreover, genes encoding inflammasome sensor proteins (Nlrp3) and proteolytic enzymes (Casp1, Casp8, and Gzmb) were significantly upregulated in multiple immune cell types in the KO sample. These enzymes are essential not only for processing pro-IL-1β into an active form but also for pro-IL-1α to further increase its bioactivity, even though IL-1α is active in its full-length form (53). Upregulation of other cytokines, such as Tnf and Il6, was also observed in multiple immune cell types in the KO sample (Figure 3E). Moreover, chemoattractant proteins were highly expressed by Langerhans cells, macrophages/dendritic cells, neutrophils, and B cells in the KO sample, while chemokine receptors were upregulated in T cells and B cells (Figure 3E), suggesting their recruitment by these chemotactic axes.

A heatmap of single cells within the Langerhans cell population of WT and KO samples demonstrated that the novel cluster 4 in the KO sample exhibited significant upregulation of antiapoptosis genes, such as Bcl2a1a, and genes involved in actin filament binding, lamellipodium assembly, podosome formation, and cell projection and adhesion, such as Abracl, Wdr1, and Fscn1 (Figure 3F). Together with the significant upregulation of Ccr7, induced only at later stages of Langerhans cell maturation and endowing migratory capacity (54), our results suggested that this novel cluster could represent mobilized Langerhans cells migrating to lymph nodes in RDEB mice. The migration of Langerhans cells is an essential mechanism for the induction of Tregs and can be induced by IL-1 or TNF (55).

The heatmap of single cells within clusters 2 and 3 also revealed neutrophil activation in KO samples, exhibited by the upregulation of genes involved in neutrophil chemotaxis and chemokine-mediated signaling, IL1 signaling, TNF signaling, NF-кB signaling, C-type lectin receptor signaling, Toll- and NOD-like receptor signaling, integral components of membrane composition, as well as antiapoptotic regulation gene activation (Figure 3G). Moreover, this activation appears to have occurred with a gradual, three-stage progression. One part of cluster 3 contains lower gene expression and is present in both WT and KO samples. The other part of cluster 3 contains a slightly elevated gene expression and is uniquely expressed in the KO sample. In contrast, cluster 2 is almost exclusively found in the KO sample and has the highest expression of these genes (Figure 3G).

Overall, these data suggested pervasive recruitment and activation of immune cells in RDEB skin that supported the induction of inflammatory fibroblast activation trajectory in the dermal microenvironment.

Since skin compartments are tightly interconnected, we next explored signaling at the epidermal level. Differential expression analysis between WT and KO keratinocytes revealed that multiple keratin genes, Krt2, Krt5, Krt10, Krt15, and Krt24, were significantly downregulated, and Krt6a, Krt6b, Krt16, Krt17, Sprr1a, and Sprr1b were significantly upregulated in KO keratinocytes as compared to the WT (Supplementary Tables S2, S3). Subclustering and UMAP analysis of keratinocytes (Figures 4A–C) demonstrated that most of the differential expression occurred within interfollicular epidermal basal cells (clusters 7 and 8), suprabasal cells (clusters 9 and 10), and germinative layer cells (cluster 5), based on the expression patterns of each subcluster (Data File S6) (56). The imbalanced expression of Krt and Sprr genes in RDEB mouse keratinocytes was consistent with a previous transcriptome analysis of RDEB patients’ skin (57). Importantly, the downregulation of Krt2/15 and upregulation of Krt6/16/17 have been associated with inflammation and the proliferative stage of keratinocyte differentiation under pathological conditions such as psoriasis (58–60). Supportively, strong KRT16 immunostaining was observed in the suprabasal epidermal layer of KO skin, particularly in the area with hyperkeratosis (Figure 4D). The KO keratinocytes also had upregulated expression of Sprr1a and Sprr1b (Figure 4C), suggesting their role in protecting against the damaged skin barrier in RDEB. However, Sprr1b has also been demonstrated to be overexpressed in multiple types of cancer with poor prognosis (61, 62).

Importantly, the epidermal basal cells in cluster 8 of KO samples were also enriched with the expression of inflammatory cytokines and chemokines, including Tslp, Il1f9 (IL-36γ), and Cxcl16 (Figure 4E), all of which have been demonstrated to contribute to the pathogenesis of psoriasis (63–65). Thymic stromal lymphopoietin (TSLP) is well known for its association with type 2 immunity at the barrier surface, allergic conditions, chronic inflammatory diseases, and cancer (66). It has been demonstrated that following the binding of TSLP to its receptor TSLPR on immune cells, the TSLP-TSLPR complex recruits IL-7R with high affinity to activate TSLP signaling (67). This event thus sequesters IL-7R from IL-7R signaling and may further influence IL-2R signaling, as IL-7R and IL-2R share a common γ-chain for their respective activation (68). Therefore, the upregulation of Tslp may be one of the mechanisms underlying the crosstalk between keratinocytes and immune cells in the RDEB skin. Importantly, not only was Tslp upregulated in keratinocytes, but Il7r was also significantly upregulated in both T cells and macrophages/dendritic cells of the KO sample (Figures 4B, E). Moreover, recent studies demonstrated that following inflammatory stimulation, monocytes and macrophages were more likely to produce soluble IL-7R (sIL7R), a decoy for IL-7R signaling, than full-length IL-7R, which is normally produced by T cells (69). Therefore, there appeared to be multiple levels of immunomodulation in the KO sample, which warrants further investigation.

In line with specific immune regulation within the major cell types described above, the ingenuity pathway core analysis on globally differential gene expression demonstrated that immune responses, centered on IL-1 and regulated by other cytokines, constitute a major part of the upregulated signaling network in RDEB samples (Supplementary Figure S4). Therefore, we next profiled key cytokines in the skin lysate of KO mice and WT controls. Specimens from C7^hypo mice were included with KO as an additional validation of our findings and to increase the robustness of our data analysis.

Our results suggest that IL-1α, IL-10, TNF, IL-6, and IFN-γ were significantly elevated overall in the RDEB skin lysate compared the WT control (Figure 5A). The mean concentration of IL-1ra was high, but it was not significantly different between RDEB and WT mice. In addition, although scRNAseq showed a significantly higher expression of Il1b in RDEB samples, quantitation of IL-1β through two independent methods revealed similarly low levels in both groups.

Importantly, dynamic profiling of individual cytokines suggests that there are likely multiple waves of immune responses in RDEB skin from birth to early adulthood (Figure 5B; Supplementary Figure S7). At birth, RDEB skin contained significantly higher levels of both proinflammatory (IL-1α, IL-1β, TNF, IFN-γ, and IL-17A) and anti-inflammatory (IL-10) cytokines than the WT skin (Figure 5B). Among these, IL-1α exhibited the highest magnitude of elevation, with its concentration > 2.8-fold higher in the RDEB than WT skin. Meanwhile, the level of IL-33 was significantly lower in the RDEB skin than in the WT skin at birth. At postnatal days 1–5, IL-1α was still significantly elevated in the RDEB samples, while the other cytokines had lost the significant difference between the RDEB and WT samples. Interestingly, at around a week of age (d6–d10), the level of IL-1α came down in RDEB skin. The decreased level of IL-1α was accompanied by a significant decrease in IL-1ra and a significant increase in IL-6. The level of TNF was also much higher in the RDEB skin than in the WT, although it was not statistically significant due to the high variability caused by spikes in TNF concentrations in some KO samples. Importantly, from postnatal day 11 to day 20, the level of IL-1α was again significantly upregulated in RDEB, not only compared to the WT but also to the RDEB lysate at earlier time points (Figure 5B; Supplementary Figure S5). Moreover, the concentrations of IFN-γ, IL-6, TNF, and TGF-β1 were all elevated. In contrast, the level of IL-33 significantly decreased. The levels of all the proinflammatory cytokines appeared to diminish when the RDEB mice reached young adulthood; however, the concentrations of IL-1α, TNF, and IFN-γ were still significantly higher in the RDEB than in the WT skin.

The levels of IL-1α in the plasma also reflected the dynamics of IL-1α in the RDEB skin lysate: they were not significantly different between the RDEB and WT mice at around 1-week of age; however, were significantly elevated at later days (Figure 5C). Of note, there was no significant difference between the KO and C7^hypo mice in the dynamics of all the measured cytokines except IL-1α (Supplementary Figure S6A). While the C7^hypo skin still had significantly higher levels of IL-1α than the WT skin, its level was significantly lower than the IL-1α in KO skin, suggesting that the complete lack of C7 resulted in a higher magnitude of IL-1α production (Supplementary Figure S6A). Importantly, we also detected suppression of tumorigenicity 2 (ST2) in the plasma of RDEB mice, resembling the dynamics of IL-1α (Figure 5D). ST2 is a member of the IL-1-receptor superfamily encoded by Il1rl1, a gene that was upregulated predominantly in the inflammatory fibroblast subset (cluster 0) in the KO samples (Figure 2). There are two forms of ST2, a membrane-bound form (STL) expressed primarily in hematopoietic cells and mediating signaling upon IL-33 binding and a soluble form (sST2) largely inducible in fibroblasts, serving as a decoy receptor for IL-33/ST2L and a serological biomarker for fibrosis (70). The observed decrease of IL-33 in RDEB skin lysate and increase of sST2 in RDEB plasma were thus suggestive of a downregulation of the IL-33/ST2L signaling axis in RDEB. However, the full extent to which sST2 might contribute to the pathological condition of RDEB remains to be determined.

As expected from the function of IL-1α and/or TNF on inducing NF-κB signaling, phosphorylated Iκ-Bα (pIκBα) was identified in the newborn RDEB skin, suggesting that NF-κB signaling is initiated before or at birth (Figure 5E). Upregulation of this signaling became apparent at a week of age, was further enhanced at 11 days of age, determined by the significant elevation in the levels of NF-κB subunits, p50/p105 and C-Rel and nuclear localization of p50, and appeared to be lower in the following timepoints (Figures 5E; Supplementary Figure S6B). Interestingly, phosphorylated STAT3 (pSTAT3) was absent in the newborn KO skin, sporadically detected at a week of age, became significantly positive at 11 days, and stayed high at 3 weeks of age (Figure 5F; Supplementary Figure S6C), suggesting that JAK/STAT signaling was activated later than NF-κB signaling. Moreover, p50/p105 and pSTAT3 were positive in the skin of 3-month-old C7^hypo mice (Supplementary Figure S6D), suggesting the persistent activation of both signaling pathways. In contrast, all these factors were at a low level or undetectable in the WT skin (Figures 5E, F). The results were congruent with the later peak of IL-6 and IFN-γ, two potent activators of JAK/STAT signaling in the RDEB skin, and with the reported IL-1-induced signaling cascade that leads to JAK/STAT activation in iCAFs (15).

We next evaluated the dynamics of inflammatory fibroblasts in RDEB mouse skin based on the expression of PDPN. PDPN is a mucin-type transmembrane protein that is normally expressed on lymphatic endothelial cells and fibroblastic reticular cells of lymphoid organs. PDPN is also utilized as a biomarker for CAFs, and its overexpression in CAFs, epithelial tumor cells, and inflammatory macrophages is associated with worse outcomes in many cancers (71). In WT mouse skin, PDPN expression was mainly associated with the vascular cells (LYVE⁺) (Figure 6A), whereas in RDEB mouse skin, in addition to its positive staining on the vascular cells, PDPN was associated with both fibroblasts (VIM⁺) and CD68⁺ immune cells (Figures 6A–C). Moreover, in line with the dynamics of JAK/STAT signaling, PDPN expression on fibroblasts was rare in the skin of RDEB mice less than a week of age and was substantially increased in the skin of older RDEB mice (Figure 6C). Importantly, PDPN⁺ fibroblasts were distinct from the αSMA⁺ myofibroblasts (Figure 6D), consistent with the scRNAseq analysis that suggested distinct trajectories of inflammatory fibroblast and myofibroblast differentiation in RDEB mouse skin.

It has been well established that inflammatory cytokines such as TNF, IL-1, and IL-17 are potent stimulators of fibroblast activation in inflammatory diseases (19, 72). We, therefore, asked whether stimulation of fibroblasts with the cytokines that were prominently elevated in RDEB mouse skin lysates would lead to upregulations of inflammatory fibroblast signature genes identified in our scRNAseq analysis, i.e., expression of PDPN and secretion of ST2, and whether there is any differential response between RDEB-derived and normal fibroblasts.

Here, we isolated fibroblasts from RDEB mice and WT controls at a week of age and analyzed their responses to mouse IL-1α, TNF, IFN-γ, and TGF-β individually following 24-h starvation in a serum-free (basal) medium. We demonstrated that PDPN expression on the cell surface and ST2 release in the conditioned medium (CM) were detected at similar levels in both RDEB and WT mouse fibroblasts under starved conditions and could be significantly enhanced under various stimulus conditions, confirming that they are fibroblast inflammatory responders (Figures 6E, F). TNF, IL-1α, and IFN-γ treatments all individually significantly enhanced PDPN expression and ST2 secretion in both RDEB and WT mouse fibroblasts, while TGF-β1 had no stimulatory effects on their expression (Figures 6E, F). Conversely, the addition of TGF-β1, but not IL-1α or IL-1α/IFN-γ, increased the fibroblast expression of αSMA, consistent with the effect of TGF-β1 on myofibroblast differentiation (Figure 6G). There was a trend of higher PDPN expression in RDEB fibroblasts compared to WT mouse fibroblasts under stimulatory conditions; however, it was not statistically significant (Figure 6E). ST2 was secreted at a significantly higher amount by RDEB mouse fibroblasts than WT following stimulation with TNF, IL-1α, and IFN-γ, respectively (Figure 6F).

We also analyzed RDEB human patient-derived fibroblasts (n = 6; Table 1) and normal control fibroblasts (n = 3) under basal conditions and following stimulation with recombinant human IL-1α, TNF, IFN-γ, IL-6, or TGF-β (Figure 7). A notable discrepancy between human and mouse fibroblasts is that under basal conditions, the expression of PDPN and ST2 was almost undetectable in human fibroblasts, whereas mouse fibroblasts expressed a detectable amount of both under basal conditions. However, consistent with the postulated trajectory of inflammatory fibroblast stimulation, they are significantly upregulated in response to inflammatory cytokine stimulation, similarly in both mouse and human fibroblasts. Moreover, TNF, IL-1α, or IFN-γ enhanced PDPN expression to a significantly higher level in RDEB patient-derived fibroblasts compared to normal controls (Figures 7A, C; Supplementary Figure S5). As for ST2, normal control fibroblasts did not produce ST2 under any condition (Figure 7B). While IL-1α stimulated certain RDEB patient-derived fibroblasts to secrete ST2, IFN-γ was the most potent ST2 stimulator among all RDEB patient-derived fibroblasts (Figure 7B). Membrane-bound ST2 (STL) was not detectable within all the stimulatory conditions (Supplementary Figure S7B), confirming the production of sST2 from the IL1RL1 gene in RDEB patient-derived fibroblasts. In addition, immunocytochemistry staining demonstrated that expression of the receptor to IL-1α, i.e., IL-1R1, can be induced specifically by IL-1α, and the level of its induction was more robust in RDEB patient-derived fibroblasts than normal controls (Figure 7C).