- 1Department of Biological Sciences and The RNA Institute, University at Albany, State University of New York, Albany, NY, United States
- 2Molecular, Cellular, Developmental, and Neural Biology Graduate Program, Department of Biological Sciences, University at Albany, State University of New York, Albany, NY, United States
Extracellular matrix remodeling is a natural response to injury but, excessive extracellular matrix accumulation, or fibrosis, is a causative factor in hundreds of diseases that limit organ function, regenerative responses, and can interfere with regenerative therapies. Fibrosis is closely related to inflammation, both of which occur in the salivary glands of patients treated with radiation for head and neck cancers and in patients suffering from autoimmune conditions, such as Sjögren's Disease. Despite the known involvement of fibrosis in disease and the inhibitory effects of fibrosis on tissue regeneration, the mechanisms through which extracellular matrix is elaborated in the salivary gland are poorly understood. Stromal fibroblasts are the primary matrix-producing cells and are known to drive both fibrosis and inflammation. To define the temporal responses of fibroblasts to injury, we induced a temporary obstructive injury though ligation of the primary submandibular and sublingual salivary gland ducts and then performed single-cell RNA sequencing and pathway analysis at timepoints immediately following the injury. Using bioinformatic approaches, we identified three unique fibroblast groups that dynamically respond to the injury. We characterized the changes in matrisomal and inflammatory gene expression over a 7-day time course and identified one group of fibroblasts to be the primary injury-responsive fibrogenic cell type. Understanding how fibroblasts respond at the early and later injury timepoints, along with defining signaling pathways regulated by fibroblasts, could lead to a better understanding of the contribution of fibroblast to acute injury responses to facilitate the development of therapeutics that minimize fibrosis and promote regenerative gland responses in chronic disease states.
1 Introduction
The term fibroblast is often used to describe a heterogenous cell population that performs a multitude of roles including extracellular matrix (ECM) elaboration, tissue maintenance and repair, and crosstalk with various cells present in organs (1). Indeed, fibroblasts are the primary cell type that elaborates and remodels the structural and regulatory components of the ECM, known collectively as the matrisome (2). Fibroblasts are known to change their gene expression in response to the environment, and this plasticity is critical to allow fibroblasts to contribute to changing tissue physiology and alter their communication with other cell types. In response to injury or disease, a fibroblast population, referred to as myofibroblasts or fibrogenic fibroblasts, produce the increased and sometimes excessive ECM proteins that make-up fibrosis (3, 4). Fibroblasts are heterogeneous, and recent scRNA-Seq studies highlight the diversity of fibroblast subpopulations present in different organs and in disease states (5–7, 58).
In salivary glands, fibroblasts are essential for organ development and maintaining organ homeostasis in adults, as they exhibit temporal and spatially dynamic deposition and modification of the matrisome. A hallmark of salivary gland damage and disease-associated hyposalivation is the development of a fibrotic response. Importantly, a well-orchestrated fibrotic response is an essential aspect of normal tissue repair and wound healing, while an aberrant and excessive response contributes to progressive organ dysfunction and degeneration (8). This response can be caused by radiation damage, chronic inflammation, or other causes (9–13) and likely contributes to reduced gland function. Recent work has highlighted the heterogeneity in the severity of fibrosis in the submandibular glands of patients with sclerosing sialadenitis (14). Although there have been many recent advances in the development of therapeutic approaches to treat salivary gland hypofunction (15), how to limit salivary gland fibrosis has not been determined. How fibroblast plasticity and function contribute to healing verses pathological salivary gland fibrosis remains unclear, but this knowledge could inform the development of improved therapeutics.
In this study, we examined the temporal changes in fibroblast subpopulation transcriptomes using the mouse submandibular salivary gland ductal ligation injury model, which exhibits a progressive Transforming Growth Factor Beta (TGFβ)-dependent fibrotic response concomitant with gland degeneration (16–19). This injury model was first developed to induce secretory dysfunction in rats that is reversible upon removal of the ligature (20). In prior work, we examined the fibrotic response that occurs 2 weeks after ligation injury in 12 week old female mice (18). Here, we performed single-cell RNA sequencing (scRNA-seq) to identify changes in the fibroblast cell populations immediately following ligation from 1 to 7 days to examine the temporal changes in the construction of the ECM following this acute injury.
2 Materials and methods
2.1 Animal husbandry
All animal husbandry, surgical procedures, and tissue collection were performed in accordance with protocols approved by the University at Albany, SUNY IACUC committee. Mice were housed in 12-hour light/dark cycle with access to water and dry food. C57BL/6J (JAX #000664) mice were purchased from The Jackson Laboratory. The founders for the B6.129S-Pdgfratm1.1(cre/ERT2)Blh/J (JAX #032770) and B6.Cg-Gt(ROSA)26Sor-tm9(CAG-tdTomato)Hze/J (JAX #007909) mouse colonies were obtained from Jackson Laboratories. To create double transgenic Pdgfrα reporter mice, The PdgfrαCreERT2 +/wt (PdgfrαCreERT2) males were crossed with ROSA26TdTomato +/+ females to generate PdgfrαCreERT2; R26tdT +/wt (PdgfrαCreERT2; R26tdT) mice. Mice were assigned a unique identifier between postnatal day 7 and 10 and were genotyped with PCR to detect Cre and tdT in PdgfrαCreERT2 and PdgfrαCreERT2R26tdT mice.
2.2 Tamoxifen induction
For lineage tracing experiments, PdgfrαCreERT2; R26tdT young female mice were induced at 10–12 weeks old with 3 × 100 µg/g bodyweight tamoxifen (Sigma T5648) injections, as previously described (21). The tamoxifen was first dissolved in 100% ethanol at 55°C then added into corn oil (Sigma C8267). The final concentration was 10% ethanol in corn oil. Injections were performed via i.p. starting 7 days before surgery with three consecutive injections, one every other day.
2.3 Salivary gland ductal ligation
We performed pre-operative and operative procedures, as previously described (18). Mice were continuously monitored under anesthesia and post-operatively for pain, distress, and changes in weight for a minimum of 48 h following surgery, in applicable cases. Wharton's and Bartholin's ducts were ligated for 1, 3, or 7 days. As negative controls for the surgical experiments, induced PdgfrαCreERT2; R26tdT mice underwent no surgical manipulations and were euthanized at 12-weeks old. All mice were euthanized at the desired time points under CO2 with secondary cervical dislocation.
2.4 Single-cell isolation for scRNAseq
For enrichment of stromal cell populations, three mice were used for each timepoint with the mouse's right submandibular and sublingual gland both being harvested together. Excess fat and interstitial tissue were removed. The glands were then transferred to a dish containing 1X phosphate buffered saline (1XPBS), liberase TL Research Grade low Thermolysin (Roche 05401020001) at 0.25 mg/ml, DNase I (Stemcell Technologies 07900) at 15 mg/ml, and dispase (Gibco 17105-041) at 0.53 U/ml and micro dissected for 10 min. The sample was then incubated in a 37°C incubator for 15 min then triturated 100 times. After trituration, the sample was incubated for an additional 5 min in a 37°C incubator and triturated once more. A 15 ml conical tube was placed on ice and the entire sample was transferred to the conical tube and incubated for 10 min. The supernatant was isolated and transferred to a fresh conical tube containing an equal volume of DMEM/F12 (Gibco 11039-021) with 10% fetal bovine serum (Gibco 10082-147) and was centrifuged for 5 min at 150 × g 4°C. The supernatant was removed and discarded, and the cell pellet was resuspended in 2 ml of Dynabead isolation buffer composed of sterile-filtered Ca2+- and Mg2+-free 1X PBS, 0.1% w/v bovine serum albumin, and 2 mM ethylenediaminetetraacetic acid. 5 μg of EpCAM (Invitrogen 14-5791-81) and Ter119 (Invitrogen 14-5921-82) monoclonal antibodies were added to the cell suspension and incubated at 4°C on a rocker set on low speed for 10 min. The cell suspension was washed using isolation buffer and centrifuged at 150 × g for 5 min at 4°C. The supernatant was discarded, and the cell pellet was resuspended in 1 ml of isolation buffer.
For depletion of epithelial and red blood cells, 25 μl of sheep anti-rat Dynabeads, which bind to EpCAM- and Ter119 antibody-labeled cells, (Invitrogen 11035) were added. The cell suspension was transferred to a 35 mm dish before incubating at 4°C for 20 min on a rocker set at low speed. The sample was placed on a microcentrifuge tube magnet for 2 min, and the supernatant was collected. One additional epithelial and red blood cell depletion was performed. The sample was then depleted of dead cells using a dead cell removal kit (Miltenyi Biotec 130-090-101), following manufacturer's instructions. Cells were resuspended to a concentration of 1,000 cells/μl and counted using a TC20 automated cell counter (BioRad). Following the manufacturer's protocol for the Chromium Next GEM Single Cell 3' Reagent kits v3.1 (Dual Index) user guide, scRNA-seq libraries were generated.
2.5 Single-cell RNA-sequencing analysis
All harvested samples were sequenced on an Illumina NextSeq 2000 at the Center for Functional Genomics - University at Albany, Albany, NY, RRID:SCR018262. Initial processing steps included aligning to the GRCm39 genome with the addition of the tdTomato transcript (22) and generating counts files, which was completed using CellRanger version 8.0.0. Background noise was removed using CellBender (23). Data files were imported using Seurat v5.0.1 in R v4.3.1 (24, 25). Dead or apoptotic cells were removed if >25% of unique molecular identifiers (UMIs) mapped to mitochondrial genes. Any cells with less than 200 or over 9,000 genes were excluded. Potential doublets were detected by DoubletFinder (26) and scDblFinder (27). Droplets that both packages agreed on were removed as doublets. Integration was performed on 12 datasets to create one dataset. The dataset was normalized with SCTransform (28) after integration. Clusters were calculated following the default pipeline (29). All analysis of the integrated dataset was performed using the Seurat package in combination with other R packages. Annotated codes used for analysis are available on GitHub (https://github.com/MLarsenLab).
3 Results
3.1 Single cell RNA sequencing of ligated submandibular glands
To induce a progressive fibrotic response, we performed the ductal ligation reversable injury model to female 10–12-week-old mouse submandibular and sublingual salivary glands by applying a metal clip to Wharton's and Bartholin's ducts, respectively, proximal to the glands (Figure 1A). We used PdgfraCreERT2; R26tdT mice and collected glands from the homeostatic pre-injury state and at timepoints shortly after injury (e.g., 1-, 3- and 7-days post-injury) to reveal early dynamics in the Pdgfra-expressing fibroblasts that presumably drive the fibrotic response (Figure 1A). After gland dissociation and enrichment, the stromal cell populations in the tissue were subjected to drop-seq single-cell RNA sequencing. After using Seurat to perform unsupervised clustering, the expected cell clusters were detected in homeostatic glands and ligated glands across all surgical timepoints (Figure 1B,C and Supplementary Figure S1A). Significantly, we detected a progressive increase in the relative contribution of the fibroblast population during the injury time course, that peaked at 7 days. Inflammatory cells, including T cells and macrophages, peaked at 3 days, and NK cells peaked at 7 days (Figure 1D,E). The fibroblasts were then isolated with the subset command, selecting cells that were labeled as fibroblasts by their Pdgfra expression or cells that expressed tdTomato, creating a new Seurat object. After performing unsupervised clustering, the Seurat object contained 13 distinct subclusters (Figure 2A, and Supplementary Figure S2). Differential gene expression analysis revealed transcriptomic similarities and differences among specific fibroblast subclusters that were then binned into 3 distinct groups based on these similarities (Figure 2B,C).

Figure 1. Generation of scRNA-seq dataset in salivary glands 1-, 3-, and 7-days post-ligation. (A) Schematic describing the ductal ligation surgery procedure where a clip is placed on ducts exiting the SMG and SLG of young female mice. Tissue was harvested at 4 timepoints (n = 3) and subjected to stromal enrichment through gravity sedimentation and MACS after enzymatic dissociation. (B) A UMAP displays the cell populations present after processing the datasets and integration. (C) A violin plot shows the expression of marker genes used to identify cell populations in the Seurat object. (D) A split UMAP shows the respective contributions of each timepoint to the Seurat object. (E) Percentages of cells which contribute to cell populations at each timepoint.

Figure 2. SMG fibroblasts are a heterogeneous population. (A) UMAP displays the fibroblast subset obtained by selecting for cells labeled as fibroblasts or expressing tdTomato. (B) UMAP shows the 3 groups used to bin clusters that express similar transcriptomes. (C) Heatmap shows selected representative genes which highlight the transcriptomic similarities and differences that were used to justify grouping fibroblast subclusters into groups.
3.2 Characterization of the fibroblast subgroups
Heatmaps and GO term analysis of the top 50 differentially expressed genes in each fibroblast group highlight their unique transcriptome profiles that define them as groups and are suggestive of functional differences between the groups (Figure 3, Supplementary Tables S1–S3). Group 1-enriched transcripts exhibit cellular communication and morphogenesis signatures (Figure 3A,D). Cells in this cluster express Piezo2 and Transforming Growth Factor Beta 1 (Tgfb1). Platelet Derived Growth Factor Receptor Beta (Pdgfrb) is a profibrotic mediator that is uniquely expressed in this population along with Secreted Phosphoprotein 1 (Spp1). Interleukin 34 (Il34) is also expressed in this cluster. Together, these data suggest that the group 1 fibroblasts are responsive to mechanical changes in the ECM composition and are both pro-inflammatory and pro-fibrotic.

Figure 3. Differential gene analysis in fibroblast groups. (A–C) Heatmaps show the top 50 differentially expressed transcripts by each group. Group orders are from top to bottom; Group 1, Group 2, and Group 3. (D–F) GO term analysis performed using ShinyGO with the top 500 differentially expressed genes. Order from top bottom is Group 1, Group 2, and Group 3. Refer to Supplementary Tables S1–S3 for complete differentially expressed gene lists.
Group 2-enriched transcripts exhibit an over 20-fold enrichment in cytoplasmic translation and an extracellular matrix signature with enrichments ranging from 5 to 10 fold in ECM-related categories (Figure 3B,D). In this cluster, Collagen Triple Helix Repeat Containing 1 (Cthrc1) is expressed. Secreted Frizzled Related Protein 2 (Sfrp2) and Dermatopontin (Dpt) are also expressed. We see increased expression of Loxl2 expression which is a member of the lysyl oxidase family proteins that crosslink collagens and elastin to stabilize ECM structure during fibrosis (30, 31). Thus, group 2 fibroblasts are likely to be the primary mediators of ECM deposition and structural modification during the fibrotic response to injury.
Group 3-enriched transcripts exhibit metabolic and synthetic signatures, suggesting that these cells may be critical for directing and maintaining tissue structure and function (Figure 3C,F). Noteworthy genes expressed include Leukemia Inhibitory Factor Receptor (Lifr), which has recently been described as a master amplifier of fibrosis in idiopathic pulmonary fibrosis (32). Fibroblast Growth Factor 10 (Fgf10), which is a critical mediator of salivary gland development and adult progenitor cell function, are expressed by cells in this cluster (33, 59). These changes are consistent with a primarily a homeostatic function for group 3 fibroblasts. In summary, examination of differential gene expression highlights the likely functional separation of the fibroblast subgroups.
3.3 Fibroblast populations change in abundance and transcriptomes over injury time course
The abundance of cells in each of the fibroblast subgroups changes over the injury time course (Figure 4A,B). Group 1 fluctuates slightly but is largely stable in abundance over the time course, while group 3 declines in abundance over the time course. Group 2, however, expands dramatically in response to injury. To determine how the gene expression of the fibroblast groups changed in response to injury, we compared the three groups at the three surgical timepoints (e.g., 1, 3, and 7 days post-injury) to the homeostatic timepoint to identify the genes that were differentially regulated.

Figure 4. Fibroblast groups exhibit altered transcriptomes at different timepoints. (A) UMAPs highlighting the cells from each timepoint and where they are located within the Seurat object. (B) Line graph depicts how the populations present in each group change over the course of the injury. (C) Heatmaps show the 50 top differentially expressed genes in each group separated by timepoint. All genes shown are upregulated when compared to the homeostatic population in the group. Order from top to bottom is 1-, 3- and 7-days post-ligation. Order from left to right is Group 1, Group 2, Group 3. Refer to Supplementary Table S4 for complete list of differentially expressed genes separated by groups and timepoint.
At day 1 after ligation, all 3 groups exhibit an upregulation of heat shock protein transcripts and cytoplasmic translation, indicative of an acute response to injury (Figure 4C, Supplementary Table S4). There is ubiquitous expression of the chaperone, Calreticulin (Calr), in all 1-day ligated fibroblasts. At day 3 after ligation, group 1 has a translation signature and includes selective expression of transcripts such as Fibronectin (Fn1) along with Spp1, whereas groups 2 and 3 exhibit a more robust fibrogenic signature with large increases in expression of fibrillar ECM genes encoding molecules such as Collagens 1, 3, and 5. We also see the expression of Postn and Cthrc1 which persists at 7-days in group 2. Group 3 includes Fn1, Col1a1, Col3a1, Col5a1, Thbs1, and Transforming Growth Factor Beta Receptor 2 (Tgfbr2). At day 7 after ligation, groups 1 and 3 exhibit a cell migration and morphogenesis signature, suggestive of resolution of acute injury response, while group 2 retains a strong fibrogenic signature. These data indicate that group 2 is the primary fibrogenic population in the injury response.
3.4 Dynamics of the inflammatory and matrisome components of the fibrotic injury time course
We then investigated the inflammatory and matrisome transcript dynamics in the fibroblast cells in the fibrotic injury time course. Using Gene Set Enrichment Analysis (GSEA) gene lists curated by Molecular Signatures Database (MSigDB), matrisome components and the hallmark inflammatory response genes (2, 34), we quantified the average number of transcripts that were expressed in these categories over the time course to get a glimpse into the order in which tissue remodeling occurs (Figure 5A). Proteoglycans are expressed at the highest levels at the homeostatic timepoint. 1 day after ligation, we see an increase in inflammatory gene expression and ECM regulators. At 3 days, we see the highest expression of collagen genes by total transcripts, and glycoprotein expression peaks 7 days after ligation. As the diversity of genes expressed in these categories was also of interest, we used a dotplot to examine the number of different genes expressed in the categories (Figure 5B). We see similar patterns to gene abundance with an enriched diversity of inflammatory genes present at 1 day post ligation and secreted factors at 3 days post ligation, but also some differences, such as collagen transcript diversity being the highest 7 days post ligation. ECM regulators also show some difference from total genes expressed with the highest diversity of ECM regulator genes present 3 days after ligation.

Figure 5. Matrisome gene expression shows changes in abundance and diversity by timepoint. (A) A violin plot shows the total transcripts expressed by cells at the different timepoints for gene lists made from the matrisome master list. The gene categories are on the right Y axis. Cells are removed from the plot for easier interpretation. (B) A dotplot shows the diversity of genes for the same lists.
We then looked at the individual transcripts related to the matrisome in the time course data to identify contributions of each fibroblast group at each timepoint (Supplementary Table S5). Using this approach, we can see not only which group contributes the most genes to each category, but also at which timepoint these genes are upregulated or downregulated. To normalize for differences in cell numbers, we used transcripts per cell, which is the number of transcripts divided by the number of cells present in each population.
When we examined the collagen genes, we found that fibrillar collagens (e.g., Col1a1, Col1a2, and Col 3a1) are the main contributors to transcript abundance showing 3-to-4-fold increases relative to the homeostatic population at day 3 and persisting through day 7, with all three genes expressed most abundantly by group 2 cells (Figure 6A). Glycoprotein analysis showed that Matrix Gla Protein (Mgp), Spp1, and Insulin-like Growth Factor Binding Protein 7 (Igfbp7) contribute the most to glycoprotein transcript abundance, and glycoproteins are dominant in group 1 (Figure 6B). Their expression increases steadily beginning at day 1 and continuing to day 7. Compared to the homeostatic state, we detected a 4-fold increase in Mgp expression, a 452-fold increase in Spp1 expression, and a 2-fold increase in Igfbp7 in group 1 at 7 days post injury. For inflammatory transcripts, Tissue Inhibitor of Metalloproteinases 1 (Timp1), and C-X-C Motif Chemokine ligands 5 and 10 (Cxcl5 and Cxcl10) were the top contributors. They were most highly expressed at day 1, and these transcripts were not dominated by a single fibroblast group (Figure 6C). Proteoglycans were the only category to show a decrease in abundance after injury, with Decorin (Dcn), Biglycan (Bgn), and Lumican (Lum) having the highest transcript numbers, and these transcripts were expressed by all fibroblast groups (Figure 6D). Examination of ECM regulators showed that Matrix Metalloproteinase 3 (Mmp3), Cathepsin L (Ctsl), and Serpin Family F Member 1 (Serpinf1) have the highest transcript numbers, with only Mmp3 being injury-responsive and predominantly expressed by group 1 (Figure 6E). Examination of secreted factors showed that Cxcl14, Follistatin-Related Protein 1 (Fstl1), and C-C Motif Chemokine 11 (Ccl11) contributed the most to the transcript abundance, but these transcripts were not dominated by a single fibroblast group (Figure 6F).

Figure 6. Most abundantly expressed matrisome genes in fibroblast groups 7 days after ligation injury. (A) The three most abundantly expressed collagens are shown with the number of transcripts expressed per cell. Each graph is divided into the 4 timepoints. (B) Top 3 most expressed glycoproteins. (C) Top 3 most expressed inflammation genes. (D) Top 3 most expressed proteoglycans. (E) Top 3 most expressed ECM regulators. (F) Top 3 most expressed secreted ECM proteins. Refer to Supplementary Table S5 for data used to generate graphs.
This data reveals a strong fibrogenic injury response in the gland ligation time course, with complex dynamic patterns of these important matrisome genes in the fibroblast subgroups. Group 1 fibroblasts are the main contributors to glycoprotein and enzyme secretion, indicating that this population is the main signaling population present during an injury response. Group 2 fibroblasts are the main contributors of fibrillar proteins making them the main contributors to the fibrogenic phenotype observed. Group 3 fibroblasts may serve primarily a homeostatic function.
4 Discussion
The mechanisms regulating a productive, healing fibrotic response and how they differ from pathological fibrosis is topic of considerable current interest for development of improved therapeutics for fibrotic pathologies. Using the reversable ductal ligation injury model in young female mouse salivary glands, we describe a time course of the transcriptomic signatures of this productive fibrotic response. As there are few studies that report transcriptomic changes over time, our time course transcriptomic data will serve as a resource for comparison of the transcriptomic signatures of this productive fibrotic response with resolution of the fibrotic response and with pathological fibrotic responses. Similar to the well-studied wound healing response of skin, we observe an early inflammatory response together with a developing fibrotic response that dominates the fibrotic phase of this reversable injury. We observe a progressive increase in lymphocyte populations such as the macrophage/monocyte and T-cell populations, which peak 3 days after injury. Importantly, the injury-responsive group 2 fibroblasts are the primary fibrotic cells, and these cells show progressive production of ECM proteins from day 3 to day 7.
Different organs and different injuries have distinct fibroblast trajectories and mediators; however, we do detect many universal fibroblast markers (5) expressed by the fibroblasts in our transcriptome dataset, indicating that fibroblasts within salivary glands are similar to other described organs and highlighting the potential for conserved therapeutic approaches. An example of this can be found in group 2 fibroblasts which express Secreted Protein Acidic and Cysteine Rich (Sparc), and Col3a1 which has been reported to be expressed by fibroblasts found in a perivascular niche in multiple inflamed human tissue samples including salivary glands of humans with Sjogren's Disease (35). Notably, group 1 fibroblasts have expression of genes associated with vasculature development and yet group 2 fibroblasts express Sparc and Col3a1 which are associated with perivascular fibroblasts. Perhaps these two fibroblast populations function synergistically to promote vasculature development. Related to Korsunsky et al, fibroblasts associated with T lymphocyte niches that express Cxcl10 and Ccl19 were not found in our dataset, which is not surprising, as this is a simple injury model.
Here, we continued previous work which highlighted the transcriptomic differences in salivary gland fibroblasts 14 days after ligation (18). At 14 days after injury, we detected upregulation of similar genes that we here detect increased 7 days after injury. These changes observed in our groups are indicative that the earlier changes persist throughout the injury. Because of the larger numbers of cells in this current study, we can identify distinct fibroblasts subpopulations which seem to have distinct functions. Of note, these different fibroblast subgroups may communicate as ligands and cognate receptors can be found expressed by different fibroblast groups as is the case with Tgfb1, which is expressed by group 1, and Tgfbr2, which is expressed by group 3. This regulation may enable fibroblasts to dynamically respond to ECM changes allowing secretion of proteins as the need arises.
Our 1-day timepoint sheds light on the initial inflammatory response seen after damage to tissue. Earlier in the injury, we see expression of chemoattractants like Csf1, which increase at day 1 but are reduced as the injury progress, we also see widespread expression of heat shock proteins which serve as chaperones, assisting in protein folding but also have functions in activating degradation of aberrant proteins via the endoplasmic reticulum-associated protein degradation pathway (36, 37).
The abundance of fibroblasts allows us to not only see the changes over time, but also the unique subpopulations that are present and the genes that lend to their uniqueness. Genes expressed by group 1 include Piezo2 which is a mechanosensitive ion channel that responds to matrix rigidity (38), a known hallmark of fibrosis (39, 41) and Pdgfrb which mediates PDGF-BB-induced proliferation and migration (42–44). Cthrc1 expression has been used to identify fibroblasts producing elevated levels of fibrillar ECM proteins (6) which can be seen in group 2 fibroblasts. These same cells express Sfrp2 and Dpt proteins which have been found in fibroblast progenitor populations that can give rise to myofibroblasts (5, 45). These findings point to group 2 fibroblasts serving as the fibrogenic population in this model.
It is interesting to see the changes in expression of ECM proteins at different timepoints after surgical injury. Here we see that not only are genes upregulated in response to injury, as is the case with the glycoproteins, but we see a large decreases in Dcn. Dcn is of interest due to its interaction with collagen fibrils and ability to regulate the diameter of collagen fibrils (46–48). It is possible that the decrease in Dcn expression leads to the altered collagen organization detected after ductal ligation injury in previous studies (17, 18). Another notable change is in the group 1 fibroblasts where we see a large increase in Spp1. Spp1 has been implicated in attraction of lymphocytes that infiltrate injured tissue (49). Spp1 is also a potent mediator of cell communication that is dysregulated in many diseases and can promote a myofibroblast transition leading to expression of Periostin (Postn) (50). Interestingly, Postn is expressed by group 2 fibroblasts suggesting possible signaling from group 1 to group 2.
Spp1 steadily increases and may mediate communication between cells at later stages as the early response cytokines decrease in abundance. In-vitro assays performed using Spp1 knock-down and overexpression showed that Spp1 is a potent mediator of fibroblast motility and myofibroblast differentiation (40). We also see the expression of other cytokines such as Tgfb1 and Il34. Il34 is known to promote inflammation via the Colony Stimulating Factor receptor (Csfr), suggesting that Il34 may be a mediator of the inflammation observed in the ductal ligation injury (51). Tgfb1 is a pleiotropic driver of fibrosis and a potent mediator of tissue repair (52). We will look to see the role that these genes play in a more 3D context in mouse salivary gland organoids where we can add genetically manipulated salivary gland fibroblasts, as in our prior work (53).
We see progressive increases of collagens, such as Col1a1, Col1a2, Col3a1, Col5a2, and Col8a1 starting at day 3, which succeeds the expression of Spp1 at day 1. Interestingly, Collagen VIII has been found to be expressed in tumor microenvironments and higher levels of expression results in a lower survival for patients (54). We should also note the differences in ratios of fibrillar collagens type I, II, and V as compared to basement membrane collagen type IV. It has been shown that the relative amount of these collagens can regulate cell polarity and function (55). With salivary gland tissue fibrosis, we see a decrease in saliva production, which may be caused by epithelial cell polarity being disrupted in response to a disrupted basement membrane. While the differential gene expression predicts unique functions of fibroblast subclusters, functional studies are required to validate our findings. Investigating the protein levels and disruption of the genes increased in fibroblast groups 1 and 2 will reveal the genes driving tissue repair that may be altered leading to diseased states.
As this acute experimental injury model is reversible, it will be interesting to evaluate the mechanisms involved in fibrosis remediation in future work. The question of whether a fibrotic response in a patient is reversible or irreversible has great clinical importance. Thus, it will be important to compare this reversible response to irreversible fibrotic responses, as that which occurs with partial gland resection in our prior work (56) and with the long-term fibrotic response that occurs following irradiation for head and neck cancers. As fibrosis is said to be responsible for 45% of all deaths in the industrialized world (57), there is a need to recognize fibrotic responses as irreversible or potentially reversible and to develop effective therapeutics to reverse fibrosis in salivary glands and other organs.
Development of effective therapeutics for fibrotic disease has been difficult and often disappointing. Interestingly, a recent report investigating novel therapeutic targets in idiopathic pulmonary fibrosis (IPF) identified Leukemia Inhibitory Factor receptor (LIFR) as an “autocrine master amplifier” of multiple upstream activators of lung fibroblasts (32). TGFβ1, IL-4, and IL-13 stimulation of fibroblasts required the LIF-LIFR axis to evoke a strong fibrogenic effector response in fibroblasts. Thus, in IPF, LIFR drives an autocrine circuit that amplifies and sustains pathogenic activation of IPF fibroblasts. They suggest that targeting a single, downstream master amplifier of fibroblasts, like LIFR, is an attractive alternative therapeutic strategy that can attenuate the profibrotic effects of multiple upstream stimuli. LIFR is expressed primarily in the group 3 fibroblasts in our dataset, which is not the primary ECM-producing group; however, LIF is strongly upregulated in in the group 3 fibroblasts in response to the ligation injury in these young female mice. Communication between the group 3 fibroblasts and the group 2 fibroblasts may be an important component of the fibrotic response. Future studies will be able to investigate whether the dysregulation of the LIF/LIRF signaling pathway is a driver of pathological fibrosis in salivary gland fibrotic diseases and if this is a universal amplification circuit that may be amenable to therapeutic modulation.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The animal study was approved by Institutional Animal Care and Use Committee. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
JT: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. JK: Data Curation, Software, Writing – review & editing. SG: Investigation, Writing – review & editing. DN: Writing – original draft, Writing – review & editing. ML: Funding Acquisition, Resources, Supervision, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. Research reported in this publication was supported by an NIH Training Grant (T32GM132066, RNA Science and Technology in Health and Disease) to JT and the National Institute of Dental & Craniofacial Research of the National Institutes of Health under Award Numbers R01DE027953, R01DE030626 and R56DE033253 to ML and by the Williams-Raycheff Endowment maintained by the University at Albany Foundation. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the Williams-Raycheff Endowment.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fdmed.2025.1581376/full#supplementary-material
Supplementary Figure 1 | UMAP shows all integrated replicates. (A) A split UMAP shows the 12 replicates that were integrated into the Seurat object. Cell are grouped based on assigned identities.
Supplementary Figure 2 | Differentially expressed genes in fibroblast clusters. (A) A heatmap shows the top 5 differentially expressed genes in the fibroblast clusters.
References
1. Younesi FS, Miller AE, Barker TH, Rossi FMV, Hinz B. Fibroblast and myofibroblast activation in normal tissue repair and fibrosis. Nat Rev Mol Cell Biol. (2024) 25:617–38. doi: 10.1038/s41580-024-00716-0
2. Naba A, Clauser KR, Hoersch S, Liu H, Carr SA, Hynes RO. The matrisome: in silico definition and in vivo characterization by proteomics of normal and tumor extracellular matrices*. Mol Cell Proteomics. (2012) 11:M111.014647. doi: 10.1074/mcp.M111.014647
3. Micallef L, Vedrenne N, Billet F, Coulomb B, Darby IA, Desmoulière A. The myofibroblast, multiple origins for major roles in normal and pathological tissue repair. Fibrogenesis Tissue Repair. (2012) 5:S5. doi: 10.1186/1755-1536-5-S1-S5
4. Pakshir P, Noskovicova N, Lodyga M, Son DO, Schuster R, Goodwin A, et al. The myofibroblast at a glance. J Cell Sci. (2020) 133:jcs227900. doi: 10.1242/jcs.227900
5. Buechler MB, Pradhan RN, Krishnamurty AT, Cox C, Calviello AK, Wang AW, et al. Cross-tissue organization of the fibroblast lineage. Nature. (2021) 593:575–9. doi: 10.1038/s41586-021-03549-5
6. Tsukui T, Sun K-H, Wetter JB, Wilson-Kanamori JR, Hazelwood LA, Henderson NC, et al. Collagen-producing lung cell atlas identifies multiple subsets with distinct localization and relevance to fibrosis. Nat Commun. (2020) 11:1920. doi: 10.1038/s41467-020-15647-5
7. Deng C-C, Hu Y-F, Zhu D-H, Cheng Q, Gu J-J, Feng Q-L, et al. Single-cell RNA-seq reveals fibroblast heterogeneity and increased mesenchymal fibroblasts in human fibrotic skin diseases. Nat Commun. (2021) 12:3709. doi: 10.1038/s41467-021-24110-y
8. Martin P, Pardo-Pastor C, Jenkins RG, Rosenblatt J. Imperfect wound healing sets the stage for chronic diseases. Science. (2024) 386:eadp2974. doi: 10.1126/science.adp2974
9. Nelson DA, Kazanjian I, Melendez JA, Larsen M. Senescence and fibrosis in salivary gland aging and disease. J Oral Biol Craniofac Res. (2024) 14:231–7. doi: 10.1016/j.jobcr.2024.02.009
10. Ma D, Feng Y, Lin X. Immune and non-immune mediators in the fibrosis pathogenesis of salivary gland in Sjögren’s syndrome. Front Immunol. (2024) 15:1421436. doi: 10.3389/fimmu.2024.1421436
11. Jasmer KJ, Gilman KE, Muñoz Forti K, Weisman GA, Limesand KH. Radiation-induced salivary gland dysfunction: mechanisms, therapeutics and future directions. J Clin Med. (2020) 9:4095. doi: 10.3390/jcm9124095
12. Lombaert I, Movahednia MM, Adine C, Ferreira JN. Concise review: salivary gland regeneration: therapeutic approaches from stem cells to tissue organoids. Stem Cells. (2017) 35:97–105. doi: 10.1002/stem.2455
13. Chibly AM, Nguyen T, Limesand KH. Palliative care for salivary gland dysfunction highlights the need for regenerative therapies: a review on radiation and salivary gland stem cells. J Palliat Care Med. (2014) 4:1000180. doi: 10.1371/journal.pone.0107893
14. Yin H, Pranzatelli TJF, French BN, Zhang N, Warner BM, Chiorini JA, NIDCD/NIDCR Genomics and Computational Biology Core. Sclerosing sialadenitis is associated with salivary gland hypofunction and a unique gene expression profile in Sjögren’s syndrome. Front Immunol. (2021) 12:699722. doi: 10.3389/fimmu.2021.699722
15. Takashi I, Kanai R, Hasegawa K, Ogaeri T, Tran SD, Sumita Y. Recent progress in regenerative therapy for damaged salivary glands: from bench to bedside. Oral Dis. (2024) 30:38–49. doi: 10.1111/odi.14692
16. Cotroneo E, Proctor GB, Paterson KL, Carpenter GH. Early markers of regeneration following ductal ligation in rat submandibular gland. Cell Tissue Res. (2008) 332:227–35. doi: 10.1007/s00441-008-0588-6
17. Woods LT, Camden JM, El-Sayed FG, Khalafalla MG, Petris MJ, Erb L, et al. Increased expression of TGF-β signaling components in a mouse model of fibrosis induced by submandibular gland duct ligation. PLoS One. (2015) 10:e0123641. doi: 10.1371/journal.pone.0123641
18. Altrieth AL, O’Keefe KJ, Gellatly VA, Tavarez JR, Feminella SM, Moskwa NL, et al. Identifying fibrogenic cells following salivary gland obstructive injury. Front Cell Dev Biol. (2023) 11:1190386. doi: 10.3389/fcell.2023.1190386
19. Altrieth AL, Kenney J, Nelson DA, Suarez EG, Gellatly V, Gabunia S, et al. Single-cell transcriptomic analysis of salivary gland endothelial cells. J Dent Res. (2024) 103:269–78. doi: 10.1177/00220345231219987
20. Martinez JR, Bylund DB, Cassity N. Progressive secretory dysfunction in the rat submandibular gland after excretory duct ligation. Arch Oral Biol. (1982) 27:443–50. doi: 10.1016/0003-9969(82)90082-6
21. Kimura K, Ramirez K, Nguyen TAV, Yamashiro Y, Sada A, Yanagisawa H. Contribution of PDGFRα-positive cells in maintenance and injury responses in mouse large vessels. Sci Rep. (2021) 11:8683. doi: 10.1038/s41598-021-88126-6
22. Madisen L, Zwingman TA, Sunkin SM, Oh SW, Zariwala HA, Gu H, et al. A robust and high-throughput cre reporting and characterization system for the whole mouse brain. Nat Neurosci. (2010) 13:133–40. doi: 10.1038/nn.2467
23. Fleming SJ, Chaffin MD, Arduini A, Akkad A-D, Banks E, Marioni JC, et al. Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nat Methods. (2023) 20:1323–35. doi: 10.1038/s41592-023-01943-7
24. R Project for Statistical Computing. R: the R project for statistical computing (2021). Available at: https://www.r-project.org/ (Accessed June 1, 24).
25. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive integration of single-cell data. Cell. (2019) 177:1888–902.e21. doi: 10.1016/j.cell.2019.05.031
26. McGinnis CS, Murrow LM, Gartner ZJ. Doubletfinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. (2019) 8:329–37.e4. doi: 10.1016/j.cels.2019.03.003
27. Germain P-L, Lun A, Meixide CG, Macnair W, Robinson MD. Doublet identification in single-cell sequencing data using scDblFinder (2022).
28. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. (2019) 20:296. doi: 10.1186/s13059-019-1874-1
29. Seurat-Guided Clustering Tutorial. Satijalab. (2023). Available at: https://satijalab.org/seurat/articles/pbmc3k_tutorial.html (Accessed August 1, 2024).
30. Puente A, Fortea JI, Cabezas J, Arias Loste MT, Iruzubieta P, Llerena S, et al. LOXL2-A new target in antifibrogenic therapy? Int J Mol Sci. (2019) 20:1634. doi: 10.3390/ijms20071634
31. Yang N, Cao D-F, Yin X-X, Zhou H-H, Mao X-Y. Lysyl oxidases: emerging biomarkers and therapeutic targets for various diseases. Biomed Pharmacother. (2020) 131:110791. doi: 10.1016/j.biopha.2020.110791
32. Nguyen HN, Jeong Y, Kim Y, Kamiya M, Kim Y, Athar H, et al. Leukemia inhibitory factor (LIF) receptor amplifies pathogenic activation of fibroblasts in lung fibrosis. Proc Natl Acad Sci U S A. (2024) 121:e2401899121. doi: 10.1073/pnas.2401899121
33. Mauduit O, Aure MH, Delcroix V, Basova L, Srivastava A, Umazume T, et al. A mesenchymal to epithelial switch in Fgf10 expression specifies an evolutionary-conserved population of ionocytes in salivary glands. Cell Rep. (2022) 39:110663. doi: 10.1016/j.celrep.2022.110663
34. Howe DG, Blake JA, Bradford YM, Bult CJ, Calvi BR, Engel SR, et al. Model organism data evolving in support of translational medicine. Lab Anim. (2018) 47:277–89. doi: 10.1038/s41684-018-0150-4
35. Korsunsky I, Wei K, Pohin M, Kim EY, Barone F, Major T, et al. Cross-tissue, single-cell stromal atlas identifies shared pathological fibroblast phenotypes in four chronic inflammatory diseases. Med. (2022) 3:481–518.e14. doi: 10.1016/j.medj.2022.05.002
36. Stricher F, Macri C, Ruff M, Muller S. HSPA8/HSC70 chaperone protein: structure, function, and chemical targeting. Autophagy. (2013) 9:1937–54. doi: 10.4161/auto.26448
37. Rosenzweig R, Nillegoda NB, Mayer MP, Bukau B. The Hsp70 chaperone network. Nat Rev Mol Cell Biol. (2019) 20:665–80. doi: 10.1038/s41580-019-0133-3
38. Xiao B. Mechanisms of mechanotransduction and physiological roles of PIEZO channels. Nat Rev Mol Cell Biol. (2024) 25:886–903. doi: 10.1038/s41580-024-00773-5
39. Lu P, Takai K, Weaver VM, Werb Z. Extracellular matrix degradation and remodeling in development and disease. Cold Spring Harb Perspect Biol. (2011) 3:a005058. doi: 10.1101/cshperspect.a005058
40. Ding H, Xu Z, Lu Y, Yuan Q, Li J, Sun Q. Kidney fibrosis molecular mechanisms Spp1 influences fibroblast activity through transforming growth factor beta smad signaling. iScience. (2024a) 27:109839. doi: 10.1016/j.isci.2024.109839
41. Ding J-F, Tu B, Song K, Liu Z-Y, Lin L-C, Liu Z-Y, et al. Epitranscriptomic regulation of cardiac fibrosis via YTHDF1-dependent PIEZO2 mRNA m6A modification. Cardiovasc Res. (2024b) 120(17):2236–48. doi: 10.1093/cvr/cvae239
42. Gao Z, Sasaoka T, Fujimori T, Oya T, Ishii Y, Sabit H, et al. Deletion of the PDGFR-β gene affects key fibroblast functions important for wound healing*. J Biol Chem. (2005) 280:9375–89. doi: 10.1074/jbc.M413081200
43. Buhl EM, Djudjaj S, Klinkhammer BM, Ermert K, Puelles VG, Lindenmeyer MT, et al. Dysregulated mesenchymal PDGFR-β drives kidney fibrosis. EMBO Mol Med. (2020) 12:e11021. doi: 10.15252/emmm.201911021
44. Yamamoto S, Fukumoto E, Yoshizaki K, Iwamoto T, Yamada A, Tanaka K, et al. Platelet-derived growth factor receptor regulates salivary gland morphogenesis via fibroblast growth factor expression*. J Biol Chem. (2008) 283:23139–49. doi: 10.1074/jbc.M710308200
45. Tabib T, Huang M, Morse N, Papazoglou A, Behera R, Jia M, et al. Myofibroblast transcriptome indicates SFRP2hi fibroblast progenitors in systemic sclerosis skin. Nat Commun. (2021) 12:4384. doi: 10.1038/s41467-021-24607-6
46. Reese SP, Underwood CJ, Weiss JA. Effects of decorin proteoglycan on fibrillogenesis, ultrastructure, and mechanics of type I collagen gels. Matrix Biol. (2013) 32(7–8):414–23. doi: 10.1016/j.matbio.2013.04.004
47. Dong Y, Zhong J, Dong L. The role of decorin in autoimmune and inflammatory diseases. J Immunol Res. (2022) 2022:1283383. doi: 10.1155/2022/1283383
48. Weber IT, Harrison RW, Iozzo RV. Model structure of decorin and implications for collagen fibrillogenesis *. J Biol Chem. (1996) 271:31767–70. doi: 10.1074/jbc.271.50.31767
49. Lund SA, Giachelli CM, Scatena M. The role of osteopontin in inflammatory processes. J Cell Commun Signal. (2009) 3:311–22. doi: 10.1007/s12079-009-0068-0
50. Hoeft K, Schaefer GJL, Kim H, Schumacher D, Bleckwehl T, Long Q, et al. Platelet-instructed SPP1+ macrophages drive myofibroblast activation in fibrosis in a CXCL4-dependent manner. Cell Rep. (2023) 42:112131. doi: 10.1016/j.celrep.2023.112131
51. Monteleone G, Franzè E, Troncone E, Maresca C, Marafini I. Interleukin-34 mediates cross-talk between stromal cells and immune cells in the gut. Front Immunol. (2022) 13:873332. doi: 10.3389/fimmu.2022.873332
52. Kim KK, Sheppard D, Chapman HA. TGF-β1 signaling and tissue fibrosis. Cold Spring Harb Perspect Biol. (2018) 10:a022293. doi: 10.1101/cshperspect.a022293
53. Hosseini ZF, Nelson DA, Moskwa N, Sfakis LM, Castracane J, Larsen M. FGF2-dependent mesenchyme and laminin-111 are niche factors in salivary gland organoids. J Cell Sci. (2018) 131:jcs208728. doi: 10.1242/jcs.208728
54. Zhang J, Liu J, Zhang H, Wang J, Hua H, Jiang Y. The role of network-forming collagens in cancer progression. Int J Cancer. (2022) 151:833–42. doi: 10.1002/ijc.34004
55. Karsdal MA, Nielsen SH, Leeming DJ, Langholm LL, Nielsen MJ, Manon-Jensen T, et al. The good and the bad collagens of fibrosis – their role in signaling and organ function. Adv Drug Delivery Rev. (2017) 121:43–56. doi: 10.1016/j.addr.2017.07.014
56. O’Keefe KJ, DeSantis KA, Altrieth AL, Nelson DA, Taroc EZM, Stabell AR, et al. Regional differences following partial salivary gland resection. J Dent Res. (2020) 99:79–88. doi: 10.1177/0022034519889026
57. Henderson NC, Rieder F, Wynn TA. Fibrosis: from mechanisms to medicines. Nature. (2020) 587:555–66. doi: 10.1038/s41586-020-2938-9
58. Gao Y, Li J, Cheng W, Diao T, Liu H, Bo Y, et al. Cross-tissue human fibroblast atlas reveals myofibroblast subtypes with distinct roles in immune modulation. Cancer Cell. (2024) 42(10):1764–83.e10. doi: 10.1016/j.ccell.2024.08.020
Keywords: salivary glands, extracellular matrix (ECM), fibroblasts, injury response, inflammation, single cell RNA-seq
Citation: Tavarez JR, Kenney J, Gabunia S, Nelson DA and Larsen M (2025) Temporal evolution of fibroblast responses following salivary gland ductal ligation injury. Front. Dent. Med. 6:1581376. doi: 10.3389/fdmed.2025.1581376
Received: 22 February 2025; Accepted: 15 April 2025;
Published: 1 May 2025.
Edited by:
Bruno Špiljak, University of Zagreb, CroatiaReviewed by:
Blake Matthew Warner, National Institutes of Health (NIH), United StatesDanielle Wu, University of Texas Health Science Center at Houston, United States
Wenpeng Song, Capital Medical University, China
Copyright: © 2025 Tavarez, Kenney, Gabunia, Nelson and Larsen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Melinda Larsen, bWxhcnNlbkBhbGJhbnkuZWR1