Single-Cell Analysis Reveals Cxcl14+ Fibroblast Accumulation in Regenerating Diabetic Wounds Treated by Hydrogel-Delivering Carbon Monoxide

Nonhealing skin wounds are a problematic complication associated with diabetes. Therapeutic gases delivered by biomaterials have demonstrated powerful wound healing capabilities. However, the cellular responses and heterogeneity in the skin regeneration process after gas therapy remain elusive. Here, we display the benefit of the carbon monoxide (CO)-releasing hyaluronan hydrogel (CO@HAG) in promoting diabetic wound healing and investigate the cellular responses through single-cell transcriptomic analysis. The presented CO@HAG demonstrates wound microenvironment responsive gas releasing properties and accelerates the diabetic wound healing process in vivo. It is found that a new cluster of Cxcl14+ fibroblasts with progenitor property is accumulated in the CO@HAG-treated wound. This cluster of Cxcl14+ fibroblasts is yet unreported in the skin regeneration process. CO@HAG-treated wound macrophages feature a decrease in pro-inflammatory property, while their anti-inflammatory property increases. Moreover, the TGF-β signal between the pro-inflammatory (M1) macrophage and the Cxcl14+ fibroblast in the CO@HAG-treated wound is attenuated based on cell–cell interaction analysis. Our study provides a useful hydrogel-mediated gas therapy method for diabetic wounds and new insights into cellular events in the skin regeneration process after gas-releasing biomaterials therapy.

Streptozotocin was purchased from Meilun Biotech.(Dalian, China).Silver nitrate was purchased from Tianjin Kemiou Chemical Reagent Co. Ltd. (Tianjin, China).PBS buffer and paraformaldehyde fixative (4%) were purchased from Servicebio Service BioTechnology Co. Ltd. (Wuhan, China).Optimum Cutting Temperature (OCT) compound was purchased from Applygen Technologies Inc. (Beijing, China).Above commercial reagents and solvents were used as received without further purification.RAW264.7 cell, NIH3T3 cell and cell culture medium were purchased from Procell Life Science and Technology Co, Ltd (Wuhan, China).3,3'-Dithiobis (propionic hydrazide) was synthesized according to the previous literature. 1 Preparation of HAase cleaved HA-SH solution: 3 mL of 5 mg/mL HA-SH solution was mixed with 10 μL of 100 mg/mL HAase solution, and the mixture was incubated at 50 °C for 45 min to cleave HA-SH.Subsequently, the HAase-cleaved HA-SH solution was incubated at 70 °C for 1 h to inactivate HAase.

Synthesis of HA-SH derivative
S3 400 mg of HA (corresponding to 1 mmol carboxyl group) were dissolved with 60 mL H2O followed by addition of 16.7 mg (0.07 mmol/L) of 3'3-Dithiobis(propionic hydrazide), 135 mg of HoBt (1 mmol) in 3 mL of 1:1(V/V) acetonitrile: deionized water, 57.5 mg of EDC (0.3 mmol).The above mixture was reacted at room temperature overnight under pH 4.5-5 environment, thereafter 54.5 mg of DTT (0.25 mmol) was added into the reaction mixture to cleave disulfide bond and obtain HA-SH.The structure of resulting derivative after dialysis and lyophilization were measured with 1 H-NMR spectra (400 MHz, D2O).The new appearance of peaks at 2.80 and 2.65 ppm in the 1 H-NMR curves of HA-SH derivative were attributed by protons of methylene group (-CH2CH2SH) side chain on HA backbones (Figure 1b).The thiol group modification density in HA-SH polymer was around 10%, calculated by the integration of the above methylene protons compared with the acetamide group of the acetylglucosamine unit in the HA backbone.

Fabrications and characterizations of CO@HAG and HAG
CO@HAG was formed by simply mixing HA-SH sol (2% (w/v)), Mn2[CO]10 dispersion (1.5% (w/v)), and AgNO3 solution (15 mM).Particularly, to prepare 1 mL of CO@HAG hydrogel, firstly 20mg HA-SH was dissolved into 700 μL H2O and 15 mg Mn2[CO]10 was dispersed in 200 μL H2O.Thereafter HA-SH solution and Mn2[CO]10 dispersion were quickly mixed and 100 μL of 150 mM of AgNO3 solution was immediately added into the above mixture under vortexing.The hydrogel without Mn2[CO]10 was represented as HAG.Clinical used Ag + ion (AgNO3) solution of concentration is 29 mM and above 60 mM of Ag + ions are reported to be toxic to normal tissue. 2,3 he 15 mM Ag + ions inside the presented CO@HAG is four-ford less than the toxic concentration, and the same concentration was also safely used in our previous paper. 4Rheological properties of CO@HAG and HAG were investigated using MCR-92 advanced rheometer (Anton Paar, Austria) with a 15-mm-diameter parallel plate.
Storage modulus (G') and loss modulus (G") of hydrogels containing different Mn2[CO]10 concentrations were measured by frequency sweep experiments (from 1 to 100 rad/s) at a fixed strain (0.5%).Moreover, CO@HAG and HAG of G' and G" were monitored under the strain's gradual increase from 1% to 100% at a fixed 1 Hz frequency.SEM (TESCAN, Czech Republic) was used to observe the microstructure of the lyophilized CO@HAG and HAG.Mass-remaining ratio curves for CO@HAG and HAG were investigated in PBS medium with pH 7.4, and the mass-remaining ratio (%) (n = 3) was calculated by dividing hydrogel mass after each time point of incubating to the initial mass before incubation.

Hydrogels anti-bacterial and H2O2-scavenging experiments
E.coli (ATCC 25922, Gram-negative strains) and S.aureus (ATCC 25923, Grampositive strains) were used to evaluate the anti-bacterial properties of CO@HAG and HAG. 1 mL of bacterial solution with 10 7 CFU/mL was co-cultured with 200 µL of CO@HAG and HAG on a 24-well plate for 12 hours under 37℃.200 µL of PBS instead of hydrogel as control group was cultured with the above bacterial solutions.CO@HAG, HAG and PBS-treated bacterial solutions were diluted with the same dilution ratio and further cultured on Luria-Bertani (LB) agar medium for 12 hours.
Thereafter the remaining colonies were observed using optical photographs.To quantify hydrogels-treated remaining bacterial viabilities, the OD values at 450 nm were recorded by SpectraMax M2 instrument (Molecular Devices, US) according to CCK8 assay.We investigated hydrogels of non-radical ROS (H2O2) scavenging abilities using the Fe 2+ -Fe 3+ transition of Fenton reaction method.Particularly, CO@HAG and HAG were incubated with 1 mL 100 μM H2O2 solution and the residual H2O2 amounts after 3 and 6 hours of incubating were measured using a H2O2 quantitative assay kit (Sangon Biotech, Shanghai).

CO gas release measurement
The released CO amount was measured using the hemoglobin (Hb)carboxyhemoglobin (HbCO) conversion method.Particularly, 4.2 μM of bovine hemoglobin was dissolved in PBS buffer and reduced with 1.6 mg of sodium dithionite (SDT) under a nitrogen protection environment.200 μL of CO@HAG was incubated with 4 mL reduced Hb solution containing 5 μM of H2O2.Immediately, the above reaction mixture was sealed in a UV quartz cuvette and the UV absorbances were recorded from 350 nm to 600 nm using a UV-Vis spectrophotometer (UV-6100A, Metash Instruments Co.Ltd., China) at the certain time-points.200 μL of HAG under 5 μM H2O2 medium was used as negative control.Moreover, the UV absorbance at 410 and 430 nm for the CO@HAG under different H2O2 concentrations of medium (0.05, 2.5, 3.75, 5 μM) were recorded to calculate the CO releasing amount using the Figure S1f equation. 5I410nm and I430nm represent the OD value of the collected spectrum at 410 nm and 430 nm, respectively.CCO and CHb represent releasing CO amount and the initial Hb amount (4.2 μM), respectively.The membranes furtherly were incubated with HRP-conjugated Goat anti-Rabbit IgG (H+L) (Abiowell, AWS0002a, 1:10000) for 1 hour.The targeted protein bands were visualized on Automatic Gel Imaging Analysis System with the presence of enhanced chemiluminescence (ECL) reagents (Thermo Scientific, 34580).As for cytometry measurements, the treated cells were also analyzed using flow cytometry (CytoFLEX, Beckman, USA) after incubating with the PE anti-Nos2 (iNOS) Antibody (1:500, Biolegend, 696806) overnight.

In vivo diabetic wound healing study
SD male rats (6 -8 weeks of age; Hunan Sta Laboratory Animal Co., Ltd., China) were intraperitoneally injected streptozotocin (65 mg/kg).The rats with blood glucose over 16.7 mM were used to prepare skin full-thickness defects with 15 mm of diameter under anesthesia.The prepared 300 μL of CO@HAG or HAG was painted onto each wound area and the wounds were protected with sterilized gauze.The same volume of PBS was used as untreated group.Digital images of the wound area were captured with photo camera at post-therapy of 0, 3, 7, 10 and 16 days, and the wound closure rate was analyzed using Image J software.Wound tissues were obtained and embedded in the OCT compound after fixation with 4% paraformaldehyde.For the histological analysis,

Skin immunofluorescence staining
All sections were washed three times with PBS and permeabilized in 0.1% TritonX-100 for 15 min at room temperature, followed by being blocked with 5% FBS for 1 hour at room temperature.The blocked sections were incubated with primary antibody at 4°C overnight.After washing three times with PBS, all sections were incubated with corresponding secondary antibodies for 1 hour and DAPI dye for 15 min at room temperature.The images of skin tissue sections were captured with an Ultra-high resolution confocal microscope (LSM980, ZEISS, Germany).Herein used primary antibodies included anti-iNOS antibody (ab178945, 1:400) and anti-CXCL14 antibody (ab264467, 1:100).Secondary antibody (from Invitrogen) used in the study was Alexa Fluor 546-conjugated Goat anti-Rabbit Secondary Antibody (A-11035, 1:1000).

Skin real-time qPCR analysis
Harvested skin tissues were washed with PBS and total RNA was extracted with Trizol (Solarbio, China) according to the protocols.The concentration of extracted total RNA was determined using Nano Drop spectrophotometer (Nanodrop2000, Thermo Fisher Scientific).Complementary DNA (cDNA) was synthesized by reverse transcription of the total RNA (1 μg) using reverse transcriptase (Bimake, USA).The relative mRNA expression levels (Il1b, Il6, Tnf, Il4, and Il10) were quantified using Hieff® qPCR SYBR Green Master Mix (Yeasen, China) by real-time fluorescence quantitative PCR instrument (CFX96, Bio-rad, USA).Data were analyzed using the 2 -ΔΔCt method.The cycle threshold (Ct) values of target gene were normalized using the Ct values of Gapdh gene.Primer sequences were listed in Table S3.

Single-cell capture and scRNA-Seq library preparation
Fresh wound skin samples were obtained from 7 days of untreated, HAG and CO@HAG groups, cleaned with PBS, and digested using the Skin Dissociation Kit (Miltenyi biotec, 130-101-540) for ~1 hours at 37℃.Then the filter was used to purify the above cell suspension, and ACK buffer (Lonza, 10-548E) was used to remove red blood cell.The resulting cells were necessary to reach >90% of cell viabilities and resuspended in PBS containing 0.04% Bovine Serum Albumin (BSA).The cells were immediately captured by Single Cell 3' Library and Gel Bead Kit V3.1(10x Genomics, 1000121) and Chromium Single Cell G Chip Kit (10x Genomics, 1000120).The cell suspension with 300-600 living cells/μL was loaded onto the Chromium single cell controller (10x Genomics) to generate single-cell gel beads in the emulsion.The captured cells were lysed inside gel bead-in-emulsions (GEMs) and the released RNA was barcoded through reverse transcription using S1000TM Touch Thermal Cycler (Bio Rad) at 53℃ for 45 min, followed by 85℃ for 5 min, and hold at 4℃.The generated cDNA was amplified, and quality estimated using Agilent 4200 (performed by CapitalBio Technology, Beijing).The sequences of libraries were finally obtained using the Illumina Novaseq6000 sequencer with a sequencing depth of 25,000~40,000 reads per cell and pair-end 150 bp (PE150).

Reads mapping and gene expression quantification
Single-cell cDNA library samples were sequenced as described above.Demultiplexing, alignment, and quantification of single-cell fastq files were performed using the Cell Ranger Pipeline v7.0.0 (10x Genomics) 6 with default settings against the Rattus norvegicus reference genome (v.6.0, release 104, downloaded from the Ensemblwebsite). 7

Data quality control and normalization
Low-quality cells were removed based on the number of expressed genes (nGene) and the expression level of mitochondrial genes (percent.mito).Specifically, cells with percent.mitoless than 0.15 and 750 < nGene < 6000 were retained.Doublet discrimination was implemented using the doubletFinder_v3 function in the DoubletFinder package (v.2.0.3) for each sample individually. 8Following the selection of high-quality cells, we further filtered based on genes.Genes that excluded mitochondrial genes were considered truly expressed if they contained one or more counts in at least five cells (assessed for each sample separately).Log-normalized counts were calculated using the deconvolution strategy implemented by the computeSumFactors function in the scran package (v.1.14.6). 9 Rescaled normalization was performed using the multiBatchNorm function in the batchelor package (v.1.2.4)   to ensure that the size factors were comparable across samples. 10The log-normalized expression after rescaling was used in marker gene detection and differential gene expression analysis.

Data integration, dimensionality reduction, and clustering
We integrated the filtered count matrices from 3 samples using the sctransform approach implemented in the Seurat package (v.4.2.0) on the 3000 anchor features. 11ter integration, principal component analysis (PCA) was performed on the integrated data followed by embedding into low-dimensional space with Uniform Manifold Approximation and Projection (UMAP) based on the top 40 dimensions.Clusters were generated by graph-based method using the FindClusters function from the Seurat package and assigned to cell types by consulting the expression of known marker genes and automatic annotation from the SingleR package (v.1.4.1). 12Proliferating cells were identified by the high-level expression of Top2a, and Stmn1 and excluded from further analysis.To identify detailed types of fibroblast cells, cells from the 2 big fibroblast clusters were further extracted and integrated based on 2000 anchor features and top 30 dimensions.Sub-fibroblast clusters were further identified using the FindClusters function based on fibroblast-only PCA dimensions and by requiring at least one highly expressed marker (power > 0.4 using 'roc' test) (Figure 5; Figure S12, Supporting Information).Raw counts for cells belonging to the same fibroblast subclusters were aggregated into a pseudo-bulk sample.Normalization of expression for pseudo-bulk samples was calculated using the estimateSizeFactors function, followed by the estimateDispersions function from the R package DESeq(v.1.38.0). 13The logtransformed and normalized gene expression values of expressed genes from pseudobulk data were used to calculate the Spearman correlation among different subclusters, followed by hierarchical clustering by the hclust function in R.

Marker gene detection and differential gene expression analysis
Marker genes for each cell type were identified using the FindAllMarkers function with the "roc" test from the Seurat package.The top 50 marker genes with an average power of at least 0.4 were selected (Table S1).Differential gene expression analysis among samples in the macrophage clusters was performed using the "MAST" test implemented in the FindMarkers function from the Seurat package. 14Genes with an FDR less than 0.05, log2(fold change) greater than 0.25, and expressed in more than 15% of the cells were considered differentially expressed.The scaled gene expression in the untreated, HAG, and CO@HAG samples were checked by the radial ggplot function in the volcano3D package. 15The gene ontology (GO) enrichment analysis for differentially expressed genes was conducted using the TopGO package (v.2.42.0). 16Adrian Alexi's improved weighted scoring algorithm and Fisher's test were used to define the significance of GO term enrichment.Functional enrichment analysis and gene set enrichment analysis (GSEA) were performed using the "enricher" and "GSEA" function from the clusterProfiler package(v.3.18.1), 17 respectively.9][20] Significantly enriched GO and functional terms were identified as those with a p-value and FDR less than 0.05, respectively.

RNA velocity, trajectory interference, and pseudotime analysis
RNA velocity was performed using velocyto package(v.0.17.17) 21by running "velocyto run10x" and further processed by scvelo package (version 0.2.4) 22 based on spliced and unspliced transcript reads, as previously reported. 22The dynamical modeling (scv.tl.velocity) method was used, and only cells belonging to fibroblasts were included in the analysis.The "scv.pl.velocity_embedding_stream" function was used to project RNA velocities onto UMAP plots.All default parameters were used unless noted otherwise.The "scv.tl.paga" function was used to evaluate the relationship between different cell clusters using the partition-based graph abstraction (PAGA) analysis. 23Trajectory analysis of CO@HAG dermal (Dpt + and Plac8 + ) fibroblast cells was performed using the R monocle package (v.2.18.0) 24,25 with the DDRTree method and default parameters.Raw UMI count was used as the input for Monocle2.After default normalization, the top 500 variable genes were selected by the FindVariableFeatures function, followed by reconstructing the single-cell trajectory, estimating functional states of cells, and calculating pseudotime using the reduceDimension and orderCells function.Variation from the effects of "number of expressed genes" was subtracted by setting the residualModelFormulaStr argument to "~nGene".

Analysis of cell-cell interactions
To identify and visualize cell-cell interactions, we employed the R package CellChat (v.1.5.0) 26 to perform cell-cell communication analysis.Briefly, we followed the official workflow((http://www.cellchat.org/)),loaded the normalized counts into CellChat, and applied standard preprocessing steps, including those using the functions "identifyOverExpressedGenes," "identifyOverExpressedInteractions," and "projectData" with a standard parameter set.We selectively used a total of 2,017 precompiled mouse ligand-receptor interactions (named "CellChatDB.mouse")as a priori network information.In each individual sample, we calculated the potential ligand-receptor interactions among cells based on the functions "computeCommunProb," "computeCommunProbPathway," and "aggregateNet."The distance of signaling networks between untreated and CO@HAG samples and network centrality scores were calculated using the function "rankSimilarity" and "netAnalysis_computeCentrality," respectively.The function "netAnalysis_signalingRole_network" was applied to the network to determine the senders and receivers.Significant interactions associated with the TGF-β signaling pathway were extracted using the function "extractEnrichedLR".We used all standard parameters unless otherwise noted.

Statistical analysis
Results are presented as the mean ± standard deviation (SD) or mean ± standard error of the mean (s.e.m.).Statistical significance between the different groups were evaluated using one-way ANOVA with posthoc multiple comparisons test in IBM SPSS Statistics 26.0 software.All data were collected from at least three repeated samples(n≥3) unless otherwise noted.
the 10 μm of frozen slices were stained with Haematoxylin & Eosin (H&E) and Masson's Staining Kit.Images of stained slices were captured using a Digital Slice Scanning System (Pannoramic MIDI, 3DHISTECH, Hungary).All experimental procedures involving animals were approved by the Institutional Animal Care and Use Committee of College of Biology, Hunan University (Approval NO.HNUBIO202101006).

Figure S5 .
Figure S5.Gas chromatography analysis of CO release from CO@HAG and HAG in the presence of H2O2 (5 μM).

Figure S6 .
Figure S6.Flow cytometry analysis results.(a, b) Flow cytometry analysis of iNOS expression of LPS-induced macrophage treated by HA-SH and Mn2[CO]10.(c) Based on flow cytometry analysis, the value of average inhibition rate of HA-SH and Mn2[CO]10 against iNOS-positive macrophage (n = 3; mean ± s.e.m.).

Figure S8 .
Figure S8.CO@HAG regulated regenerating wound of cellular heterogeneity.(a) Uniform Manifold Approximation and Projection (UMAP) plot of all cells from the untreated, HAG, and CO@HAG samples, respectively.(b) Expression of fibroblast and macrophage marker genes (RGD1559482, Cd14, Dcn, and Pdgfra) overlaid onto the UMAP plot, as shown in Fig.3h.(c) Barplots showing the proportion of different cell types in the untreated, HAG, and CO@HAG samples, respectively.

Figure S9 .
Figure S9.CO@HAG and HAG differentially regulated wound macrophage of inflammatory responses.(a) Venn diagrams showing overlapping up-regulated and down-regulated DEGs in the CO@HAG sample for four types of macrophages compared with untreated sample.(b) Volcano plot showing the fold change in expression and adjusted P-values (FDR) of DEGs in macrophage cells when comparing CO@HAG and untreated samples.The DEGs that are also differently expressed in

Figure S13 .
Figure S13.DEGs between CO@HAG and untreated in Fib_SC6 cells (a), along with the corresponding top 10 enriched biological pathway Gene Ontology (GO) terms (b).

Figure S14 .
Figure S14.Cell-cell interactions in HAG group.(a) Hierarchical plot showing