Introduction

Extracellular matrix (ECM)-producing myofibroblasts are a key therapeutic target to combat tissue fibrosis, one of the biggest unresolved clinical problems across most major chronic diseases. Several recent single-cell RNA-sequencing (scRNAseq) studies described distinct subsets of collagen-producing stromal cells in mouse and human lungs with distinct spatial locations and different functions in supporting epithelial repair [1–5]. How this heterogeneity of fibroblast identities leads to different cell states and functions in fibrotic disease is unclear.

Genetic lineage tracing in mouse models provides evidence for an alveolar lipofibroblast-to-myofibroblast switch after lung injury that is reversible during the resolution of transient fibrosis upon completed epithelial regeneration [6, 7]. The heterogeneity of cellular sources for myofibroblasts in lung fibrosis remains the subject of current investigations. Apart from lipofibroblasts, alveolar pericytes are also potential sources of myofibroblasts [8]. Proliferation, evasion of apoptosis and invasive capacity of myofibroblasts are key hallmarks of fibrotic disease, and transforming growth factor (TGF)β is a known master regulator of these processes [9, 10]. A vast amount of literature demonstrates the TGFβ-induced induction of α-smooth muscle actin (ACTA)2 in fibroblasts, the most widely used marker of mature myofibroblasts in tissues. However, fate mapping and immunofluorescence analysis of fibrotic tissues and myofibroblast foci also show that a substantial fraction of fibroblasts is ACTA2− [11, 12], suggesting additional molecular complexity and heterogeneity among injury-activated fibroblasts. A recent single-cell analysis of collagen-producing cells in lung fibrosis revealed CTHRC1 as a specific marker of highly invasive ACTA2+ myofibroblasts. These cells occur in both mouse and human lung and occupy fibroblastic foci in idiopathic pulmonary fibrosis (IPF) [1].

The regulation of injury-activated fibroblast states is not well understood and of high clinical relevance. Our study reveals the spatiotemporal evolution of distinct fibroblast states towards Cthrc1+ myofibroblasts, highlighting early events of injury-induced fibroblast activation that preceded myofibroblastic differentiation. We discovered a novel SFRP1+/ACTA2− transitional state that is initially noninvasive and only becomes invasive upon TGFβ-driven differentiation towards the CTHRC1+ myofibroblast state. We show that SFRP1 modulates TGFβ1-induced fibroblast invasion and RHOA pathway activity, which constitutes a novel pathway with potential for targeting fibrotic disease mediated by myofibroblasts.

Material and methods

For details of the materials and methods, please refer to the supplementary material.

Experimental study design

C57B6 (males or females) mice were treated with 3.5 units·kg−1 body mass of bleomycin administered intratracheally. Up to two bleomycin-instilled mice were sacrificed at time points (days 0, 2–14, 21, 28, 35, 56) after instillation and subject to scRNAseq using the Drop-seq platform. Tcf21-lineage labelled (Tcf21m-Crem-R26R-tdTomato) mesenchymal cells were analysed using the 10x scRNAseq platform. Two healthy C57B6 mouse lungs were used for SCRINSHOT. Primary human lung fibroblasts (pHLFs) and micro-computed tomography (CT) staged IPF tissues, and data from an integrated IPF cell atlas were used for human validation.

Human tissue and ethics statement

pHLFs of non-chronic lung disease donors were obtained from the CPC-M bioArchive at the Comprehensive Pneumology Center (Munich, Germany). Participants provided written informed consent to participate in this study, in accordance with approval by the local ethics committee of the Ludwig Maximilians University (Germany) (project 333-10). Micro-CT-staged IPF samples were provided from the KU Leuven lung biobank (ethical approval S52174). Samples were derived from explanted lungs, after written informed consent from all patients. Unused donor lungs were included as controls, following Belgian legislation. Three IPF lungs and three controls were included.

Data availability

RNA-seq data were deposited to the Gene Expression Omnibus (GEO) database. The high temporal mesenchymal enriched Drop-seq data can be found with the accession code GSE207851, and the Tcf21-lineage labelled mesenchymal 10x data with the accession code GSE207687. Microarray data of pHLFs (Sfrp1-siRNA knockdown) can be found with the accession code GSE207561.

ResultsHeterogeneity of mesenchymal cells at distinct spatial localisations in the lung

The Tcf21+ lineage constitutes the lung lipofibroblast population [13]. To study heterogeneity within the Tcf21+ and Tcf21− fibroblast lineages, we tamoxifen-labelled Tcf21+ cells in lungs at 11 weeks of age using Tcf21m-Crem-R26R-tdTomato mice (supplementary figure S1a). Tcf21 lineage-positive (lin+) cells were indeed found near alveolar type 2 (AT2) cells and alveolar capillaries (supplementary figure S1b). We flow-sorted Tcf21 lin+ and lin− stromal cells and performed scRNAseq (n=4 mice, k=12.068 cells), identifying six distinct cell types with different Tcf21-lineage proportions (figure 1a and supplementary figure S1c and d). Previously established marker genes [1, 14, 15] were used to annotate clusters (supplementary table S1). All cell type identities expressed type 1 collagen (Col1a2) (supplementary figure S2c) and were consistent with previous work [1].

FIGURE 1

Stromal cell heterogeneity in the adult lung is associated with distinct spatial locations. a) Uniform manifold approximation and projection (UMAP) clustering depicts the seven distinct mesenchymal cell types (12 068 cells) identified in the healthy adult mouse lung. Circles mark Pdgfra- and Pdgfrb-expressing cells. b) The seven mesenchymal cell types are characterised with distinct gene expression profiles, as depicted in the matrixplot. c) Marker genes of the indicated mesenchymal subtypes are enriched for characteristic gene ontology (GO):biological process (BP) terms. False discovery rate <10%. Top 500 marker genes considered. d) Single mRNA multiplexed fluorescence in situ hybridisation (SCRINSHOT) with two cell type specific marker genes per cell type, as well as general markers allowed to identify preferential spatial locations of the mesenchymal cell types in the mouse lung. Scale bars=50 μm (overviews); 10 µm (insets). e) Localisation analysis shows stereotypic localisation of the mesenchymal cell types to peribronchiolar, adventitial and alveolar regions. The spatial location of 976 cells was quantified over five regions with 508 pericytes, 70 peribronchiolar cells, 206 lipofibroblasts, 41 adventitial fibroblasts and 151 smooth muscle cells (SMCs) counted. Graphs show cell density as percentage of cell type specific cells compared to all the cells (4′,6-diamidino-2-phenylindole (DAPI) stain-positive nuclei) in the area. Error bars represent sem. f) Colocalisation analysis of 1269 alveolar type 2 (AT2) cells revealed preferred interaction partners among the mesenchymal cell types. A total of 291 AT2 interactions with mesenchymal cell types in the alveolar space were found: 167 pericytes; 93 lipofibroblasts, 15 peribronchiolar fibroblasts; three adventitial fibroblasts; 13 SMCs. g) SCRINSHOT pictures exemplifying the colocalisation of AT2 cells with distinct mesenchymal subtypes. Scale bars=10 μm. ECM: extracellular matrix.

We identified three main Pdgfra+ populations (lipofibroblasts, adventitial fibroblasts, Adh7+ fibroblasts), two Pdgfrb+ populations (smooth muscle cells and pericytes) with highly distinct marker genes and pathway enrichments (figure 1b and c), and Pdgfra/Pdgfrb double-positive cells, as recently found in the kidney [16]. These cells expressed Lgr5 and Lgr6 (figure 1b), markers for a peribronchiolar fibroblast population [4]. Importantly, a recent study identified similar LGR5+ fibroblast populations in distal human airways [17]. Subclustering of this Hhip+ peribronchiolar fibroblast population revealed additional complexity with Lgr5/Lgr6 single- and double-positive populations (supplementary figure S1e–g). All cell types were lineage-labelled in the Tcf21m-Crem-R26R-tdTomato mouse, with exception of the Hhip+ peribronchiolar fibroblasts (supplementary figure S1c and d). Thus, Tcf21lin− peribronchiolar fibroblasts may constitute a developmentally distinct lineage from Tcf21+ stromal cells.

Next, we used targeted spatial transcriptomics [18] to multiplex the mRNA localisation of 18 cell type marker genes in six representative regions of adult murine lungs (n=2) along the proximal distal axis of the airway tree (figure 1d and supplementary figure S2a–c). Hhip+/Aspn+ peribronchiolar fibroblasts were enriched around airways, with some cells also in alveoli. Myh11+/Acta+ smooth muscle cells and Serpinf1+/Clec3b+ adventitial fibroblasts were enriched around airways and large vessels. Gucy1a3+/Postn+ pericytes localised preferentially to alveolar space and around larger vessels. The Tcf21hi/Npnt+ lipofibroblasts were localised preferentially to alveolar space (figure 1e and supplementary figure S2d and e). Consequently, the number of cells in close physical proximity (direct cell–cell contact) to Sftpc+ AT2 cells was highest for pericytes and lipofibroblasts with some Hhip+/Aspn+ cells also participating in the AT2 cell niche (figure 1f and supplementary figure S2f). Our data highlight the complexity of the AT2 cell niche, which we here demonstrate to contain at least three distinct stromal cell types.

An activated fibroblast state characterised by high Sfrp1 and Col28a1 expression

To follow the fate of fibroblasts during injury and repair we sorted cells from Tcf21m-Crem-R26R-tdTomato mice 14 days after bleomycin-induced lung injury (supplementary figures S1 and S3). Three injury-induced clusters became apparent (figure 2a and b), which were a mixture of Tcf21-lineage-positive and -negative cells. Interestingly, Tcf21-lineage-negative Hhip+ peribronchiolar fibroblasts expanded (supplementary figure S1j, p and q), and some markers (e.g. Lgr5 and Lgr6) were expressed in a myofibroblast subset (supplementary figure S1l), indicating that both Tcf21-lineage negative and positive cells converged into myofibroblasts.

FIGURE 2

Three distinct activated mesenchymal cell types emerge after bleomycin (bleo)-induced mouse lung injury. a) Uniform manifold approximation and projection (UMAP) depicts mesenchymal cell types (n=12 254 cells) identified at day 14 (d14) after bleomycin injury. b) The marker genes of activated mesenchymal subtypes are enriched for characteristic gene ontology (GO):biological process (BP) terms. False discovery rate <10%. Top 500 marker genes considered. c) The dotplot depicts a gradient of marker gene expression across the activated cell types. d) The mesenchymal cell types are characterised with distinct gene expression profiles, as depicted in the matrixplot. e, f) Immunofluorescence analysis of lung tissue sections from bleomycin-treated mice demonstrates colocalisation of SFRP1 (red) and COL28A1 (green) 14 days after injury. Arrowheads in the magnified insets indicate SFRP1/COL28A1 double-positive cells. SFRP1/COL28A1 double-positive cells did not colocalise with α-smooth muscle actin (ACTA)2 (cyan). No expression of SFRP1 was detected in PBS controls. Scale bars=100 µm (PBS), 100 µm (Bleo d14) and 50 µm (inset). f) At d14 after injury, ACTA2 (green) and CTHRC1 (red) double-positive cells were found (indicated by arrowheads in the magnified insets for ROI1 and ROI2 in the very right panel). Nuclei in blue colour (4′,6-diamidino-2-phenylindole DAPI). Scale bars=100 µm and 20 µm (inset). g, j) Distinct markers from the bleomycin mouse model were analysed in a recently published integrated atlas of human lung fibrosis at single-cell resolution [20]. g) SFRP1 showed specific expression on human COL1A2+ positive mesenchymal cells. h) Marker gene expression across all mesenchymal cell types identified in the human lung. i, j) Differential gene expression analysis between interstitial lung disease (ILD) patient samples and healthy control samples: i) reveals significant regulation of SFRP1 in adventitial fibroblasts (AdvF) and inflammatory fibroblasts cluster 2; j) this regulation is consistent across all three cohorts in the atlas. k) Immunofluorescence analysis of micro-CT staged idiopathic pulmonary fibrosis (IPF)-patient tissues demonstrating an increase of SFRP1 (red) in regions of mild fibrosis (IPF stage I) compared to healthy controls. In regions of severe fibrosis (IPF stage III) and tissue remodelling, fibroblastic foci (white dashed lines, FF) stained positive for α-smooth muscle actin (ACTA2) (green), but mostly not for SFRP1. Adjacent to fibroblastic foci (A-FF) SFRP1 expressing fibroblasts were still present. Scale bars=200 µm. Representative images are shown. l) Quantification of SFRP1's mean fluorescence intensity (MFI) of various regions-of-interest (ROIs) from three different patients (n=3). Statistics: one-way ANOVA. m) Quantification of SFRP1 MFI in late-stage IPF comparing SFRP1 expression within fibroblastic foci (FF) and adjacent to fibroblastic foci (A-FF) from three different patients (n=3). Statistics: unpaired t-test. ****: p<0.0001. SMC: smooth muscle cell; ECM: extracellular matrix.

We identified two Cthrc1+/Acta2+ myofibroblast types, with one subcluster showing enhanced Spp1 expression (figure 2c and d). One cluster coexpressed myofibroblast genes with different fibroblast markers, especially from lipofibroblasts; hence, we named these “transitional fibroblasts” (figure 2b–d). Several genes, including the secreted frizzled-related protein 1 (Sfrp1), showed highest expression in this putative intermediate cell population (figure 2c and d; supplementary table S2). Gene–gene correlation of Sfrp1 with other genes detected a transitional fibroblast core gene set, including Col28a1 (supplementary figure S4d). Tissue proteomics revealed transient induction of COL28A1 and SFRP1 proteins after bleomycin injury (supplementary figure S4c) [19]. Co-staining SFRP1 and COL28A1 confirmed coexpression in ACTA2-low or -negative transitional cells (figure 2e), while CTHRC1+ cells expressed higher levels of ACTA2 (figure 2f).

To address clinical relevance, we compared Col1a2+ cells after bleomycin injury with COL1A2+ cells in ILD patients [20]. SFRP1 was specific to COL1A2+ mesenchymal cells and was weakly expressed in some CLDN5+ endothelial cells (figure 2g). Expression was restricted to adventitial fibroblasts and disease enriched “inflammatory” fibroblast subsets (figure 2g and h). Increased expression of SFRP1 in disease-induced fibroblast states was consistent in all three study cohorts (figure 2i and j). Importantly, in IPF tissues, more SFRP1+ fibroblasts were present in mildly affected (early-stage) regions, while end-stage regions had more ACTA2+ myofibroblasts localised to fibroblast foci (figure 2k–m). These data suggests a similar trajectory of fibroblast states in human lung fibrogenesis as seen in the bleomycin model.

Sfrp1+ transitional fibroblasts precede the appearance of Cthrc1+ myofibroblasts

To follow transcriptional dynamics of Col1a2+ stromal cells during inflammatory, fibrogenic and resolution phases of lung regeneration, we collected Epcam−/Pecam1−/Lyve1−/CD45− stromal cells at 18 time points after bleomycin injury (figure 3a–c). This longitudinal scRNAseq dataset of Col1a2+ cells (69 185 cells) (figure 3b) featured three distinct injury-induced cell states (supplementary figure S4a) similar to the Tcf21-lineage tracing dataset (figure 3c and supplementary table S3). Sfrp1+ transitional fibroblasts peaked at day 3, preceding Spp1+ and Cthrc1+ myofibroblasts at day 9 to day 21 (figure 3d). Proliferation analysis showed a transient increase in proliferation rates, mainly in Spp1+ and Cthrc1+ myofibroblasts, returning to baseline from day 28 onwards (supplementary figure S4b). This suggests the early post-injury increase of Sfrp1+ cells was not from expanding pre-existing cells, but from differentiating baseline fibroblast states. Using immunofluorescence (figure 3e–l), we found that in healthy lungs, signals for SFRP1, CTHRC1 and SPP1 were mostly absent, except for nonspecific signals in the airway epithelium's luminal areas (figure 3e). Filamentous COL28A1 localised around bronchovascular cuffs, whereas ACTA2 was primarily observed in smooth muscle cells near airways and blood vessels (figure 3e). Localisation of adventitial and peribronchial fibroblasts next to bronchovascular cuffs has been described [1] and concurs with our single-molecule fluorescence in situ hybridisation (smFISH) analysis (figure 1).

FIGURE 3

Sfrp1+ transitional fibroblasts emerge in adventitial and alveolar space upon injury preceding the emergence of Cthrc1+ myofibroblasts. a) A high-resolution longitudinal dataset was generated by subjecting magnetic-activated cell-sorting (MACS)-sorted cells from the mesenchymal compartment to single-cell RNA-sequencing (scRNAseq) at the 18 indicated time points. Uniform manifold approximation and projection (UMAP) embedding displays cells coloured by b) time point and c) cell type identity. d) The distribution of cell type frequencies across time points. e) Immunofluorescence analysis of lung tissue sections from PBS control mice. In the left panel, cell nuclei were immunostained with 4′,6-diamidino-2-phenylindole DAPI (blue) and smooth muscle cells (SMCs) with α-smooth muscle actin (ACTA)2 (cyan). SFRP1 (red) staining was absent apart from nonspecific signals in bronchial epithelial cells. COL28A1 (green) was mostly found as a filamentous staining in “cuffs” surrounding blood vessels (bv). In the right panel, cell nuclei were stained with DAPI (blue) and SMCs with ACTA2 (cyan). CTHRC1 (red) staining was absent apart from nonspecific signals in bronchial epithelial cells. SPP1 is depicted in green. Scale bars=100 µm. f) Immunofluorescence analysis of lung tissue sections from bleomycin-treated mice at day 3 after injury (Bleo d3) demonstrating appearance of SFRP1 (red) positive cells surrounding blood vessels and airways (aw) (arrowheads). Scale bar=100 µm. g) Arrowheads in the magnified insets taken from region of interest (ROI)1 and ROI2 of figure 4f indicate SFRP1/COL28A1 double-positive cells surrounding blood vessels in ROI1 reminiscent of adventitial fibroblasts, as well as surrounding airways (aw) in ROI2 reminiscent of peribronchiolar fibroblasts. Scale bars=20 µm. h) Yellow dashed line indicates a fibrotic region of lung tissue after 14 days of bleomycin treatment. ACTA2 staining in cyan exhibits the appearance of myofibroblasts, and concomitant expression of SFRP1 (red) and COL28A1 (green) in this fibrotic region. Cell nuclei are stained with DAPI in blue. Scale bar=100 µm. A magnified view of fibrotic ROI1 is displayed in (i), in the left panel demonstrating mutually exclusive appearance of SFRP1+ (red and dashed white lines) and ACTA2+ (green and dashed yellow lines) cells. Here, for better detection of colocalisation signals between ACTA2 and SFRP1, we switched colours of ACTA2 from cyan as depicted in (h) to green. The right panel denotes SFRP1 (red) and COL28A1 (green) double-positive cells encircled with white dashed lines. Cell nuclei stained with DAPI in blue. Scale bars=20 µm. j) Iterative immunofluorescence staining of parenchymal lung tissue (4i-FFPE) at day 3 after injury indicating COL28A1 (green)/SFRP1 (red) double-positive cells in the left image. The same cells, as indicated by white dashed lines, were found to be CTHRC1- (red) and ACTA2- (cyan) negative. Scale bar=50 µm. A larger overview image can be found in supplementary figure S7a. k) Stacked bar-plot denoting relative frequencies of ACTA2-, CTHRC1-, SPP1-, COL28A1- and SFRP1-positive tissue segments at day 3 and day 14 after bleomycin treatment and compared to PBS controls from software-based segmented images as exemplified in supplementary figure S7c. Three different ROIs (each 1.1 mm2 in size) from two different mice for each condition were analysed. l) Quantification of double-positive SFRP1+/COL28A1+ and CTHRC1+/ACTA2+ tissue segments were analysed by fluorescent-signal colocalisation (supplementary figure S7d), its segmentation and software-based quantification. The number of colocalising tissue segments relative to the total cell count (20 691 cells) is shown. 10–20 different 0.1–0.25-mm2 ROIs (peribronchial/perivascular/parenchymal on day 3 compared to fibrotic on day 14) from two different mice per condition were analysed. Statistics: two-way ANOVA with Tukey's multiple comparison: **: p<0.01.

3 days post-injury, SFRP1/COL28A1 double-positive cells emerged around airways, blood vessels and alveoli, with ACTA2 expressed in bronchovascular smooth muscle cells (figure 3f, g and j). In contrast, at day 14 post-injury, ACTA2-myofibroblasts and SFRP1+/COL28A1+ cells coexisted in fibrotic dense areas (figure 3h and i). Multiplexed immunofluorescence-analysis confirmed the increase of SFRP1+/COL28A1+ transitional cells at day 3, and until day 14 post-injury expanding SPP1+ and CTHRC1+ myofibroblasts (figure 3k and l; total analysis of 20 691 cells, supplementary figure S7c), thus verifying our data in figure 3d on protein level.

Convergence of multiple mesenchymal cell types towards myofibroblast identity

To infer lineage relationships between the different cell types and activation states we used CellRank (figure 4) [21]. We analysed three phases as “early injury” (figure 4a–c), “fibrogenesis” (figure 4d–f) and “resolution” (figure 4g–i). CellRank computed fate probabilities for each cell, showing their specific differentiation potential towards a “terminal state” end-point in the time-course data. Heatmaps demonstrate fate probability (figure 4a, d and g), and violin plots visualise the probability distribution among the cells for selected fates in each cell type (figure 4b, e and h). Key lineage driver genes correlating with different fate trajectories are depicted in scatter plots (figure 4c, f and i). Early injury phase analysis (days 2–5) showed that adventitial fibroblasts and lipofibroblasts had a high fate probability of transitioning to the Sfrp1+ state (figure 4a and b). Probabilities towards myofibroblasts were low, represented only by few Cthrc1+ myofibroblasts (figure 4a and b). Top driver genes towards transitional cells included various collagens, chemokines and notably Sfrp1, while the myofibroblast lineage showed marker genes including Tnc, Spp1 and Thbs1 (figure 4c). In the fibrogenesis phase (day 6 to day 21), fate probabilities predicted Sfrp1+ cells transitioning to myofibroblasts (figure 4d,e). Top driver genes towards myofibroblasts represented classical myofibroblast-associated genes including Lgals1, Sparc and Spp1 as well as Cthrc1 (figure 4f). In the resolution phase, differentiation probabilities suggested a reversion of Cthrc1+ myofibroblasts towards lipofibroblast, peribronchiolar fibroblast and pericyte states (figure 4g and i). This prediction correlates with previous observations of a two-way conversion between lipogenic and myogenic fibroblastic phenotypes [6].

FIGURE 4

Lineage convergence and two-way conversion of fibroblasts in lung regeneration. a, d, g) CellRank's coarse-grained and directed transition matrix was calculated for the given three phases of the bleomycin lung injury time course, of early phase (day 2 to day 5), fibrogenesis phase (day 6 to day 21) and resolution phase (day 28 to day 54). Terminal macrostates were set manually and annotated according to their overlap with the underlying gene expression clusters. The heatmap represents the mean absorption probabilities of every cell in the given cell types to the terminal fates. b, e, h) Violin plots indicating the absorption probability towards the specified fate for every cell in all the cell types. c, f, i) Scatter plots show top driver genes, predicted to facilitate transition to the given terminal fate, visualised as gene-correlations between two terminal fate lineages. The top genes of a third lineage are visible by the third colour of the gene names. Dots are coloured according to their mean gene expression. SMC: smooth muscle cell.

Next, we used a spline regression model revealing genes with differential expression in at least one cell type over time (supplementary table S4). We found 25 different collagens with dynamic expression patterns after injury or between various stromal cells (figure 5a). Additionally, 90 secreted matrisome [22] genes with significant expression changes along the injury time-course were identified (figure 5b). All cell types displayed transiently increased expression of collagen type-I (Col1a2) and the ECM protein Sparc, a marker for myofibroblasts (figure 5c) [23]. We also observed common early fibroblast activation events, like the early post-injury upregulation of S100a6 (figure 5c). Peribronchiolar fibroblasts specifically featured important secreted morphogens like Hhip, Wif1, Fgf18 and Wnt5a (figure 5d). This matters as fibroblast-derived Wnt-signaling, especially Wnt5a, defines a specific AT2 cell niche in normal homeostasis [24], is part of a distinct mesenchymal niche in human distal airways, and is secreted by LGR5+ fibroblasts [17]. Adventitial fibroblasts specifically produced the Wnt modulator Sfrp4, besides injury-induced chemokines like Cxcl13 for B-cell recruitment, and Cxcl12 and Ccl8 to attract T-cells, monocytes and neutrophils. Interestingly, the chemokine eotaxin (Ccl11) specifically marked adventitial fibroblasts, unaffected by injury (figure 5e). Lipofibroblasts specifically expressed the morphogens Wnt2 and Bmp3, plus stem cell factor (SCF) (Kitl) (figure 5f). The SCF-c-Kit pathway is activated in bleomycin-injured lungs, with potential profibrotic effects via recruitment of Kit+ immune cells to the lung [25].

FIGURE 5

Highly specific expression changes in distinct fibroblast lineages. a, b) Heatmaps display genes that showed differential expression along the time-course in at least one cell type using a spline regression model for a) 25 different collagens and b) 90 secreted matrisome [22] genes. c–g) Line plots show average expression of indicated genes for each cell type along the bleomycin time course, for c) general mesenchymal marker genes, d) genes regulated over time in peribronchial fibroblasts, e) adventitial fibroblast, f) lipofibroblasts and g) Cthrc1+ myofibroblasts. SMC: smooth muscle cell.

Conclusively, our data suggest that early after injury multiple mesenchymal cell types converge towards Sfrp1+ transitional cells, which ultimately give rise to Cthrc1+ myofibroblasts, with a potential to revert towards lipofibroblasts, peribronchiolar fibroblasts and pericytes during the resolution phase.

TGFβ mediates differentiation of Sfrp1+ transitional fibroblasts into myofibroblasts

To infer potential regulators of the fibroblast cell state transitions after injury we used NicheNet [26], which predicted the ligands with highest probability to induce expression of the top driver genes from the CellRank outputs (figure 4) for Sfrp1+ transitional cells and Cthrc1+ myofibroblasts (figure 6a and b). The top ligand upstream of Sfrp1+ transitional cells was the Notch ligand Jag1, predominantly expressed in secretory airway epithelial cells, and to a lesser degree in alveolar epithelial and vascular endothelial cells (figure 6a). Consistent with our finding here, Notch deficiency in mesenchymal cells reportedly reduced fibrotic remodelling and myofibroblast differentiation in the bleomycin model [27]. Notably, upstream of Sfrp1+ transitional cell state driver genes the top ligands excluded Tgfβ, which was the top ranked ligand for Cthrc1+ myofibroblast driver genes. Tgfβ1 was the top ligand, primarily expressed by myeloid lineage immune cells, along with Tgfβ3 expressed by activated AT2 cells and lymphatic endothelial cells (figure 6b).