Distinct roles of KLF4 in mesenchymal cell subtypes during lung fibrogenesis

During lung fibrosis, the epithelium induces signaling to underlying mesenchyme to generate excess myofibroblasts and extracellular matrix; herein, we focus on signaling in the mesenchyme. Our studies indicate that platelet-derived growth factor receptor (PDGFR)-β+ cells are the predominant source of myofibroblasts and Kruppel-like factor (KLF) 4 is upregulated in PDGFR-β+ cells, inducing TGFβ pathway signaling and fibrosis. In fibrotic lung patches, KLF4 is down-regulated, suggesting KLF4 levels decrease as PDGFR-β+ cells transition into myofibroblasts. In contrast to PDGFR-β+ cells, KLF4 reduction in α-smooth muscle actin (SMA)+ cells non-cell autonomously exacerbates lung fibrosis by inducing macrophage accumulation and pro-fibrotic effects of PDGFR-β+ cells via a Forkhead box M1 to C-C chemokine ligand 2—receptor 2 pathway. Taken together, in the context of lung fibrosis, our results indicate that KLF4 plays opposing roles in PDGFR-β+ cells and SMA+ cells and highlight the importance of further studies of interactions between distinct mesenchymal cell types.

Lung cells were subjected to scRNA library construction (10X Genomics), sequencing and then annotating and clustering. Uniform manifold approximation and projection (UMAP) of 720 non-immune lung cells is shown. a, Cells are colored by their categorized cell-type identity. b-g, The same UMAP, colored by normalized expression levels of Pdgfrb, ZsGreen1 (Zs) reporter, Gli1, Pdgfra, Axin2, Acta2, Tagln, Des and Klf4 as indicated. Unique molecular identifiers are scaled to ten thousand transcripts per cell and then natural log transformed with a pseudocount of one.  Figure 6. Cell marking with Acta2-CreERT2. Acta2-CreERT2, ROSA26R(mTmG/+) mice were induced with tamoxifen, rested, injected with a single orotracheal dose of PBS (a, b) or bleomycin (c) and euthanized fourteen days later. Lungs were stained for SMA, GFP (fate marker) and nuclei (DAPI). In a, boxed regions are shown as close-ups below. In c, fibrotic and non-fibrotic regions are indicated by yellow and orange demarcations, respectively. n=3 mice. Scale bars, 50 μm. Pdgfrb-CreERT2 mice were or were not induced with tamoxifen as indicated and then rested and exposed to normoxia or hypoxia (FiO2 10%) for 21 days. Lung vibratome sections were analyzed. In a, mice also carrying ROSA26R(mTmG/+) were induced with tamoxifen, and sections were stained for SMA, GFP (fate marker) and MECA-32 (ECs). Per treatment group, n=2 mice, 8 sections per mouse were analyzed. In b, sections from mice also carrying Klf4(flox/flox) were stained for SMA, MECA-32 and nuclei (DAPI). For each experimental group, n=3 mice and at least 10 sections per mouse were stained. Scale bars, 25 μm.  Figure 8. Following bleomycin treatment, patches of clones with more cells cover a larger area than clones with fewer cells. Pdgfrb-CreERT2, ROSA26R(Rb/+) were labelled with tamoxifen, rested, subjected to a single orotracheal bleomycin dose and euthanized 14 days later. Lung cryosections were stained with DAPI and directly imaged for Rb colors (mCherry, mOrange, Cerulean). Among clones consisting of 5-7, 10-12 or 20-22 cells, clones were chosen randomly and the clonal patch area was measured. n=3 mice, 5 patches per each clone cell size per mouse. One-way ANOVA with Tukey's multiple comparisons test was used. Patch size 5-7 vs. 10-12 (p=0.011), 10-12 vs. 20-22 (p=0.0002), *** vs. 5-7 (p<0.0001). Data are averages +/-SD. Source data are provided as a Source Data file.   Pdgfrb-CreERT2, ROSA26R(YFP/+) mice also carrying Klf4(flox/flox) or wild type for Klf4 were induced with tamoxifen, rested, treated with a single orotracheal dose of bleomycin and then euthanized five days later. EdU was administrated 12 hours before euthanasia. In b, lung cryosections were stained for YFP (fate marker), EdU and nuclei (DAPI) with arrowheads indicating proliferating YFP+ cells. Scale bar, 25 μm. In c, percent of YFP+ cells that are EdU+ is shown. For each genotype, n=3 mice, 3 sections per mouse and an average of 48 YFP+ cells were analyzed per section. Two-tailed Student's t-test; ns, not significant. Data are averages +/-SD. Source data are provided as a Source Data file.

Supplementary Figure 12. Isolation of Zs+ cells from the lungs of Pdgfrb-CreERT2, ROSA26R(Zs/+) mice.
Lungs were harvested from wild type and tamoxifen-injected Pdgfrb-CreERT2, ROSA26R(Zs/+) mice, and a single cell suspension was generated. Cells were stained with DAPI, and Zs+DAPI-cells were isolated by FACS. a, FACS plots of wild type (left panel) and Pdgfrb-CreERT2, ROSA26R(Zs/+) (right panel) mice. Lung cells from wild type mice were used to control for auto-fluorescence and set gating for Zs+ cell isolation. Box within quadrant #1 (Q1) indicates the isolated Zs+DAPI-cells. n=5. b, Isolated Zs+ cells were cultured for two weeks and then stained for nuclei (DAPI) and directly imaged for Zs. n=5. Scale bar, 25 μm.  Figure 13. With bleomycin administration, KLF4 is dynamically upregulated in SMA+ lung cells but Klf4 deletion in SMA+ cells does not alter proliferation of this lineage. Acta2-CreERT2, ROSA26R(Zs/+) mice were induced with tamoxifen and rested. a, Mice were or were not treated with a single orotracheal dose of bleomycin, and at 1, 3 and 5 days later, Zs+ cells were isolated by FACS and subjected to qRT-PCR for Klf4 and Gapdh. Transcript levels of Klf4 relative to that of Gapdh for each time point and normalized to control (no treatment) are shown. For each time point, n=3 mice in triplicate. One-way ANOVA with Tukey's multiple comparisons test was used. Control vs. day 3 (p=0.0018), day 3 vs. day 5 (p=0.0033). ns, not significant vs. control. b, c, Mice were also carrying Klf4(flox/flox) or were wild type for Klf4. EdU was administrated 12 hours before mice were euthanized on day 14 after bleomycin, and lung sections were stained for EdU and nuclei (DAPI) and directly imaged for Zs. In b, boxed areas are shown as close-ups on the right, and arrowheads indicate EdU+Zs+ cells. Scale bar, 25 μm. In c, the percentage of Zs+ cells that are EdU+ is shown for Klf4 wild type and floxed groups. For each genotype, n=3 mice, 3 sections per mouse and an average of 67 Zs+ cells per section were analyzed. Two-tailed Student's t-test was used. Data are averages +/-SD. Source data are provided as a Source Data file.

increases Pyruvate dehydrogenase kinase (Pdk) 1 transcript levels in SMCs but not PDGFR-β+ cells. Murine lung SMCs or Zs+ cells isolated from the lungs of tamoxifen-induced
Pdgfrb-CreERT2, ROSA26R(Zs/+) mice were subjected to Scr RNA or siKlf4 treatment. Transcript levels of Klf4 and Pdk1 were determined by bulk RNA-seq (a) or qRT-PCR (b). In b, transcript levels were relative to that of Gapdh.
In a, b, results of siKlf4 treatment were normalized to results of Scr treatment. n=3 per condition, each n in b was done in triplicate. FC, fold change. FDR, false discovery rate. For a, DESeq2 was used which employs two-sided Wald test to determine the p-value and Benjamini-Hochberg FDR for p-adjusted value. For b, two-tailed Student's t-test was used; *, p=0.011; ***, p=0.0004; ****, <0.0001; ns (not significant, p=0.068) vs. Scr. Data are averages +/-SD. Source data are provided as a Source Data file.

Supplementary Methods
Sample preparation for single-cell RNA sequencing

c) Sequencing libraries
The single cell 3' library comprised standard Illumina paired-end constructs which begin and end with P5 and P7. The single cell 3' 16 bp 10x Barcode and 10 bp UMI were encoded in Read 1, and Read 2 was used for sequencing the cDNA fragment. Sequencing the single cell 3' library produced standard Illumina Binary Base Call data which includes the paired-end Read 1 and Read 2 and the sample index in the i7 index read.
Processing scRNA-seq data scRNA-seq data was processed with Cell Ranger v3.1.0. 62.7% of reads were confidently aligned to a modified version of the mouse transcriptome mm10 that includes the sequence for the gene ZsGreen1. The top cell barcodes selected by Cell Ranger were then utilized for downstream analysis.

Cell barcode clustering and annotation
All analyses were performed in R version 3.6.1. Graph embedding and clustering were performed using the R package Seurat v3.1.0. UMI counts were scaled to 10,000 UMIs per cell, then natural log transformed with a pseudo-count of one (log((TPM/100) + 1)). in suspension between tissue digestion and droplet-encapsulation for scRNA-seq.
Following cell annotation, immune cells and multiplets were discarded; the remaining clusters of mesenchymal, endothelial and epithelial cells were then used to generate a Uniform Manifold Approximation and Projection (UMAP) embedding to visualize expression of specific transcripts. Cell hash-tagging was not included in the analysis due to poor correlation between tagging and ZsGreen1 expression.

UMAP embedding
The top 2,500 variable genes were selected with the Seurat FindVariableFeatures implementation and then scaled and used to generate principal components with the ScaleData and RunPCA implementations respectively. The final UMAP was generated by using the top 21 principal components with the Seurat RunUMAP implementation, with a neighborhood size of 10, a minimum distance parameter of 0.8 and 5,000 epochs with a learning rate of 0.5. Following first-strand synthesis with random primers, second strand synthesis and A-tailing were performed with dUTP to generate strand-specific sequencing libraries. Adapter ligation with 3' dTMP overhangs were ligated to library insert fragments. Library amplification of fragments carrying the appropriate adapter sequences at both ends was undertaken. Strands marked with dUTP were not amplified. Indexed libraries were quantified by qRT-PCR using a commercially available kit (Roche, KAPA Biosystems) and insert size distribution was determined with the Agilent Bioanalyzer. Samples with a yield of ≥ 0.5 ng/µl and a size distribution of 150-300 bp were used for sequencing.

c) Flow cell preparation, sequencing and data analysis
Samples at a concentration of 1.2 nM were loaded onto an Illumina NovaSeq6000 flow cell to yield 25 million passing filter clusters per sample. Samples were sequenced using 100 bp paired-end sequencing per Illumina protocols. Reads were trimmed to remove low quality base calls. HISAT2 was used to align trimmed reads to the reference genome mm10 with GENCODE annotation for mouse. Gene counts and transcript abundance were estimated using ballgown/stringTie and differential gene expression analysis was performed using R-based DESeq2. The log2 fold change was calculated by dividing knockdown values by Scr values and significance is determined from 3 independent RNA-seq experiments.

Bioinformatic analysis of Ccl2 binding TFs and cross reference to bulk RNA-seq data
The sequence spanning from 5 kb upstream of the mouse Ccl2 transcription start site to 100 bp downstream was retrieved from the UCSC Genome Browser (https://genome.ucsc.edu/ UCSC Genomics Institute, CA). This sequence was used for a TF analysis using TRANSFAC (geneXplain GmbH, Germany). TRANSFAC MATRIX TABLE, Release 2020.2 was used as matrix library, profile was set to vertebrate_non_redundant_minFP.prf, only high-quality matrices and cut-off set to minimize false positives. The resulting TFs were compared to those differentially expressed genes (DEGs) in the bulk RNA-seq data, for SMCs and PDGFR-b + cells with Klf4 knockdown, having log2 fold-change ≥ ±1.0 and adjusted p value < 0.05 (Benjamini-Hochberg false discovery rate) using Ingenuity Pathway Analysis (Version 52912811, Ingenuity Systems, QIAGEN). Overlapping TFs between TRANSFAC and the most upregulated genes in SMCs were assessed for subsequent analysis.