Single-cell transcriptomics and a boundary wound model reveal temporal decoupling and microenvironmental heterogeneity in radiation-induced fibrosis

Figures

• 

  Figure 1. Dynamic remodeling of the whole-cell community and functional evolution of fibroblasts during different disease stages. (A) t-SNE clustering and cell-type annotation map of single-cell transcriptome data. (B) Stacked bar chart of the proportion of each cell type at different time points (Control, 1M, 5M). (C) Violin plot of the expression levels of fibrosis-related genes in different cell types. (D, E) Dynamic changes in the TGF-β signaling score (D) and ECM score (E) at different time points. Pairwise two-sided Wilcoxon rank-sum tests were adjusted within each panel using the Benjamini–Hochberg method (ns, adjusted P >= 0.05; **** adjusted P < 0.0001). (F) Expression level distribution of representative genes related to fibrosis and the cytoskeleton at different time points. (G) Changes in the Myofibroblast Function Score of fibroblasts at different time points. Pairwise two-sided Wilcoxon rank-sum tests were adjusted within the panel using the Benjamini–Hochberg method (ns, adjusted P >= 0.05; **** adjusted P < 0.0001).

• 

  Figure 2. Identification of fibroblast subpopulations and analysis of their functional characteristics. (A) UMAP dimensionality reduction map of fibroblast re-clustering and subgroup annotation. (B) Composition diagram of the proportion of each fibroblast subpopulation at different time points. (C, D) Nebulosa UMAP embedding-density maps of the ECM and TGF-β functional modules; the panels show density in reduced-dimensional expression space rather than anatomical tissue space. (E, F) Quantitative comparison of ECM score (E) and TGF-β signaling score (F) among fibroblast subpopulations. Pairwise two-sided Wilcoxon rank-sum tests were adjusted within each panel using the Benjamini–Hochberg method (ns, adjusted P >= 0.05; ** adjusted P < 0.01; **** adjusted P < 0.0001).

• 

  Figure 3. Inferred fibroblast differentiation trajectories, time-point-specific gene expression, and GO Biological Process enrichment. (A) Projection of seven computationally inferred fibroblast-state branches in UMAP space using Slingshot; these paths represent model-based state relationships rather than direct lineage tracing. (B) Heatmap of specifically highly expressed genes in fibroblasts at different time points (Control, 1M, 5M). (C) Volcano plot of differentially expressed genes in the 1M stage compared to the Control group, with the top 10 significantly upregulated and downregulated genes annotated. (D) Volcano plot of differentially expressed genes in the 5M stage compared to the Control group, with the top 10 significantly upregulated and downregulated genes annotated. (E) GO Biological Process enrichment of all significant fibroblast DEGs at 1M compared with Control. (F) GO Biological Process enrichment of all significant fibroblast DEGs at 5M compared with Control. For panels E-F, DEGs were defined using adjusted P < 0.05 and |avg_log2FC| > 0.25; upregulated and downregulated genes were analyzed together; all genes tested in the corresponding comparison were used as the…

• 

  Figure 4. Identification of cell communication networks and SPP1 signaling pathway at different disease stages. (A) Bubble plot showing macrophage-to-fibroblast ligand–receptor communication across Control, 1M, and 5M; detected 5M SPP1 interactions include Spp1–(Itga8+Itgb1) and Spp1–(Itga9+Itgb1). (B) Relative contribution of ligand–receptor pairs in the complete 5M SPP1 signaling pathway across all analyzed cell populations; this panel is not restricted to macrophage-to-fibroblast communication. (C) Network plot depicting the overall number and intensity of intercellular communication in the microenvironment at different time points (Control, 1M, 5M). (D) Comparison of dynamic changes in cell communication networks of the SPP1 signaling pathway at different time points (Control, 1M, 5M). (E) Violin plots showing normalized Spp1 expression in macrophages and Cd44 expression in fibroblasts across Control, 1M, and 5M. Small points represent individual cells. Pairwise annotations are based on two-sided Wilcoxon rank-sum tests with Benjamini–Hochberg correction within each expression panel. ** adjusted P < 0.01; **** adjusted P < 0.0001.

• 

  Figure 5. Impact of localized irradiation on wound healing in mice. (A) Experimental workflow: Mice received localized irradiation (12 Gy, 1 cm diameter), followed by 6-mm full-thickness wounds at 1, 5, 9, or 13 days after irradiation (R1, R5, R9, R13). In the boundary-spanning H configuration, the wound center was placed on the field circumference so that the wound approximately crossed irradiated and non-irradiated tissue; this was not an exact half-and-half division. In the W configuration, the wound was fully contained within the irradiated field. The within-animal validation model (V group) combined boundary-spanning H and fully irradiated W wounds in the same mouse to validate the contrast between the two configurations. Wounds were collected 14 days after surgery for histology. (B) Representative wound images recorded every 2 days (days 0–14). Wounds were covered with 3M Tegaderm for the first 2 days. (C, D) Wound closure kinetics over time: wound area ratios (current/initial) in NC versus boundary-spanning H wounds (C) and NC versus fully irradiated W wounds (D). Data are mean ± SEM; these longitudinal panels are descriptive. (E) Wound closure at day 14 (1 − wound area/day 0),…

• 

  Figure 6. Histological analysis of wound healing following localized irradiation in mice: collagen deposition and COL1A1 expression. (A–C) HE panoramic images of NC, R9W, and R9H wounds; red boxes indicate dermal regions analyzed. (D) Representative images at different magnifications: 10× HE and Masson (epidermis) and 40× HE (dermis). (E) Masson trichrome staining and COL1A1 immunohistochemistry in NC, boundary-spanning H, and fully irradiated W wounds. (F, G) Quantification of collagen deposition and COL1A1 expression. Multiple regions of interest were quantified per tissue section; n = 6 animals per irradiated configuration at each irradiation-to-wounding interval. R1, R5, R9, and R13 indicate wounds created 1, 5, 9, and 13 days after irradiation, respectively; all tissues were collected 14 days after wound creation. Statistical methods are described in the Statistical Analysis section.

• 

  Figure 7. Spatial and temporal distribution of macrophages and immunohistochemical analysis of wound healing. (A, B) CD68 immunostaining and quantification showing the regional and temporal distribution of macrophages. Regions: NC = non-irradiated control, W_C = central area of a fully irradiated wound, W_P = peripheral area of a fully irradiated wound, H = boundary-spanning wound. (C) Representative immunohistochemistry images for CD31, Ki67, TGF-β, SMAD2/3, and α-SMA. (D–H) Quantitative analysis of CD31, Ki67, TGF-β, SMAD2/3, and α-SMA in NC, boundary-spanning H, and fully irradiated W wounds over time. Multiple regions of interest were quantified per tissue section; n = 6 animals per irradiated configuration at each irradiation-to-wounding interval. R1, R5, R9, and R13 indicate wounds created 1, 5, 9, and 13 days after irradiation, respectively; all tissues were collected 14 days after wound creation. * adjusted P < 0.05; ** adjusted P < 0.01; *** adjusted P < 0.001.