Data selection and DEG analysisThe selection and details of the mRNA datasets are described in Table 1. Newborn mice were exposed to hyperoxia from the first day after birth to induce BPD-like lung injury, and mice placed under room air served as controls. Due to the fact that the BPD mouse models in GSE25286 and GSE209664 were constructed under similar oxygen saturation conditions for more than 14 days, they were chosen as the subsequent exploration set and validation set, respectively. The GSE25286 contained lung tissues from 6 hyperoxia-induced BPD mice and 4 wild-type (WT) mice under room air. The GSE209664 contained lung tissues from 6 BPD mice and 6 WT mice under room air. When considering the possible incomparability in lung pathology due to different interventional methods, GSE193187 and GSE211356 were not selected as the validation sets. The flowchart of the study design and analytical procedure are shown in Fig. 1. According to the selection criteria, a total of 788 genes were found to be differentially expressed between BPD and control mice, including 456 upregulated genes and 332 downregulated genes. Volcano plots and heatmapss of the DEGs are shown in Fig. 2A,B. The enrichment analysis revealed the top 20 GO terms and KEGG pathways associated with the most DEGs with statistical significance (Fig. 2C). GO enrichment analysis demonstrated that the DEGs were mostly related to leukocyte migration, ameboidal-type cell migration, cell substrate adhesion, mitotic cell cycle phase transition and cell chemotaxis, thus suggesting the potential participation of the immune microenvironment and matrix remodelling in the pathogenesis of BPD. KEGG enrichment analysis indicated that the primary pathways involved in the pathogenesis of BPD included PI3K-Akt signalling pathways, MAPK signalling pathways, focal adhesion, regulation of actin cytoskeleton and lipid and atherosclerosis.Table 1 The selection and details of mRNA datasets in this study.Fig. 1Flow diagram of the study design.Fig. 2Differentially expressed genes (DEGs), immune cell infiltration analysis and enrichment analysis in GSE25286. (A) The volcano plot of DEGs. Red dots represent up-regulated genes and blue dots represent down-regulated genes. (B) The heatmap of DEGs. The depth of the colour represents the expression level. The deeper the colour, the higher the expression level is. (C) The bubble plot of Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. The dot size represents the count of DEGs involved and the colour represents the size of P values in the KEGG pathways. (D) The immune cell fractions in each sample of GSE25286. (E) The comparison of proportions of different immune cells between BPD and control groups.Characteristics of immune cell infiltrationWe applied the CIBERSORT algorithm to analyse the immune cell composition based on the mRNA-seq data of GSE25286. The immune cell fractions of each sample were displayed as a bar plot (Fig. 2D). Monocytes and CD4 + memory resting T cells constituted the majority of immune cells in lung tissues. As shown in the box plot (Fig. 2E), compared to those in control mice, the infiltration of M1 macrophages and resting NK cells was significantly greater in BPD mice (P < 0.05). The conversion of M2 macrophages to M1 macrophages is associated with a proinflammatory state, and increased infiltration of M1 macrophages possibly contributes to lung inflammation and developmental abnormalities in BPD.Afterwards, M1 macrophages and resting NK cells were used as the trait data for WGCNA to retrieve immune-related genes. A soft-threshold power of 8 was selected (Fig. 3A) and 6 modules were identified after clustering and merging similar module eigengenes (Fig. 3B). Correlation analysis (Fig. 3C) revealed that both BPD and differential immune cells were significantly related to the orange and yellow modules (P < 0.05). The correlation scatter plots between modules and traits are shown in Additional file 2. The genes in the orange and yellow modules subsequently intersected with the DEGs (Fig. 3D), thus yielding 172 immune-related genes (Additional file 3). When considering the important functions of macrophages compared to other immune cells in lung tissues, we focused on the relationships between these genes and macrophages in the subsequent analyses.Fig. 3WGCNA analysis using immune cells as trait data and hub gene selection. (A) Scale-free network analysis and mean connectivity analysis for various soft threshold powers. (B) The clustering dendrogram of differentially expressed genes (DEGs) constructed based on topological overlap matrix. (C) The heatmap of correlations between immune traits and modules. The coefficients and P values are shown in the box. (D) The Venn diagram showing the selection process of immune-related genes: intersection between module genes and DEGs. (E) Protein–protein interaction network of hub genes identified by imputing the immune-related genes to Cytoscape. The depth of red colour represented the scores of EPC algorithm.Hub gene identification and validationIn the PPI network of the 172 immune-related genes, topological analysis using the cytoHubba plugin and EPC algorithm revealed 15 hub genes, including Cybb, Pfn2, Papss2, Atp5a1, F7, Fpr2, Akt, CD79a, Dhfr, Clca1, Met, Sepw1, Mmp9, Ppbp and Fdps.) Cybb, Papss2, Pfn2, Atp5a1, F7, Fpr2 and Akt3 were relatively key nodes identified by EPC algorithm in the PPI network with deeper colour (Fig. 3E).First, we explored the cell source of the 15 hub genes in the scRNA-seq dataset GSE209664. Single-cell identification analysis was performed on the total 65,733 cells (38,223 cells from the BPD group and 27,510 cells from the control group). Cluster analysis, which were visualized by uniform manifold approximation and projection (UMAP), identified 25 types of cell clusters in the lung tissues and the number of cells in each cluster is shown in Additional file 4. Among the 65,733 single cells, 2492 were classified as alveolar macrophages, and 5447 were classified as other types of macrophages. We then examined the expression levels of the hub genes among the cell clusters. Among the 15 hub genes, 6 genes (Cybb, Mmp9, Fpr2, Cd79a, Papss2 and F7) were identified in the scRNA-seq dataset GSE209664. Cell source identification analysis revealed that Cybb, Fpr2 and F7 were expressed in large quantities in macrophages whether in BPD group (Fig. 4A) or control group (Fig. 4B). The expression pattern of the 3 hub genes in different cell types were visualized by using UMAP (Fig. 4C). Among the total single cells, 12,615 cells expressed Cybb (19.19%), which was the highest percentage among the 3 genes. Intuitively, Cybb and F7 were mainly expressed in macrophage clusters, and Fpr2 was mainly expressed in macrophage and neutrophils clusters. Since Cybb was the key node with highest EPC score in PPI network and had the closest relationship with macrophages, we mainly focused on the role of Cybb in subsequent analysis. We compared the expression of the hub genes between the BPD and control groups by using the scRNA-seq data. Consistent with the prediction in the mRNA dataset GSE25286, the results in GSE209664 also indicated that Cybb was generally downregulated in the BPD group (Fig. 5A,B). Comparison of Cybb expression at the single-cell level between BPD and WT mice further demonstrated that the decrease was mainly due to the fact that the average expression and percentage of Cybb in macrophages were lower in BPD mice compared to WT (Fig. 5C) mice. These results suggested that Cybb expression predominantly decreased in macrophages under hyperoxia.Fig. 4Exploring the cell source of hub genes using scRNA-seq data GSE209664. (A) Dot plots of the expression pattern of the hub genes in different types of cells in WT group. (B) Dot plots of the expression pattern of the hub genes in different types of cells in BPD group. (C) The expression pattern of the 4 hub genes Cybb, F7 and Fpr2 in all single cells shown by UMAP algorithm. Cell clusters in red line represent macrophages. The range of colour represents the expression level of genes.Fig. 5Validating expression levels of hub genes in scRNA-seq data GSE209664 identified downregulation of Cybb inmacrophages of the BPD group. (A) Dot plots of the expression pattern of the hub genes in the WT and BPD group. The size of the circles indicates the percentages of cells expressing the gene in different cell types; the colour represents the average expression levels of the gene. (B) Comparison of the expression level of Cybb between BPD and WT groups. (C) The dot plot of the expression pattern of Cybb in the BPD and WT groups.Cell–cell crosstalk among macrophages, alveolar epithelial cells and endothelial cellsThe analysis was first conducted on the merged data of the BPD and control groups. Among the signal-sending cells, type I alveolar epithelial cells (AT1) had the most interactions with other cell types in tterms of quantity and strength; among the signal-receiving cells, CAR4 + endothelial cells were the most significant in terms of strength (Additional file 5A,B). In total, there were 11 ligand–receptor pairs among macrophages, alveolar epithelial cells and endothelial cells predicted by the “CellChat” package with statistical significance (P < 0.01) (Additional file 5C). The ligand–receptor signalling with the highest communication probability was macrophage migration inhibitory factor (MIF)-(CD74 + CD44) from type II alveolar epithelial cells (AT2) to macrophages, from AT2 to AT2, from AT1 to macrophages, and from AT1 to AT2, which suggested a close relationship between macrophages and alveolar epithelial cells. The incoming and outgoing signalling pathways in different cell clusters are shown in Additional file 5D. Macrophages mainly acted as signalling receivers involved in SEMA3, MIF, GALECTIN, CXCL, CCL and COMPLEMENT signalling pathways, whereas macrophages acted as signalling senders involved in GRN and CCL signalling pathways.Afterwards, we explored cell–cell crosstalk in the BPD group and control group, respectively. As the signalling senders, macrophages only interacted with alveolar epithelial cells in WT controls, while interacted with alveolar epithelia cells and various endothelial cells in the BPD mice (Fig. 6A). AT1 and AT2 had more interaction counts in the BPD group than in the WT group, whereas Lym-endothelial cells showed the opposite trend. The interaction strength of these cells remained unchanged in general, except for the observation that the interaction strength of AT2 was increased in the BPD group (Fig. 6B). In the BPD group, the ligand–receptor signalling pathways with the highest communication probability was still MIF-(CD74 + CD44) from AT2 to macrophages, from AT2 to AT2 and from AT1 to macrophages, the probability of which was not high in the WT group (Fig. 6C). In addition, the relative importance of each cell type in MIF signalling varied between BPD and WT mice, which indicated increased involvement of macrophages, especially as signalling receivers, in BPD mice. In contrast to the major interaction between AT2 and endothelial cells in WT mice, macrophages became the main receiver of MIF signalling in BPD mice (Fig. 6D). Collectively, these results suggested that macrophages were more active and interacted more with alveolar epithelial cells in BPD mice through MIF-(CD74 + CD44) signalling than in WT mice.Fig. 6Comparison of cell–cell crosstalk patterns between BPD and controls in the scRNA-seq data. (A) Signalling to alveolar cells and endothelial cells with macrophages as senders. (B) The incoming and outgoing interaction strengths. The dot size represents the involved ligand–receptor pairs in each cell type. (C) The ligand–receptor pairs predicted via the “CellChat” package with statistical significance (P-value < 0.01). The scale colour represents the level of communication probability from minimum to maximum. (D) The network centrality measures of macrophage migration inhibitory factor (MIF) signalling. The depth of colour represents the relative importance of each cell type in MIF signalling.Validation of the hub genes and key pathways in mouse modelsThe differences in lung morphometry and body weight between the hyperoxia and room-air groups indicated that we successfully established BPD-like mouse models (Additional file 6). To further validate the roles of macrophages-related genes in BPD, the gene expression levels between the hyperoxia group (BPD) and room-air group (WT) were further compared in our mouse models. In accordance with the bioinformatic analysis, the validation experiments revealed that the mRNA expression levels of the 3 macrophage-related genes (Cybb, Fpr2 and F7) were downregulated in the BPD group (Fig. 7A). Due to the fact that M1 macrophages was significantly upregulated in BPD mice and that Cybb was the most closely related to macrophages in our previous analysis, we focused on macrophage infiltrations and Cybb expression. At the protein level, Krüppel-like factor 4 (Klf4), the reduction of which is known to indicate macrophages M1 polarization, was downregulated in the hyperoxia group. The level of the Cybb protein was lower in the hyperoxia group as expected (Fig. 7B). The full-length blots of Cybb and Klf4 are included in Additional file 7. We subsequently detected the relationships of Cybb expression patterns and macrophage infiltration via immunofluorescence. The number of M1 macrophages (marked with CD86) was evidently greater in BPD than in WT. Interestingly, we detected colocalization of Cybb and CD86 in WT, whereas a more irregular and scattered distribution pattern of Cybb independent of CD86 was observed in BPD (Fig. 7C). On average, 60.1% of the CD86(+) cells in WT expressed Cybb, which was significantly greater than the 4.6% in BPD (P < 0.05). In addition, 86.5% of the Cybb(+) cells were also CD86(+), compared to 7.3% in BPD (P < 0.05). The colocalization of Cybb and CD206 (macrophage M2 marker) also seemed to be greater in WT than BPD; however, the difference did not reach statistical significance (P > 0.05) (Fig. 7D). Collectively, the above mentioned results suggested that the downregulation of Cybb in macrophages may be closely related to M1 macrophage infiltration in BPD.Fig. 7The validation of hub genes and key pathways in BPD mice models. (A) The mRNA expression level of Cybb, Fpr2 and F7 in lungs between hyperoxia and room-air groups evaluated by qRT-PCR. (B) The protein level of Cybb and macrophage polarization indicator (Klf4) in lungs between hyperoxia and room-air groups evaluated by WB. (C) The representative images of the expression pattern of Cybb(orange) and macrophages markers [M1 phenotype labeled with CD86 (pink) and M2 phenotype labeled with CD206 (red)] in lungs between hyperoxia and room-air groups evaluated by immunofluorescence staining (× 80 magnification). Co-labeled cells are indicated with arrows. (D) Fluorescence quantification: comparing percentages of co-labeled cells in single-labeled cells and using student t test.
