Integrated spatial multiomics reveals fibroblast fate during tissue repair

Abstract

In the skin, tissue injury results in fibrosis in the form of scars composed of dense extracellular matrix deposited by fibroblasts. The therapeutic goal of regenerative wound healing has remained elusive, in part because principles of fibroblast programming and adaptive response to injury remain incompletely understood. Here, we present a multimodal -omics platform for the comprehensive study of cell populations in complex tissue, which has allowed us to characterize the cells involved in wound healing across both time and space. We employ a stented wound model that recapitulates human tissue repair kinetics and multiple Rainbow transgenic lines to precisely track fibroblast fate during the physiologic response to skin injury. Through integrated analysis of single cell chromatin landscapes and gene expression states, coupled with spatial transcriptomic profiling, we are able to impute fibroblast epigenomes with temporospatial resolution. This has allowed us to reveal potential mechanisms controlling fibroblast fate during migration, proliferation, and differentiation following skin injury, and thereby reexamine the canonical phases of wound healing. These findings have broad implications for the study of tissue repair in complex organ systems.

Keywords: chromatin accessibility; fibrosis; multiomics; spatial epigenomics; spatial transcriptomics.

Fig. 1.

Wounding triggers polyclonal proliferation of tissue-resident fibroblasts. (A) Schematic of the Rainbow mouse construct. (B) Schematic showing wound healing model using Rainbow mice with local Cre recombinase induction using 4-hydroxytamoxifen liposomes (LiTMX). (C) Schematic showing a Rainbow wound cross-section. Black dotted line highlights wound scar area; arrows indicate the direction of cellular proliferation during wound healing. Structures are as labeled. (D) Representative confocal image of POD 14 wound cross-sections from Actin-Cre^ERT2::Rosa26^VT2/GK3 mice induced locally with LiTMX at the time of wound creation. Thick white dotted lines highlight scar boundaries. Individual Rainbow cell clones are highlighted with thin colored dotted lines. Arrows indicate direction of wound healing. n > 5. (Scale bar, 50 μm.) (E) Representative confocal images of unwounded skin from Actin-Cre^ERT2::Rosa26^VT2/GK3 mice induced locally with LiTMX. Thick white dotted lines highlight dermal boundaries. Individual Rainbow cell clones are highlighted with thin white dotted lines. n > 5. (Scale bar, 50 μm.) (F) Rainbow clone counts in wounds versus uninjured skin. n = 5 per condition. *P < 0.05. (G) Schematic of dorsal, stented, excisional wound healing in the Rainbow mouse model (whole-mount view), with polyclonal proliferation of Rainbow fibroblasts from the outer wound edge inward across time from POD 2 (Left), to POD 7 (Middle), to POD 14 (Right). Black arrows highlight the apparent direction of proliferation. (H) Representative confocal imaging of a POD 14 whole-mounted wound harvested from Actin-Cre^ERT2::Rosa26^VT2/GK3 mice showing the polyclonal proliferation of wound fibroblasts radially toward the center of the wound (dark area at center). White arrows highlight the direction of cell proliferation; individual cell clones are highlighted with thin colored dotted lines. Bottom subpanels denote individual Rainbow color contributions to merged image. mCh, membrane (m)Cherry; mOr, mOrange; mCe, mCerulean; eG, eGFP. n > 5. (I) Schematics illustrating microdissection strategy for isolation of inner and outer wound regions (Top), followed by enzymatic separation of the dermal scar from the epi- and hypodermis (Bottom). (J) Heatmap displaying expression data for genes significantly different between POD 7 inner and outer region wound fibroblasts. Legend at Right displays fold change. (K) Gene Ontology (GO) enrichment analysis comparing gene expression data from POD 7 inner and outer region wound fibroblasts. Top shows GO biological processes up-regulated in inner region fibroblasts compared with outer region fibroblasts, while the Bottom shows the same for outer region fibroblasts compared with inner. Top 10 most significant gene sets are displayed for each condition.

Fig. 2.

Single-cell transcriptomic and chromatin accessibility analyses delineate mechanoresponsive fibroblast subpopulations. (A) Schematic illustrating single-cell (sc) isolation of Rainbow wound fibroblasts from inner and outer wound regions (highlighted with black dotted lines). For scRNA-seq, mCerulean⁺ fibroblasts were arbitrarily selected from the available Rainbow colors and used for the remaining experiments in this figure. (B) (Left) Uniform manifold approximation and projection (UMAP) embedding showing scRNA-seq data from mouse wound fibroblasts FACS isolated using a lineage-negative sort strategy (29) from POD 2, POD 7, and POD 14, digitally pooled and clustered in a manner agnostic to POD and inner versus outer wound regions. Four unique fibroblast clusters were identified (clusters 1 through 4). Dotted lines highlight individual cluster distributions. (Right) Recoloring of Left UMAP plot based on fibroblast tissue region: inner (black) versus outer (orange). (C) CytoTRACE analysis of scRNA-seq data using the UMAP embedding from F. Shading indicates inner (light gray) versus outer (dark gray) wound regions. (D) Box plots showing the predicted ordering by CytoTRACE for individual cells within the four scRNA-seq clusters. Gray arrow indicates direction of predicted differentiation from scRNA-seq cluster 1 to cluster 4 (which corresponds to outer-to-inner wound region expansion). P value was derived from two-sided Student’s t test. (E) scATAC-seq evaluation of Rainbow mouse wound fibroblasts isolated in parallel with our scRNA-seq experiments (SI Appendix, Methods), integrated using the ArchR toolkit with default Louvain parameters (18) to delineate four unique multimodal fibroblast clusters. (F) Heatmap of scATAC-seq motifs highlighting key gene loci differentially open or closed in putative fibroblast subpopulations. (G) Genome tracking plots showing scATAC-seq peaks for pseudobulk replicates generated for each cluster. Associations between the peaks with fibrosis and mechanotransduction-related genes (Peak2GeneLinks) are included at the Bottom of each plot. Pale orange shading highlights differentially expressed peaks across the scATAC clusters. All highlighted peaks demonstrated statistically significant differential expression in at least one pairwise comparison (false discovery rate [FDR] <0.1 and fold change [FC] ≥2).

Fig. 3.

Clonal proliferation of injury-responsive fibroblasts is dependent on mechanotransduction signaling. (A) Representative confocal images of sectioned Rainbow mouse wound specimens treated with FAK_i (Second), FAK^fl/+ (Third), or FAK^fl/fl (Bottom) compared with vehicle control (Top). Imaris rendering in second column of images highlights individual Rainbow clones. Dermal wound area highlighted with thick white dotted line. n = 5. (Scale bars, 25 μm.) (B) Quantitation of average clone size based on Imaris rendering. (C) Wedge sections of representative whole-mount confocal images of Rainbow wound specimens embedded within surrounding wound schematics for vehicle control (Top), FAK_i-treated (Second), FAK^fl/+ (Third), and FAK^fl/fl (Bottom) samples. Corresponding vector analyses are provided to the Right of each subpanel. (D) Schematic illustrating our approach to deconvolve bulk RNA-seq data using our multimodal scRNA–ATAC construct. Transcriptionally defined cluster labels from scRNA-seq analysis were projected onto the scATAC-seq manifold using an anchor transfer–based approach in ArchR as previously described (18) (Left column) to construct four multimodal fibroblast subpopulations. Putative names were assigned to these ArchR-clusters based on integrated functional and temporospatial characteristics. Feature and peak plots, above and below, for FAK (Ptk2) are provided for illustrative purposes (Center column). Deconvolution of bulk RNA-seq specimens representing wound fibroblasts treated with FAK_i versus vehicle control (Right column) was then performed using CIBRERSORTx (19) (SI Appendix, Methods). Wound schematics (with silicone ring around the outside, and outer and inner regions indicated) are provided to represent CIBRERSORTx output identifying changes in the percentages of ArchR-cluster 1 (mechanofibrotic) cells in bulk samples over time and with/without FAK_i treatment (shown in green). Parallel schematic of corresponding changes in other ArchR-clusters are provided in yellow.

Fig. 4.

Spatial transcriptomics applied to wound healing and tracking of fibroblast subpopulations over time and space. (A) Schematic for generating spatial transcriptomics data from splinted excisional wounds using the 10× Genomics Visium protocol. Fresh Rainbow mouse wound tissue was harvested, flash frozen, embedded in optimal cutting temperature (OCT), and then sections were taken representing the complete wound radius. H&E staining and tissue section imaging were completed as described in the Visium protocol (SI Appendix, Methods). Each spot captures mRNA from 1 to 10 individual cells at that tissue location. (B) Delineation of scar layers based on underlying tissue histology at each timepoint (Top row), and UMAP plot showing that the three scar layers can easily be distinguished by their transcriptional programs, even independent of spatial information. (C) (i) Schematic of classic stages of wound healing evaluated at POD 2, 7, and 14 relative to uninjured skin. (ii) Keratinocyte activity as measured through expression of the Krt6b gene. (iii) Fibroblast activity as measured through expression of the Pdgfra gene. (iv) Immune cell activity as measured through expression of the Msr1 gene. (D) Anchor-based integration of scRNA-seq populations (defined in Fig. 2B) with Visium gene expression to project partial membership within each spot across all timepoints. These populations exhibit strong spatial preferences within the wound.

Fig. 5.

Integrated analysis permits imputation of spatial epigenomic properties. (A) Punnett square schematic summarizing the data acquired in Figs. 2 and 4; setting the stage for imputation of spatial epigenomics. (B) Schematic summarizing imputation of spatial epigenomics. Multimodal scRNA–ATAC fibroblast data were first reclustered into a higher-resolution space to generate 20 partitions, each representing between 27 and 552 cell equivalents. Gene score matrix distributions, informed by both modalities, were then extracted for each partition and subjected to SCT transformation. “Spike-in” RNA-seq data for keratinocytes, endothelial cells, granulocytes, and macrophages were obtained from pure Visium spots across all timepoints. These data were combined and subjected to a similar variance-stabilizing transformation. The resulting putative single-cell gene expression reference matrix was then used to assign initial partial set memberships for each spatial transcriptomic datapoint using an anchor transfer–based approach. Nonfibroblast contributions were subsequently regressed out, and a single-step spatial smoothing filter was applied to the resulting membership space, followed by renormalization. The resulting partial set memberships for each spatial datapoint were then treated as a topological vector space, onto which epigenomic peak, motif, and binding activity from the 20 scRNA–ATAC partitions can be projected. (C) Visium plots showing POD 0, 2, 7, and 14 (Top to Bottom) wound sections, imputed spatial epigenomics. For housekeeping genes such as Hprt (Top), gene imputed matrix (GIM) correlates with gene score matrix (GSM) epigenomic data and is fairly stable over space and time (Top). However, for Runx1, which we have shown to be very active within wound fibroblasts, GSM data show opening at the Runx1 motif at POD 2, which yields strong gene expression primarily among inner wound fibroblasts at POD 7 (Bottom). (D) Visium plots showing POD 0, 2, 7, and 14 (Top to Bottom) wound sections, motif deviations for genes of interest related to FAK-mediated mechanotransduction, and fibroblast proliferation including Runx1, Ets1, and Ehf.