Single-nucleus chromatin accessibility identifies a critical role for TWIST1 in idiopathic pulmonary fibrosis myofibroblast activity

Abstract

Background: In idiopathic pulmonary fibrosis (IPF), myofibroblasts are key effectors of fibrosis and architectural distortion by excessive deposition of extracellular matrix and their acquired contractile capacity. Single-cell RNA-sequencing (scRNA-seq) has precisely defined the IPF myofibroblast transcriptome, but identifying critical transcription factor activity by this approach is imprecise.

Methods: We performed single-nucleus assay for transposase-accessible chromatin sequencing on explanted lungs from patients with IPF (n=3) and donor controls (n=2) and integrated this with a larger scRNA-seq dataset (10 IPF, eight controls) to identify differentially accessible chromatin regions and enriched transcription factor motifs within lung cell populations. We performed RNA-sequencing on pulmonary fibroblasts of bleomycin-injured Twist1-overexpressing COL1A2 Cre-ER mice to examine alterations in fibrosis-relevant pathways following Twist1 overexpression in collagen-producing cells.

Results: TWIST1, and other E-box transcription factor motifs, were significantly enriched in open chromatin of IPF myofibroblasts compared to both IPF nonmyogenic (log₂ fold change (FC) 8.909, adjusted p-value 1.82×10^-35) and control fibroblasts (log₂FC 8.975, adjusted p-value 3.72×10^-28). TWIST1 expression was selectively upregulated in IPF myofibroblasts (log₂FC 3.136, adjusted p-value 1.41×10^{- 24}), with two regions of TWIST1 having significantly increased accessibility in IPF myofibroblasts. Overexpression of Twist1 in COL1A2-expressing fibroblasts of bleomycin-injured mice resulted in increased collagen synthesis and upregulation of genes with enriched chromatin accessibility in IPF myofibroblasts.

Conclusions: Our studies utilising human multiomic single-cell analyses combined with in vivo murine disease models confirm a critical regulatory function for TWIST1 in IPF myofibroblast activity in the fibrotic lung. Understanding the global process of opening TWIST1 and other E-box transcription factor motifs that govern myofibroblast differentiation may identify new therapeutic interventions for fibrotic pulmonary diseases.

FIGURE 1

Single-cell RNA-sequencing (scRNA-seq) and single-nucleus assay for transposase-accessible chromatin sequencing (snATAC-seq) profiling of the human idiopathic pulmonary fibrosis (IPF) and healthy lung. a) Uniform manifold approximation and projection (UMAP) plot snATAC-seq dataset (IPF n=3, control n=2) identified by cluster number and cell type. b) Merged coverage plots demonstrating pseudo-bulk chromatin accessibility (fragment coverage by frequency of Tn5 insertion) around marker gene promoters. Y-axis cluster numbers correspond to cell clusters in figure 1a. Range of normalised accessibility for fragment coverage of each gene listed on x-axis. c) UMAP plot of scRNA-seq dataset (IPF n=10, control n=8) identified by cell type. d) Dot plot of scRNA-seq dataset showing gene expression of selected cell-type specific marker genes. The diameter of the dot corresponds to the proportion of the cells expressing the gene, and the colour density of the dot corresponds to the average expression level relative to all cell types. e) Heatmap of average number of Tn5 cut sites within the differentially accessible regions (DARs) (each row is a unique DAR) for each cell type. The colour scale represents a z-score of the number of Tn5 sites within each DAR. f ) Genomic annotation of differentially accessible region locations. Mac: macrophage; AT: alveolar type; SM: smooth muscle; NK: natural killer; pDC: plasmacytoid dendritic cell; UTR: untranslated region.

FIGURE 2

Fibroblast subpopulations and transcription factor motif activity in idiopathic pulmonary fibrosis (IPF) myofibroblasts. a) Uniform manifold approximation and projection (UMAP) plot of single-cell RNA-sequencing (scRNA-seq) fibroblasts, smooth muscle cells and pericyte clusters by cell identity. b) UMAP plot of scRNAseq fibroblasts, smooth muscle cells and pericyte clusters from (a) with cells depicted by origination from IPF versus control samples. c) UMAP plot of single-nucleus assay for transposase-accessible chromatin sequencing fibroblast, smooth muscle cells and pericyte clusters by cell identity. d) Heatmap of average number of Tn5 cut sites within the differentially accessible regions when comparing IPF myofibroblasts to IPF nonmyogenic fibroblasts, depicted by individual sample. Each row is a unique differentially accessible region (DAR). e) Ingenuity pathway analysis pathways significantly enriched for genes annotated to the upregulated DARs in IPF myofibroblasts versus IPF nonmyogenic fibroblasts. f ) Ingenuity pathway analysis pathways significantly enriched for genes annotated to the upregulated DAR in IPF fibroblasts versus control fibroblasts. g) Transcription factors with the most significantly enriched motif activity when comparing IPF myofibroblasts to the IPF nonmyogenic fibroblasts. h) Transcription factors with the most significantly enriched motif activity when comparing all IPF fibroblasts to all control fibroblasts. EMT: epithelial–mesenchymal transition; BMP: bone morphogenetic protein; FC: fold change; adj: adjusted.

FIGURE 3

TWIST1 binds in areas of in vivo accessible chromatin in idiopathic pulmonary fibrosis (IPF) fibroblasts. a) Coverage plot demonstrating Tn5 insertion frequency (pseudobulk single-nucleus assay for transposase-accessible chromatin sequencing (snATAC-seq) tracks) in IPF and control fibroblasts in the region of KANK3 with ATAC-seq peaks depicted by red bars. TWIST1 Ab chromatin immunoprecipitation sequencing (ChIP-seq) peaks from IPF explant fibroblasts (grey tracks labelled IPF 1–4) in the same region are depicted below. Grey shading indicates the region of TWIST1 binding by ChIP-seq experiments. Scales for ATAC-seq tracks and ChIP-seq tracks are independent. b) Coverage plot demonstrating Tn5 insertion frequency (pseudobulk ATAC-seq tracks) in IPF and control fibroblasts in the region of DYNC1H1 with ATAC-seq peaks depicted by red bars. TWIST1 Ab ChIP-seq peaks from IPF explant fibroblasts (grey tracks labelled IPF 1–4) in the same region are depicted below. Light grey box indicates the region of TWIST1 binding by ChIP-seq experiments. Scales for ATAC-seq tracks and ChIP-seq tracks are independent.

FIGURE 4

TWIST1 expression and motif activity in idiopathic pulmonary fibrosis (IPF) and control fibroblasts. a) Uniform manifold approximation and projection (UMAP) plot of TWIST1 motif activity in the mesenchymal single-nucleus assay for transposase-accessible chromatin sequencing (snATAC-seq) clustering, depicted by scaled expression with red representing the highest motif activity. b) Violin plot of TWIST1 motif activity comparing IPF versus control fibroblasts only. c) UMAP plot of TWIST1 gene expression in the mesenchymal single-cell RNA-sequencing (scRNA-seq) clustering demonstrating high TWIST1 expression in the myofibroblasts only. d) Coverage plot demonstrating Tn5 insertion frequency by snATAC-seq fibroblast cluster in the TWIST1 gene region. Red boxes represent regions of statistically significant differential accessibility in IPF myofibroblasts versus nonmyogenic fibroblasts. e) Ingenuity pathway analysis pathways significantly enriched for genes annotated to the upregulated differentially accessible regions (DAR) in TWIST1-open chromatin myofibroblasts versus TWIST1-closed chromatin myofibroblasts. f ) Dot plot of E-box transcription factor gene expression in IPF and control mesenchymal populations by scRNA-seq with deeper colour indicating higher level of gene expression and circle size indicating the percentage of cells in the population expressing the gene. NFAT: nuclear factor of activated T-cells; EMT: epithelial–mesenchymal transition; FGF: fibroblast growth factor.

FIGURE 5

Overexpression of Twist1 in Col1a2+ cells leads to increased collagen I levels in vitro and in vivo. a) A triple transgenic animal was bred where Cre recombinase is under control of the col1a2 enhancer element (col1a2-Cre-ER(T)). In the presence of tamoxifen (TAM), the STOP signal is excised leading to expression of the reverse tetracycline transactivator (rtTA). In the presence of doxycycline (DOX) and the rtTA, the tetO7 operator is activated leading to expression of Twist1 and luciferase. b) Lung fibroblasts from Twist1-WT (wild-type) and Twist1-Luc (Twist1 overexpressors) were incubated in the presence of TAM and DOX with and without transforming growth factor (TGF)-β (2 ng·mL−1 ). Cells were lysed and subjected to immunoblotting for collagen I, α-smooth muscle actin (SMA), Twist1 and the loading control, cyclophilin A. In the presence of TAM/DOX, increased c) collagen, d) α-SMA and e) Twist1 in Twist1-LUC fibroblasts. This was amplified in the presence of TGF-β (n=3). Data were analysed by robust nonparametric two-way ANOVA. p-values for the effects of Twist1 and TGF-β and the interaction are presented in the panels. f ) Twist1-WT and Twist1-Luc mice were injured with bleomycin. Animals were sacrificed at 14 days. Lungs were excised, and the ratio of the lung mass before and after freeze-drying was determined. g) Determination of acid-soluble collagen content showed a significant increase in bleomycin-induced collagen in Twist1-Luc mice compared to Twist1-WT mice (by robust nonparametric two-way ANOVA; n=6–8). h) A comparable degree of histological injury was observed in Twist1-WT and Twist1-Luc mice. Haematoxylin and eosin (H&E) and trichrome images are presented. Incidentally noted multinucleated giant cells are identified and magnified by the black arrows. Scale bar=200 μm, inset ×100 magnification. NS: nonsignificant.

FIGURE 6

Twist1 overexpression in mouse lung fibroblasts is associated with dysregulation of several pulmonary fibrosis genes and pathways. Bulk RNA-sequencing (seq) was performed on fibroblasts isolated from lungs of wild-type (WT) and Twist1-LUC mice (n=3). Estimation of differential gene expression using CLC Genomics Workbench was performed comparing fibroblasts from knock-in with normal lungs. a) Hierarchical clustering heatmap of significant differentially expressed genes was generated using CLC Genomics Workbench using minimum absolute fold change of 3.0 and false discovery rate (FDR) p-value threshold of 0.05. Immunoblotting is shown for the individual lines subjected to RNA-seq. Densitometry normalised to β-actin is shown beneath. b) Volcano plot shows comparative analysis of differentially expressed genes between the WT and Twist1-LUC. c) List of dysregulated genes by FDR that are upregulated ranked by −log10 (FDR p-value) with p<0.05 cut-off. d) List of top downregulated genes ranked by −log10 (FDR p-value) with p<0.05 cut-off downregulated in Twist1-LUC fibroblasts compared to Twist1-WT fibroblasts. e) Ingenuity pathway analysis (IPA) of dysregulated canonical pathways between WT and Twist1-LUC (n=3) by z-score and FDR.