Single-cell RNA sequencing data were obtained in matrix format from the Gene Expression Omnibus (GEO) under accession number GSE183682. Lung tissues were harvested from mice of normal saline (NS) and SiO₂ groups at 7 days and 56 days (mouse and modeling methods as described in Chen et al.), immediately frozen in OCT on dry ice, and stored at -80°C before further processing (15). Quality control was performed using Seurat and DoubletFinder to identify and remove multiplets and low-quality cells. The Seurat package was also used to filter out foreign cells, where a gene was considered expressed if detected in more than 3 cells, and each cell was required to have at least 200 expressed genes. Following these quality control procedures, a total of 36,122 high-quality single cells were retained for downstream analysis. Gene expression data were subsequently normalized using Seurat’s “LogNormalize” method to reduce variability in gene expression counts.

To identify highly variable genes (HVGs), Seurat’s “FindVariableFeatures” function was applied with default parameters. Principal component analysis (PCA) was then performed on the HVGs, and the top 30 principal components were retained for subsequent analysis. Uniform manifold approximation and projection (UMAP) was used to reduce the dimensionality of the data and to visualize the cellular clusters in two dimensions. Clustering of the cells was performed using Seurat’s “FindClusters” function. To correct for potential batch effects, the Harmony R package was employed to integrate datasets across different conditions.

Cell-cell interactions were analyzed using the CellChat R package, which relies on known ligand-receptor pairs to explore intercellular communication. The normalized counts of the merged samples were pre-processed using the identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData functions from the CellChatDB.human database. To identify significant ligand-receptor pairs, key CellChat functions, including computeCommunProb, computeCommunProbPathway, and aggregationNet, were applied. Interaction networks were visualized using CellChat’s netVisual_bubble function to assess the strength and significance of the interactions.

The “FindMarkers” function in Seurat was used to identify differentially expressed genes (DEG) in specific cell subgroups, with the following cutoff criteria: an average log fold change (Avg_logFC) greater than 0.25 and an adjusted p-value (p_val_adj) less than 0.05. To provide further biological insights into the DEGs, functional enrichment analysis was conducted using the ClusterProfiler R package to perform Gene Ontology (GO) analyses.

Trajectory analysis was conducted using the Monocle 2.0 R package to investigate the differentiation processes of epithelial cells. The “differentialGeneTest” function was used to identify genes associated with pseudotime.

To further explore the spatial transcriptomic data in conjunction with single-cell RNA-seq, the feature-barcode matrix was downloaded from GEO (accession number GSE183683). The “traint” function in the CellTrek R package was applied to co-embed the spatial transcriptomics (ST) data with the single-cell RNA-seq data. Subsequently, the celltrek function, using default parameters, was employed to project single cells onto the ST coordinates.

To score EMT related gene sets of each single cell in the scRNA-Seq data, the ‘AddModuleScore’ method from the Seurat package was utilized. Scores were calculated by expressions of EMT related genes noted in result section.

Male C57BL/6 mice (8–12 weeks; body weight range 20–24 g) were purchased from Cavion Experimental Animal Co., LTD., Changzhou City, Jiangsu Province, China (animal license number SCXY (Su) 2011-0003). Animals were raised under the condition of no specific pathogen (SPF), with a temperature of 22-24 °C, relative humidity of 50 ± 5%, and a day-night cycle of 12 hours of light and 12 hours of darkness (lights at 08:00). After a 7–14 day adaptation period, all mice were randomly divided into two groups of 12 mice each: the NS group (NS) and the silica exposure group (SiO₂).

The mice were then anesthetized by inhaling isoflurane and intranasal instillation of silica (20ug/ul, 80ul) diluted in saline was used (SiO₂). The animals in the NS group were given the same dose of normal saline nasal drops, and the mice were sacrificed at the corresponding time points, and then the lung tissues were collected and stored for further experiments.

All procedures followed the Guide for the Care and Use of Laboratory Animals (NIH Publication No. 8023, revised 1978) and were approved by the Institutional Animal Care and Ethics Committee of Anhui University of Science and Technology.

For immunohistochemistry (IHC) staining, paraffin-embedded lung tissue sections were first cleared of paraffin with xylene and rehydrated in a series of ethanol dilutions. To expose antigens, sections were heated in Citrate Antigen Retrieval Solution (Beyotime, P0090) for 30 minutes. Endogenous peroxidase activity was blocked by treating the slides with 3% hydrogen peroxide for 10 minutes. Permeabilization was conducted with 0.3% Triton^® X-100 (BioFroxx, 1139) at room temperature for 10 minutes. After a 1-hour block with 3% BSA-PBST, primary antibodies were applied and left on the sections overnight at 4°C. The following day, HRP-linked secondary antibodies were introduced for 1 hour at room temperature, and staining was developed using 3,3’-diaminobenzidine (DAB). Hematoxylin was used to stain nuclei, after which sections were dehydrated, coverslipped, and examined under a light microscope.

For immunofluorescence staining, paraffin sections of mice lungs were dewaxed, rehydrated in a series of ethanol dilutions, immersion in hydrogen peroxide, repairing antigen, permeabilization, and then blocking for 1 hour in 10% goat serum at room temperature. Subsequently, the primary antibody was incubated overnight in a 4°C refrigerator after staining, rewarmed at room temperature, and then stained with the secondary antibody at room temperature for 1 hour. Finally, nuclei were stained with DAPI (1:1000) for 10 minutes. All tissue sections were mounted on glass slides and covered with an anti-fade mounting medium. Images were captured using a confocal microscope (Olympus, FV-3000). The antibodies used in this experiment are listed in Table 1 .

Data are presented as mean ± standard error of the mean (SEM). Homogeneity of variances across groups was examined using the Brown-Forsythe and Bartlett’s tests. Comparisons involving two groups were performed using an independent-samples t-test. For comparisons involving three or more groups, two-way ANOVA was conducted, followed by Bonferroni’s multiple comparisons test. GraphPad Prism 9.0.0 software was used to perform the statistical analyses, where p< 0.05 was set as the level of significance for determining any significant differences.

To investigate the effects of SiO₂ exposure on lung tissues, we first performed dimensionality reduction using UMAP to visualize the distribution of different cell types across various experimental groups (NS_7d, SiO₂_7d, NS_56d, SiO₂_56d) using data obtained from the GEO under accession number GSE183682 ( Figure 1a ). We utilized Seurat for preprocessing the single-cell RNA sequencing data, including quality control and normalization. UMAP plots revealed distinct cell clusters, the SiO₂-treated group at 56 days exhibited markedly greater changes compared to other groups. In particular, significant alterations were observed not only within immune cell populations but also prominently among epithelial cells ( Figure 1b ). This marked epithelial cell response prompted us to further investigate their potential roles in mediating SiO₂-induced lung injury.

The increased number of epithelial cells in the SiO₂-treated group at 56 days suggests a possible involvement of cellular proliferation and repair processes induced by SiO₂ exposure. GO enrichment analysis of differentially expressed genes in epithelial cells from the disease groups revealed significant upregulation of genes involved in epithelial cell migration, proliferation and responses to external stimuli ( Figure 1c ). These findings indicate that SiO₂ exposure enhances epithelial cell expansion, likely through heightened proliferative and stress-responsive activities.

We then combined spatial transcriptomics with single-cell RNA sequencing, which revealed distinct patterns in the distribution and gene expression of various cell types across the tissue sections ( Supplementary Figure 1 ). In the 7-day treatment groups, both NS and SiO₂ exposures exhibited sparse distribution of epithelial cells across the tissue. By day 56, however, epithelial cells markedly increased in abundance and formed distinct clusters, especially within SiO₂-exposed tissues ( Figures 1d-g ). Additionally, spatial analysis revealed that macrophages were more abundant in the SiO₂-treated tissues at 56 days, which coincided with the increase in epithelial cells ( Supplementary Figure 1 ).

To further validate these observations, we performed immunofluorescence staining on lung tissue sections from the NS_56d and SiO₂_56d groups ( Figures 1h-i ). In the NS_56d group, Sftpc⁺ epithelial cells (green) were evenly distributed with moderate intensity, accompanied by a sparse presence of CD206⁺ M2 macrophages (red). In contrast, the SiO₂_56d group exhibited a significant increase in Sftpc⁺ epithelial cells and CD206⁺ macrophages. Quantitative analysis confirmed that both Sftpc and CD206 signals were significantly elevated in the SiO₂_56d group compared to NS_56d. Interestingly, Sftpc-positive cells were found to be localized around CD206-positive macrophages in the SiO₂-treated group, suggesting a coordinated expansion of epithelial and macrophage populations in response to SiO₂-induced injury ( Figure 1h ).

Taken together, these results revealed that SiO₂ exposure leads to significantly enhanced epithelial cell proliferation and a heightened immune response. These findings imply that SiO₂ exposure may induce a complex tissue response, involving both immune activation and epithelial cell proliferation.

Next, we aimed to investigate the relationship between different subclusters. The gene expression data show that the epithelial cells were classified into 5 distinct groups, with SiO₂ exposure significantly altering their composition ( Figures 2a-c ). The UMAP analysis ( Figure 2b ) shows that at 7 days, the NS and SiO₂ groups have similar distributions of epithelial subpopulation. However, we noticed C0 cells a marked shift in the SiO₂_56d group to dominate. This was further confirmed by the proportional distribution of epithelial subpopulations ( Figure 2c ), where C0 cells were significantly more abundant in the SiO₂_56d group compared to other groups.

Violin plots and heatmap analyses revealed distinct gene expression signatures for each cell cluster ( Figures 2d, e ). Cluster C1 cells, marked by high expression of Foxj1 and Tppp3, were identified as ciliated airway epithelial cells. C2 cells exhibited elevated Scgb1a1, a canonical marker of secretory club cells. C3 cells exhibited co-expressing Krt17 (progenitor/stress marker) and Ager (AT1 marker), suggesting that this population may represent transitional epithelial cells committed toward the AT1 lineage (16). C3 cells also expressed Cldn18, Cldn4, and Bcam, indicating roles in tight junction formation and barrier maintenance within the respiratory epithelium. C4 cells were identified as Lamp3 ⁺ and Sftpc ⁺ AT2 cells. C4 cells also expressed innate immune response genes, including Lyz1, Lyz6, Cxcl15, and Sfta2. Lastly, C5 cells, another secretory population, displayed high levels of Scgb3a2 and Scgb3a1, along with Calml3, Chad, and Aldh3a1, suggesting functions in epithelial cell development and response to wounding.

In contrast, C0 cells expressed epithelial markers such as Scgb3a2 and Sftpc but were distinguished by EMT-related genes, including Spp1, Mmp12, Ctss, and Vim. Unlike other clusters, C0 cells lacked Krt18, a marker of epithelial differentiation expressed across all other subpopulations, suggesting a mesenchymal transition (17–19). Further analysis revealed high expression of inflammatory activation markers S100a8 and S100a9, alongside immune- and inflammation-related genes such as Ccl6, Cxcl2, and Il1b. These genes, linked to chemotaxis and immune cell recruitment (5, 20–22), collectively indicate that C0 cells actively contribute to immune modulation and inflammatory responses in the SiO₂_56d group.

To map these dynamics over time, we employed spatial transcriptomics. At 7 days post-exposure, both the NS and SiO₂ groups exhibited sparse C0 cell populations ( Figures 3a, b ). By 56 days, however, C0 cells were significantly more abundant in the SiO₂ group ( Figures 3c, d ). To assess C0 cell dynamics, we examined the spatial expression of Sftpc and Ccl6 ( Figure 3c ). SiO₂_56d tissue showed elevated and localized Sftpc expression, while Ccl6 expression was markedly elevated and widespread ( Figure 3d ). Immunohistochemical staining further confirmed these findings ( Figures 3e, f ). In SiO₂_56d lungs, Spp1 expression was markedly elevated, co-occurring with a significant increase in Sftpc-positive cells compared to controls. Quantitative analysis showed significant increases in both markers’ area percentages, consistent with a partial EMT and inflammatory signatures in SiO₂-induced injury.

Together, these findings illuminate the dynamic epithelial response to SiO₂ exposure. The C0 cell population expresses epithelial markers such as Sftpc and Scgb3a2 but is distinguished by co-expression of EMT-related genes like Mmp12 and Vim, indicating an ongoing transition toward a mesenchymal phenotype. Notably, C0 cells also upregulate a suite of immune-related genes—including Ccl6, S100a8, and Cxcl2—driving immune cell recruitment and amplifying inflammatory responses. In contrast, other clusters maintain distinct epithelial identities while contributing to cellular defense. For instance, the Lamp3 ⁺ and Sftpc ⁺ AT2 cells (C4) express innate immune response genes such as Lyz1 and Cxcl15, supporting antimicrobial defense and tissue protection.

Having established the expansion and phenotypic characteristics of C0 cells, we next sought to investigate their developmental origin. To achieve this, we employed CytoTRACE analysis to assess the differentiation potential of epithelial cell subpopulations. The CytoTRACE analysis plot reveals that C0 cells are at the final stage of differentiation, representing mature or post-differentiation cells, while C1 and other cells are at an earlier, progenitor-like stage ( Figure 4a ). The stemness score shows that C1 cells exhibit the highest stemness score, with C0 the lowest ( Figure 4b ). This suggests that C1 cells may be involved in early-stage differentiation or progenitor functions, while C0 cells likely participate in later stages of differentiation or mature epithelial functions. This is further supported by pseudotime analysis of the epithelial cells, both with and without SiO₂ treatment ( Figure 4c ). For this analysis, we selected C1 cells as the pseudotime starting point, revealing that C0 cells represent the most differentiated state at the end of the trajectory. These findings align with our earlier observation that Krt18 was absent in C0 cells, further suggesting that they occupy a more advanced and differentiated state.

Next, pseudotime trajectory analysis and differential expression analysis revealed a clear distinction between C0 and C1 cells, as visualized in the pseudotime heatmap ( Supplementary Figure 2 ) and volcano plot ( Supplementary Figure 4 ). C0 cells were characterized by high expression of immune response genes, such as Ccl6, S100a8, S100a9, and Il 1b, which are involved in inflammation and immune signaling. Additionally, C0 cells exhibited elevated levels of EMT markers such as Spp1 and Vim, suggesting their active involvement in tissue repair or remodeling. Over the pseudotime trajectory, we observed significant downregulation of genes associated with epithelial functions, including tight junction components (Cldn4, Cldn18), an undifferentiated epithelial marker (Krt7), canonical surfactant genes (Sftpa1, Sftpb, Sftpd), and secretory epithelial markers (Scgb1a1, Scgb3a1, Scgb3a2). These shifts indicate that C0 cells diverge from traditional epithelial roles, reinforcing their specialization in immune modulation and mesenchymal transition during SiO₂-induced injury.

To further explore the spatial distribution of C0 and C1 cells, we performed spatial transcriptomics on the NS_56d and SiO₂_56d groups. Intriguingly, the data showed that C0 and C1 cells are consistently juxtaposed in both treatment groups, regardless of SiO₂ exposure ( Figures 3c, d ). This spatial proximity, coupled with pseudotime trajectory analysis ( Figure 4d ), suggests a potential developmental relationship between these clusters. Specifically, C0 cells exhibit elevated expression of EMT markers like Vim, Spp1, and Mmp12, while C1 cells retain ciliated epithelial traits marked by Foxj1, hinting that C0 cells may differentiate from C1 cells along a trajectory involving mesenchymal transition. We propose that under SiO₂-induced stress, C0 cells emerge from C1 precursors and subsequently undergo partial EMT, as evidenced by their upregulation of immune genes (Ccl6, S100a8) and EMT markers, enabling their migration and expansion to support tissue remodeling and inflammation in response to chronic injury.

We next investigate the signaling pathways and intercellular communication profiles of C0 cells. Firstly, using the CellChat analysis, we identified a substantial increase in both the number and strength of cell-cell interactions in the SiO₂_56d group compared to the NS_56d group. C0 cells in the disease group exhibited significantly enhanced interaction strength and a greater number of interactions with immune cells, particularly macrophages and neutrophils, with a smaller portion directed towards cycling cells ( Figures 5a–d ).

Further analysis of ligand-receptor interactions showed a marked increase in relative information flow for APP, SPP1, SCGB3A2, and GRN signaling in the SiO₂_56d group ( Figures 5e–g ). These pathways were categorized by their target cells: SPP1 and APP predominantly engaged macrophages and neutrophils, while SCGB3A2 and GRN mainly interacted with cycling cells ( Figure 5h ). Expression of SPP1 and APP correlated with pro-inflammatory responses, while SCGB3A2 and GRN aligned with anti-inflammatory and proliferative functions (1, 5, 23). These patterns suggest C0 contributes to both immune activation and tissue repair in SiO₂-induced injury.

Subsequent analysis of spatial transcriptomics data between cycling cells and C0 cells ( Figures 6a-c ) revealed that, as expected, Scgb3a2 and Grn were predominantly expressed in C0 cells, while their corresponding receptors, Sort1 and Marco, exhibited moderate degree of expression in cycling cells. Notably, both pairs of ligands and receptors were significantly increased in the disease group. Similarly, the ligands SPP1 and APP, along with the receptors CD44 and CD74, were observed to be expressed in C0 cells and macrophage populations, respectively ( Figures 6d-f ). Furthermore, ligand genes Spp1 and App were significantly upregulated in the SiO₂_56d samples compared to the NS_56d groups ( Figures 6d-f ), while their ligands, CD44 and CD74, showed a milder upregulation. Additionally, spatial transcriptomics data highlighted a significant rise in CD44 expression in neutrophils ( Supplementary Figure 3 ). Collectively, these findings indicate that C0 cells amplify ligand-driven signaling, via SCGB3A2, SORT, SPP1 and APP, in SiO₂-induced injury, likely enhancing immune cell recruitment and tissue repair through interactions with macrophages and cycling cells.

Previous analyses revealed shifts in epithelial cell populations, with C0 cells emerging as key mediators of immune responses and tissue remodeling after SiO₂ exposure. Interestingly, these cells exhibit mesenchymal traits, suggesting EMT. To explore this, we analyzed EMT-related gene expression along the pseudotime trajectory.

We observed distinct temporal patterns in EMT markers. Mesenchymal genes, including Vim, Fn1, Zeb1, and Dab2, were upregulated, while canonical epithelial markers Epcam and Cdh1 were downregulated as cells progressed toward the C0 phenotype ( Figure 7a ). These changes peaked in mature C0 cells, indicating a progressive but incomplete EMT. Since C0 cells also expressed differentiated epithelial markers like Sftpc and Scgb3a2, this suggests partial EMT rather than a full mesenchymal transition.

To quantify EMT dynamics, we used the AddModuleScore function to calculate an EMT index, which increased along the pseudotime trajectory, confirming a partial EMT from C1 progenitors to C0 cells ( Figure 7b ). UMAP visualization showed elevated EMT scores in C0 cells, confirming their distinct mesenchymal-like phenotype ( Figure 7c ).

To investigate drivers of this process, we examined cell-cell interactions using CellChat. Neutrophils exhibited strong ligand-receptor communication with C0 cells via neutrophil-derived ThBS1 interacting with Cd47 and Cd36 receptors on C0 cells ( Figure 7d ). ThBS1 is known to promote EMT and epithelial plasticity in cancer and tissue regeneration (24, 25).

These findings propose a model where neutrophil-derived ThBS1 drives partial EMT in epithelial cells, promoting differentiation of C1 progenitors into inflammatory, mesenchymal-like C0 cells that retain epithelial markers. This immune-epithelial crosstalk underscores epithelial-mesenchymal plasticity as a dynamic process in SiO₂-induced lung remodeling.