Mapping scRNAseq Atlas of Multiple Lung Cell Types in Idiopathic Pulmonary Fibrosis
The demographic and clinical characteristics of the scRNAseq and scATATCseq samples are listed in Table S2. To mitigate batch effects, an anchor-based strategy was employed and no significant batch effects were observed (Fig. S1A–D). After rigorous filtering and quality control, our scRNAseq atlas contained 409,303 cells, including 58,470 epithelial cells, 14,105 stromal cells, and 305,950 immune cells (Fig. 2A). The expression levels of marker genes, which were used for cell type annotation and identification of cell clusters, are depicted in the feature plots and dotplots, aligning with existing knowledge (Fig. S2A,C–E). In accordance with the prevailing understanding of the aetiology of IPF, our findings demonstrated fewer alveolar epithelial cells in the IPF group compared to the control group (Fig. 2C). Conversely, the IPF patients exhibit an increase in epithelial secretory cells and ciliated cells, which may contribute to enhanced elimination of environmental risk factors and the promotion of fibrosis development. Additionally, the IPF group displayed a higher abundance of stromal cells than the control group (Fig. S2B). Furthermore, a phenotypic transition from fibroblasts to myofibroblasts was also evident in the IPF patients (Fig. 2B).
SPP1_RecMacs Originate CD14+_Monocytes Rather Than Tissue-Resident Alveolar Macrophages
Fifteen distinct immune cell types were identified, as shown in Fig. 2A. Among these, we identified eight monocyte and macrophage clusters, including CD14+_Monocytes, CD16+_Monocytes, Mon_macs, CCL2_RecMacs, interstitial macrophages, SPP1_RecMacs, proliferating_Macs, and TRMs. Interestingly, a decreased proportion of CD14+_Monocytes was observed, while a trend was noted towards an increase in Mon-macs, CCL2_RecMacs, and SPP1_RecMacs (Fig. 3F). Compared to Mon-macs and CCL2_RecMacs, SPP1_RecMacs exhibited a higher level of CHIT1, SPP1, SDC2, and CHI3L1 (Fig. 3A, Fig. S2E). CHIT1, a well-known chitinase, was implicated in various inflammatory lung and fibrotic diseases [29–31]. CHI3L1 was also found to contribute to the development of lung fibrosis [33]. Notably, SDC2, encoding a glycosylated integral membrane protein, was discovered to alleviate fibrosis by decreasing TGFß1 receptors in epithelial cells [32].
Considering the crucial role of Mo_AMs in the development of IPF and the unclear origin of SPP1_RecMacs, we conducted Slingshot and Monocle 2/3 analysis to capture a linear pseudotime process starting from CD14+_Monocytes and progressing towards Mon-macs, CCL2_RecMacs, and SPP1_RecMacs as the trajectory endpoint (Fig. 3B–C,E). Among the top 20 variable genes associated with this trajectory, APOE exhibited an increase, which was in line with previous research findings [33] (Fig. 3D). Furthermore, our analysis revealed that the upregulated genes within the trajectory exhibited enrichment in the “cholesterol metabolism” KEGG term (Fig. 3D). The TF analysis revealed that both PU.1 and PPARγ tend to be expressed higher in SPP1_RecMacs than in TRMs (Fig. S3A,B). Moreover, SPP1_RecMacs also displayed a significantly lower phagocytosis score and higher tissue repair score than TRMs (Fig. S3C–D, Table S2). lastly, SPP1_RecMacs exhibited greater metabolic activities than TRMs in IPF patients, including lipid and phospholipid metabolism, glucose metabolism, and energy metabolism (Fig. S3E).
Identify the Mo_AMs Subtype Associated With FVC%pred
In this study, we analysed 119 IPF bulk samples with FVC%pred information and used the Scissor package to identify fibrosis-related Mo_AMs among a total of 57,005 cells. Of these cells, 5,227 Scissor − cells were negative for FVC%pred, and 4,664 Scissor + cells were positive for FVC%pred (Fig. 4A). As previously indicated, Scissor − cells predominantly occupy the latter portion of the trajectory, specifically CCL2_RecMacs and SPP1_RecMacs, suggesting the crucial role of Mo_AMs in promoting fibrosis (Fig. 4B). This result is further enhanced by the statistically significant p-value of the reliability significance test (p < 0 .001). To shed light on the transcriptional differences of Scissor + and Scissor − cells, we compared their DEGs. Our analysis revealed an upregulation of macrophage markers (MACRO, FN1, EMP1, and CHIT1) in Scissor − cells, while macrophage markers (AREG, CD163, CD68, and LYZ) were found to be greater in Scissor + cells (Fig. 4C–F).
Fibroblasts Induce Chemotaxis of Mo_AMs by CXCL12/CXCR4 Axis in Idiopathic Pulmonary Fibrosis
Intercellular communication among epithelial cells, myofibroblasts, fibroblasts, and alveolar macrophages was observed in patients with IPF. The use of CellChat revealed the CXCL12/CXCR4 signalling pathway between Mo_AMs and fibroblasts as well as myofibroblasts (Fig. 4G,I). The reliability of the communication was confirmed by the expression of CXCL12 and CXCR4 (Fig. 4H). To further validate these findings, coculture experiments of the fibrosis models and Mo_AMs were performed. To establish the fibrosis model, the HELF cell line was selected and was induced with TGFβ1, followed by assessment of the expression of α-SMA and collagen type I as the markers of fibrosis. The mRNA expression of α-SMA and collagen type I was found to be upregulated in the model (Fig. 5A). The results were also further confirmed by immunofluorescence analysis (Fig. 5B,D). We followed established protocols to induce polarization THP1 cells into M2 macrophages and found that the cell surface marker CD206 was significantly upregulated in M2 macrophages compared to M0 macrophages (Fig. 5C,E) [34]. The CXCL12 and CXCR4 protein structures were predicted by AlphaFold2, and they formed close hydrogen bonds across amino acid residue sites such as CYS-218 and LYS-64 (Fig. 5F). Additionally, CXCL12 was observed to increase in HELF cells upon TGFβ1 induction, whereas CXCR4 was found to be increased during M2 polarization (Fig. 5G,H).
To further determine the direct involvement of CXCR4 in M2 macrophages, M2 macrophages were then cultured in media from TGFβ1-treated HELF cells, either with or without CXCR4 antagonist (AMD3100) or CXCR4 small interfering (siRNA). The results demonstrated that inhibition of CXCR4 leads to a significant reduction in M2 polarization (Fig. 6A). To investigate the necessity of the CXCL12/CXCR4 axis for M2 chemotaxis to fibroblasts, we collected the supernatant of TGFβ1-treated HELF and used it as chemoattractant for M2 in vitro. The outcomes from the transwell assays revealed that the presence of chemokines in the media derived from the fibrosis model can enhance the chemotactic capacity of M2 macrophages, and this migratory response can be suppressed when AMD3100 and CXCR4 siRNA are used (Fig. 6B). To further explore the underlying mechanisms responsible for chemotactic effects, M2 macrophages were introduced at different time points and concentrations of CXCL12, and the involvement of the ERK pathway was evaluated by WB analysis. We found the expression level of ERK, pERK, MEK, and pMEK increased in a time- and dose-dependent manner (Fig. 6C). To further elucidate the potential role of the CXCX12/CXCR4 axis in M2 macrophages, fibroblasts were then cocultured with M2 macrophages using transwell plates in vitro. The expression of α-SMA and collagen type I was significantly downregulated upon treatment with the ERK pathway inhibitors (PD98059 and U0126). This trend, as anticipated, could be reversed by the CXCR4 agonists (NUCC-390) (Fig. 6D,E).
Importance of APOE for the Differentiation of SPP1_RecMacs
Next, to trace the regulatory mechanisms associated with transcriptional processing for Mo_AMs, we projected the scATATCseq data (GSE214085) with our scRNAseq atlas. Following stringent quality filtering and exclusion of doublets in the scATACseq dataset, 13,345 cell profiles, specifically 6,054 immune cell profiles, were obtained using ArchR (Fig. 7C, S4). Moreover, no discernible batch effects were observed upon LSI integration (Fig. 7A,B,D,E). Through the CCA integration method, the cluster annotation of scATACseq can be accurately labelled consistently with our scRNAseq cell types. Interestingly, despite the presence of potential heterogeneity among the samples, the IPF group tended to exhibit a higher proportion of SPP1_RecMacs compared to the control group (SC308, SC318, and SC326 belong to the IPF group; SC316 and SC324 belong to the control group; Fig. 7F). The foot printing analysis, which mirrored and assessed the genome binding sites of TFs, revealed that JUNB, BACH2, FOSL2, JUND, SPIB (PU.1), and SMARCC1 in Mo_AMs exhibited stronger binding to the wide genome than TRMs (Fig. 7F–L). Moreover, we observed the deep chromatin accessibility and high expression of APOE in SPP1_RecMacs as well as TRMs (Fig. 8A,C). Through the imputed pseudotime analysis of the scATACseq data, we also successfully illustrated the upregulation of APOE expression in both the GeneScoreMatrix and the GeneIntegrationMatrix during the differential process of Mo_AMs (Fig. 8B–E), which was consistent with trajectory heatmaps in our scRNAseq data and scATATCseq data (Fig. 3D, Fig. 8I). Additionally, the motif heatmap also showed a change in chromatin accessibility of JUNB, JUND, FOSL2, SPIB (PU.1), and SMARCC1 during the trajectory, in line with our results from the analysis of the foot printing and TF activities of our scRNAseq atlas (Fig. 7H, 8L, S3A). Using a combination of motif matrix and GeneIntegrationMatrix, we ultimately identified that SMAD2 and PPARγ could be the potential drivers during the Mo_AM differentiation (Fig. 8J–K).