Transcriptional Characterization of Bronchoalveolar Lavage Fluid Reveals Immune Microenvironment Alterations in Chemically Induced Acute Lung Injury

Abstract

Purpose: Chemically induced acute lung injury (CALI) has become a serious health concern in our industrialized world, and abnormal functional alterations of immune cells crucially contribute to severe clinical symptoms. However, the cell heterogeneity and functional phenotypes of respiratory immune characteristics related to CALI remain unclear.

Methods: We performed scRNA sequencing on bronchoalveolar lavage fluid (BALF) samples obtained from phosgene-induced CALI rat models and healthy controls. Transcriptional data and TotalSeq technology were used to confirm cell surface markers identifying immune cells in BALF. The landscape of immune cells could elucidate the metabolic remodeling mechanism involved in the progression of acute respiratory distress syndrome and cytokine storms. We used pseudotime inference to build macrophage trajectories and the corresponding model gene expression changes, and identified and characterized alveolar cells and immune subsets that may contribute to CALI pathophysiology based on gene expression profiles at single-cell resolution.

Results: The immune environment of cells, including dendritic cells and specific macrophage subclusters, exhibited increased function during the early stage of pulmonary tissue damage. Nine different subpopulations were identified that perform multiple functional roles, including immune responses, pulmonary tissue repair, cellular metabolic cycle, and cholesterol metabolism. Additionally, we found that individual macrophage subpopulations dominate the cell-cell communication landscape. Moreover, pseudo-time trajectory analysis suggested that proliferating macrophage clusters exerted multiple functional roles.

Conclusion: Our findings demonstrate that the bronchoalveolar immune microenvironment is a fundamental aspect of the immune response dynamics involved in the pathogenesis and recovery of CALI.

Keywords: chemically induced acute lung injury; heterogeneity; immune microenvironment; phenotypes; phosgene-inhalation; single cell RNAseq.

Graphical abstract

Figure 1

Continued.

Figure 1

Annotation of cell types by scRNA-seq in CALI (Gas) and normal control (NC) samples. (A) ScRNA-seq data represented integrated UMAP in the Gas and NC groups and the individual sample distribution. (B) Dot plot showing the gene expression levels in each cell type, with brightness indicating the log-normalized average expression and circle size indicating the expression percentage. (C) Percentage bar graph showing the distribution of cell types in each sample. (D) Dot panels showing log10(p-value) from enrichment analysis of KEGG pathways in all cell types. (E) Heatmap of -log10(p-value) from the enrichment analysis of GO classification. Each column represents the expression values for each cell. (F) Heatmap of the overall differential cross-talk in the Gas and NC group. Circle network diagram of significant cell-cell interaction pathways. The arrows and edge color indicate direction (ligand: receptor), while edge thickness indicates the sum of weighted paths between populations. (G) Heatmap of the weight of receptor-ligand interactions in cells in the Gas and NC group. (H) All top 16 TF regulatory activities across all cell types. The darker the color, the stronger the regulatory activity.

Figure 2

Continued.

Figure 2

scRNA-seq of DCs and cycle_AMs in BALF from the Gas group. (A) UMAP plots showing normalized expression of known markers in DCs. (B) Volcano plot depicting log2(fold-change) and -log10(p-value) between the two sampled time points. Genes are color-coded according to -log10(p-value). Brighter colors represent a significant differential expression. (C) GO enrichment analysis illustrating increased immune activity and cell activation. (D) Bar plot showing the enriched KEGG pathways of DEGs in DCs. Bars indicate the ratio of up- and down-regulated genes. (E) UMAP plots showing normalized expression of known markers in cycle_AMs. (F) Volcano plot depicting log2(fold-change) and -log10(p-value) between the two sampled time points. (G) Dot plot showing enriched GO terms of DEGs in cyc_AMs. The size of the dots is proportional to the number of DEGs in the GO term, while the color corresponds to the –log10(p-value) of the enrichment. Selected top terms are visualized. (H) KEGG pathway analysis illustrates increased DNA replication, cell cycle spliceosome, and mismatch repair in cyc_AMs.

Figure 3

Analysis of macrophage subpopulations. (A) Violin plots of normalized expression values for canonical cell-specific marker genes for AM clusters. (B) ScRNA-seq represented the integrated UMAP along with the individual sample distribution. (C) Dot plot showing the gene expression levels in each sub-cell type, with brightness indicating the log-normalized average expression and circle size indicating the expression percentage. (D) Percentage bar graph showing the distribution of cell types in each sample. (E) Violin plots indicating between-group differences in the expression levels of canonical markers. The p-values are calculated using a two-sided Wilcoxon signed-rank test from a theoretical null distribution. (**P < 0.005, ***P < 0.001, ****P < 0.0001).

Figure 4

Continued.

Figure 4

Pseudotime analysis of macrophages in BALF. (A) Trajectory of the macrophage subtypes in the reduced dimensional space. The points are colored according to cell types, and the two branches are labeled. (B) Pseudotime of macrophages in the reduced dimensional space. The points are colored according to the gradient of the pseudotime. (C) Distribution of macrophages from different groups along the pseudotime trajectory. The points are colored according to the groups. (D) Box plot of pseudotime distribution between groups. (E) Branch heatmap visualizing changes for all genes that are significantly branch-dependent. Columns represent points in pseudotime, rows represent genes, and the pseudotime begins in the middle of the heatmap. Genes were clustered based on their expression pattern, with four gene clusters being detected. KEGG enrichment analysis of genes in cluster 1 (F), cluster 2 (G), cluster 3 (H), and cluster 4 (I).

Figure 5

Continued.

Figure 5

(A) Dot panels showing the -log10 (p-value) from the enrichment analysis of KEGG pathways for all cells. (B) Heatmap of the -log10(p-value) of GO enrichment for all cells. Violin plots indicate between-group differences in the gene-related activity of steroid biosynthesis (C), cell polar (D), and autophagy (E). (F) Cell-cell ligand-receptor network analysis. Circle network diagram of significant cell-cell interaction pathways. The arrows and edge color indicate direction (ligand: receptor), while edge thickness indicates the sum of weighted paths between populations. Heatmap of differential interactions between the Gas and NC groups. (G) Heatmap showing the regulatory activity of top TFs across the cell types. The darker the color, the stronger the regulatory activity. (H) Expression activity of selected functional gene sets across cell types was calculated using AUCell. (*P < 0.05, **P < 0.005, ****P < 0.0001, NS, no significant).