Newborn C57BL/6J mice (12 h old, both sexes) and their dams were placed in either room air (21% O₂) or 95% O₂ for 3 days using an A-chamber system (BioSpherix, Redfield, NY, USA). Dams were rotated between hyperoxia and room air every 24 hours to prevent maternal injury. Pups were subsequently maintained in room air until pnd7 or pnd60. Each pup was treated as an independent biological replicate. Lungs from 3–5 mice originating from 2 litters were collected, with 3–4 mice assigned to each experimental group. Group assignment, data collection, and analysis were performed with the investigator blinded. Sex-specific effects were not assessed. This is because there is no sex difference in lung cellular senescence in neonatal mice after the 3-day hyperoxia [6].

For a senolytic cocktail treatment, mice received quercetin and dasatinib (2.5 and 5 mg/kg, i.p.) from pnd4 to pnd14 following neonatal exposure to room air or hyperoxia. These doses were selected based on prior reports demonstrating effective senescent-cell clearance without toxicity [6]. In a separate experiment, a PDK kinase inhibitor sodium dichloroacetate (DCA, 15 and 30 mg/kg, i.p.) was injected daily between pnd4 and pnd6. Animals were anesthetized with ketamine (75 mg/kg, i.p.) and xylazine (10 mg/kg, i.p.), prepared from stock solutions of ketamine (10%, w/v) and xylazine (2%, w/v) and diluted in sterile saline to achieve the appropriate injection volume. Mice were then cervical dislocated following confirmation of deep anesthesia and complete loss of reflexes, with subsequent removal of vital organs.

Lung morphometric analyses were performed on hematoxylin and eosin (H&E)-stained mouse lung sections. Non-lavaged lungs were gently inflated with 1% low-melting point agarose at a constant pressure of 25 cm H₂O and subsequently fixed in 4% neutral-buffered paraformaldehyde. After fixation, lung tissues were paraffin-embedded and sectioned at a thickness of 4 µm using a rotary microtome (Leica Biosystems, Deer Park, IL, USA). Radial alveolar count (RAC) was determined by drawing a perpendicular line from the center of a respiratory bronchiole to the distal acinus, defined by the pleural surface or the nearest connective tissue septum, and counting the number of alveolar septa intersected by this line. At least three measurements were obtained per animal to ensure reliable quantification.

Paraffin-embedded lung sections underwent deparaffinization and heat-induced antigen retrieval prior to immunofluorescence staining. Samples were incubated overnight at 4 ℃ with primary antibodies against Pdh E1 subunit $\alpha$1 (Pdha1, ab168379, Abcam, Cambridge, MA, USA, 1:100 dilution), NOS2 (PA1-036, ThermoFisher, Waltham, MA, USA, 1:50 dilution), Syk (PA5-96063, ThermoFisher, 1:100 dilution), transferrin receptor (ab84036, Abcam, 1:100 dilution), or CD68 (ab283654, Abcam, 1:100 dilution), co-staining with lamin B1 (MA1-06103, ThermoFisher, 1:100 dilution) and p21 (MA5-31479, ThermoFisher, 1:100 dilution) as senescence biomarkers. After incubation with goat-anti mouse (A-11001, ThermoFisher) or goat-anti-rabbit (A-11012, ThermoFisher) secondary antibodies (1:5000, 1 h, room temperature), sections were mounted in 4^′,6-diamidino-2-phenylindole (DAPI)-containing hard-set medium (Vector Labs, Newark, CA, USA). To assess microvascular density, von Willebrand factor (vWF) immunostaining was performed in the lung. Images were acquired using a Nikon fluorescence microscope (Melville, NY, USA).

Pdh activity was assessed in protein lysates prepared from snap-frozen lung tissue using a commercially available kit (ab287837; Abcam). Approximately 50 mg of frozen lung tissue was homogenized in 300 µL of ice-cold assay buffer and incubated on ice for 10 min. The homogenate was then centrifuged at 10,000 $\times$g for 5 min, and the supernatant was collected. Each sample was transferred to a 96-well plate and brought to a final volume of 50 µL per well with assay buffer. For the positive control, 10 µL of the supplied control material was added to the wells and adjusted to 50 µL with the assay buffer. The enzymatic reaction was initiated by adding 50 µL of the reaction mixture to each well, followed by incubation for 30 min at room temperature in the dark. Pdh enzymatic activity was quantified by measuring absorbance at 450 nm using a Cytation 5 imaging microplate reader (BioTek, Winooski, VT, USA).

Previously published single-cell RNA-seq data (GEO: GSE207866) from senescent lung cells isolated at pnd7 (SD7) were used for analysis [6]. Single-cell datasets from age-matched air-exposed (AirD7) and hyperoxia-exposed (O2D7) mice served as reference controls.

Data processing was carried out in Seurat v4.1.1 (Satija Lab, New York, NY, USA) [10]. CellRanger count matrices (10$\times$ Genomics) were imported using Read10X() to create individual Seurat objects. Cells expressing fewer than 200 genes, or genes detected in fewer than three cells, were excluded. Additional filtering removed cells with fewer than 700 or more than 8000 detected genes or $>$5% mitochondrial transcripts. Quality control was performed as previously described [6].

Normalization and scaling were performed with SCTransform using default parameters. Integration of Seurat objects employed canonical correlation analysis with SelectIntegrationFeatures(), PrepSCTIntegration(), FindIntegrationAnchors(), and IntegrateData(). Principal component analysis (first 50 PCs) was used to build a K-nearest neighbor graph (K = 20), followed by Louvain clustering using FindClusters(). Cluster resolution was guided by clustree v0.5.1 [11]. UMAP visualization was generated using the first 50 CCA embeddings. Cell types were assigned based on established marker genes (Table 1).

The same workflow was repeated iteratively to derive subclusters within each major lineage. Immune cells were isolated from the global dataset, and macrophages were distinguished from monocytes and dendritic cells. This approach yielded seven macrophage subsets (518 cells) for downstream analyses. Marker genes for each macrophage subset were identified using the cosg() function [12]. Figures were generated with SCpubr v2.0.0, scplotter v0.1.1, and fmsb v0.7.6 in R v4.4.1 using RStudio v2023.12.0 (RStudio/Posit, Boston, MA, USA).

Gene sets representing GO biological processes, KEGG pathways, Reactome, PID, BioCarta, and WikiPathways were obtained from the Molecular Signatures Database (v7) (Broad Institute, Cambridge, MA, USA) [13]. Pathway comparisons were carried out using compare_pathways() in SCPA v1.6.2 (Broad Institute, Cambridge, MA, USA) (parameters: min_genes = 15; max_genes = 500) [14]. Macrophage populations from AirD7, O2D7, and SD7 were analyzed. Visualization and post-processing were performed using Seurat v4.4.0 [10] and ggplot2 v3.5.1 (RStudio/Posit Software, Boston, MA, USA). Pathways with $|$Fold change$|$ $>$5 and a Benjamini-Hochberg (BH)-adjusted p-value 0.01 were selected for subsequent visualization.

Single-cell trajectory analysis of seven macrophage clusters was performed using Monocle2 (v2.32.0) (Trapnell Lab, University of Washington, Seattle, WA, USA) [15]. Prior to cell ordering, the differentialGeneTest function was applied to identify differentially expressed genes (DEGs) across clusters, and genes with a q value 0.01 were selected for pseudotime ordering. Cell trajectories were then constructed using the Discriminative Dimensionality Reduction with Trees (DDRTree) algorithm with default parameters to infer transitional relationships among macrophage clusters. Branch Expression Analysis Modeling (BEAM) was performed on pseudotime-ordered data to identify genes associated with branch fate decisions (Benjamini-Hochberg (BH)-adjusted p value 1 $\times$ 10^-5). Trajectory visualizations and heatmaps were generated using the plot_cell_trajectory, plot_complex_cell_trajectory, and plot_genes_branched_heatmap functions. Upregulated genes from branch-specific BEAM analyses were subsequently subjected to pathway enrichment analysis using the DAVID database (https://davidbioinformatics.nih.gov/), and pathways with a BH-adjusted p value 0.05 were selected for downstream visualization.

AUCell (Version 1.26.0) (Bioconductor, Buffalo, NY, USA) [16] was employed to assign metabolic pathways activity scores in the single-cell RNA data. Initially, a ranking of selected pathways genes was built based on the single-cell expression matrix using the AUCell_buildRankings() function with default parameters. Subsequently, the area under the curve (AUC) was calculated using the top 5% of genes in the ranking using the AUCell_calcAUC() function. Cells expressing a higher proportion of genes within the gene set received higher AUC values. The AUC score for each cell was mapped onto violin and box plots for visualization using the ggplot2 v3.5.1. Finally, the Wilcoxon rank-sum test for multiple comparisons was performed using the stat_compare_means function from the ggpubr v0.6.0 (RStudio/Posit Software, Boston, MA, USA). A p-value 0.05 was considered statistically significant.

Data are presented as mean $\pm$ SEM. Statistical analyses were performed using GraphPad Prism 9 (GraphPad Software, San Diego, CA, USA) with blinding whenever feasible. Group comparisons were made using unpaired Student’s t-tests with Welch’s correction. For multiple comparisons, the statistical significance of the differences was evaluated by using one-way ANOVA followed by Tukey’s post-test to specifically compare indicated groups. A p-value 0.05 was considered statistically significant.

Using canonical marker genes (Table 1), we first assigned all single cells in the AirD7, O2D7, and SD7 datasets to epithelial, endothelial, immune, or mesenchymal lineages (Table 1, level 1). Immune populations corresponded to clusters 5, 6, 12, 14, 16, 18, 19, 21, 25, 26, 30, 32, and 33 (Fig. 1A). Among the 493 SD7 cells, 73.63% were immune cells, with smaller fractions of epithelial (8.11%), endothelial (9.74%), and mesenchymal (8.52%) cells (Table 2). Thus, immune populations constituted the largest senescent pool.

Fig. 1.

Immune cell and coarse macrophage distributions at postnatal day 7 (pnd7). Reanalysis of public lung scRNA-seq datasets (GSE207866) was performed using lung single-cell suspensions without C12FDG sorting from air-exposed (AirD7) and hyperoxia-exposed (O2D7) mice, as well as C12FDG-sorted cells from the hyperoxia group (SD7). (A,B) Dot plots display cell type identification, including Ptprc (CD45)-positive immune clusters (A) and coarse macrophage/monocyte subpopulations (B). Dot size indicates the proportion of cells expressing the gene within a cluster, and color intensity represents expression levels.

To further resolve these immune cells, all immune clusters were combined and reclassified using lineage-specific markers (Table 1, level 2). Eight immune cell types, such as B cells, NK cells, T cells, mast cells/basophils, neutrophils, monocytes, dendritic cells, and macrophages, were identified. Monocytes, dendritic cells, and macrophages accounted for 63.3% of immune cells overall (within clusters 0, 1, 3, 4, 6, 7, 11, and 15), and 97.5% of immune cells in the SD7 group (Fig. 1B; Table 3). Clusters 3 and 6 expressed overlapping signatures of macrophages with monocytes or dendritic cells. These clusters were combined with other myeloid clusters and refined using additional markers (Table 1, level 3).

This analysis produced eleven initial subclusters (Fig. 2A). Clusters 2, 3, 5, 6, 7, 8, 9, and 10 showed clear macrophage identity, while clusters 1 and 4 contained mixed dendritic cell and macrophage markers. Subclustering of clusters 1 and 4 allowed extraction of macrophages, which were integrated with the eight pre-defined macrophage clusters (Fig. 2B). This produced seven final macrophage subsets (clusters 0–6), totaling 518 macrophages: 116 in AirD7, 77 in O2D7, and 325 in SD7 (Tables 4,5). In SD7, macrophages represented 65.90% of all senescent cells, confirming that macrophages are the dominant senescent population after neonatal hyperoxia. We then used these seven clusters of macrophages for further analyses.

Fig. 2.

Identification of macrophage populations. (A,B) Dot plots depict cell type identification of macrophages from the total coarse macrophage population (A) and the selection of macrophages from mixed subpopulations in clusters 1 and 4 (B). Dot size reflects the percentage of cells expressing each gene, and color intensity represents expression level.

Neonatal hyperoxia drives inflammatory injury with a shift toward M1 polarization [17], but its effects on senescent macrophages are unclear. Polarization analysis of the 325 SD7 macrophages using canonical markers (Table 1) showed that clusters 0, 1, 5, and 6 displayed M1 signatures; cluster 4 aligned with M2; and clusters 2 and 3 exhibited mixed M1/M2 characteristics (Fig. 3A–C). Although total macrophage numbers were reduced in O2D7 compared with AirD7, the distribution of M1, M2, and mixed phenotypes remained similar across groups (Table 6). M1 cells were the dominant subtype in all conditions. Lamin b1 loss is a senescence-associated biomarker [18]. Immunostaining showed that NOS2 (M1 marker) was enriched in lamin b1 negative macrophages at pnd7, whereas the M2-associated protein Syk was reduced in the hyperoxic group compared to air control (Fig. 3D). These results suggest that senescent macrophages formed after hyperoxia predominantly adopt an M1-like inflammatory state.

Fig. 3.

Classification of M1, M2, and mixed M1/M2 macrophages. (A) Dot plot shows expression of curated M1, M2, and mixed M1/M2 genes in total macrophages. The size of the dot represents the percentage of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. Colored boxes are used to highlight marker gene sets that define different macrophage phenotypes and to guide interpretation of the dot plot. (B) Cluster assignment of M1, M2, and mixed M1/M2 macrophages. (C) UMAP plot illustrates the spatial distribution of these macrophages, with cells colored by type; ring plots indicate proportions in each group. (D) Dual immunofluorescence staining of lamin b1 with NOS2 or Syk at pnd7 in hyperoxia-exposed lungs. lamin b1^–/NOS2⁺ and lamin b1^–/Syk⁺ signals were quantified between air and hyperoxia groups. White arrows indicate representative macrophages showing altered NOS2 or SYK expression and reduced lamin b1 signal. Inset depicts a magnified view of the boxed region to highlight the immunofluorescent signal. Scale bar: 50 µm for the main panel, and 10 µm for the inset. N = 3; ^*p 0.05.

Because lung macrophages occupy both alveolar and interstitial compartments, we next assessed whether senescence preferentially affects specific subsets. Among the 325 SD7 macrophages, marker-based classification identified alveolar, interstitial, recruited/monocyte-derived, proliferative, and non-classical interstitial macrophages (Table 1). Clusters 0, 1, and 3 corresponded to alveolar macrophages, while clusters 2, 4, 5, and 6 mapped to interstitial, recruited/monocytes-derived, proliferative, and non-classical interstitial macrophages, respectively (Fig. 4).

Fig. 4.

Functional- and tissue-specific macrophage subsets. (A) Dot plot shows expression of curated genes for functional and tissue-specific macrophage subsets. The size of the dot represents the percentage of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. Colored boxes are used to highlight marker gene sets that define different macrophage phenotypes and to guide interpretation of the dot plot. (B) Cluster allocation of alveolar macrophages (AlvMac), interstitial macrophages (IMac), non-classical interstitial macrophages (ncIMac), recruited/monocyte-derived macrophages (recMac/Mo-Mac), and proliferative macrophages (pMac). (C) UMAP plot illustrates the distribution of these macrophages, with ring plots showing subtype proportions in each group.

In O2D7 mice, alveolar macrophage proportions were decreased and interstitial macrophages increased relative to AirD7 (Table 7). The proportions of recruited/monocyte-derived macrophages, proliferative macrophages, and non-classical interstitial macrophages were largely comparable between the AirD7 and O2D7 groups, with no substantial differences observed. In contrast, senescent macrophages in SD7 were predominantly alveolar (62.8%), with reduced representation of interstitial and proliferative subsets compared with the AirD7 group (Table 7). These results indicate that neonatal hyperoxia reshapes the macrophage landscape and that alveolar macrophages are the most susceptible to senescence.

To investigate potential transcriptional state transitions among the seven macrophage clusters, we performed Monocle2 pseudotime trajectory analysis using all clusters (Fig. 5A). The inferred trajectory revealed two major branch points, partitioning cells into five distinct transcriptional states. State 1 mainly comprised clusters 0, 1, 3, and 5; State 2 predominantly consisted of cluster 4; States 3 and 4 were primarily composed of cluster 2; and State 5 mainly included clusters 0, 1, and 5 (Fig. 5B). Based on biological relevance and gene expression characteristics, recMac/Mo-Mac (cluster 4) was annotated as the root state for pseudotime ordering [19, 20]. Accordingly, State 2, which predominantly comprised cluster 4 cells, was designated as the putative early transcriptional state in the inferred trajectory (Fig. 5C).

Fig. 5.

Pseudotime analysis identifies divergent fate programs in senescent macrophages. (A–D) Cell ordering along the differentiation trajectory displayed by clusters (A), branch states (B), pseudotime states (C), and experimental groups (D). (E) Trajectory tree structure illustrating the distribution and relative abundance of macrophage clusters across each branch among the AirD7, O2D7, and SD7 groups. S denotes State. (F) Heatmap of all significantly altered genes identified by branched expression analysis modeling (BEAM) in Monocle at branch point 1. C1–C3 are self-defined gene modules, where “C” denotes Cluster. Specifically, C1 represents a group of genes that are gradually upregulated during differentiation from branch point 1 toward State 3; C2 represents a group of genes that are gradually upregulated during differentiation from branch point 1 toward State 4; and C3 represents genes expressed in the initial cells prior to differentiation toward States 3 and 4. (G) Bar plot of curated significantly enriched pathways (adjusted p 0.05) derived from genes associated with the differentiation from branch 1 towards State 4.

Cells from the AirD7 and O2D7 groups were primarily distributed toward later pseudotime positions, whereas SD7 cells progressing along branch 2 were enriched at early or intermediate pseudotime stages (Fig. 5D). Fig. 5E depicts the branched trajectory, illustrating the distribution of each cluster across distinct States. The cluster 2 (IMac) cells from the AirD7 and O2D7 groups were predominantly distributed in States 3 and 4, whereas cluster 2 cells from the SD7 group remained in State 4 only along branch 1. The right side of branch 2 corresponds to State 5, which was largely occupied by cells from the SD7 group. State 5 was predominantly composed of clusters 0 and 1 (alvMac) and cluster 6 (ncIMac) in the SD7 group, whereas State 1 was primarily occupied by clusters 0, 1, 3, 5, and 6 in the AirD7 and O2D7 groups. In addition, cluster 5 (pMac) in the SD7 group was mainly distributed at early or intermediate stages of differentiation (i.e., State 2). Together, these distinct cellular distributions indicate divergent fate-associated transcriptional programs in senescent macrophages.

No significant differences in differentiation trajectories along branch 1 were observed between the O2D7 and AirD7 groups (Fig. 5D,E). Therefore, we focused on characterizing differentiation-associated and functional disparities in State 4 (IMac) cells in the SD7 group. To identify genes associated with fate bifurcation, BEAM was applied at branch point 1, revealing genes whose expression dynamics were significantly associated with branch-dependent fate decisions of cluster 2 cells. The C2 gene module (Table 8) was progressively upregulated during differentiation from branch point 1 toward State 4 (Fig. 5F), reflecting gene activation during the transition from cluster 4 (recMac/Mo-Mac) to cluster 2 (IMac) in the SD7 group. In contrast, the C2 module was downregulated during differentiation from cluster 4 to cluster 2 in the AirD7 and O2D7 groups, corresponding to the trajectory from branch point 1 toward State 3. Accordingly, genes in the C2 module were imported into the DAVID database for functional enrichment analysis. Representative significantly enriched pathways associated with differentiation toward State 4 were visualized in a bar plot (Fig. 5G). This includes inflammatory and immune-related pathways such as CCR chemokine receptor binding, chemokine-mediated signaling, TNF, NF-$\kappa$B, IL-17, and Toll-like receptor signaling, as well as cellular responses to TNF, eosinophil chemotaxis, and MAPK signaling. Collectively, these findings indicate that senescent macrophage-derived interstitial macrophages adopt a pro-inflammatory, fate-associated transcriptional program.

To identify defining features of each senescent macrophage subset, we used cosg() [12] to extract the top upregulated genes (Table 9). The predominant functional categories for each cluster were as follows: Cluster 0—innate immunity, inflammation, lipid metabolism, pentose phosphate pathway; Cluster 1—inflammatory and DNA repair programs; Cluster 2—lipid-associated macrophage markers; Cluster 3—downregulated immune activity and thermogenic genes; Cluster 4—migration and recruitment signatures; Cluster 5—cell-cycle regulation; Cluster 6—lysosomal pathways (Fig. 6A,B). Clusters 0, 1, and 3 were expanded in SD7 relative to AirD7 and O2D7 (Fig. 6C), indicating that senescent macrophages are enriched for inflammatory, metabolic, and innate immune programs.

Fig. 6.

Characterization of macrophage subclusters. (A) Dot plot displays the top ten genes identified in each subcluster using the Cell-type-specific One-vs-all Statistical Gene identification (COSG) algorithm. The size of the dot represents the proportion of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. (B) Biological processes associated with each subcluster’s marker genes. (C) Radar plot depicts the proportion of each macrophage subcluster across the three groups.

Because clusters 0 and 2 expressed metabolic signatures, we evaluated metabolic gene expression across groups. Clusters 0 and 1 expressed high levels of metabolic genes, whereas cluster 2 showed comparatively lower expression (Fig. 7A). Relative to AirD7 and O2D7, SD7 macrophages displayed higher glycolysis genes (Slc2a1, Pkm, Pfkfb3, Ldha, Hk2, Gapdh), higher pentose phosphate pathway genes (Tkt, Taldo1, Pgd, G6pdx), and lower $\beta$-oxidation genes (Hadh, Eci1/2, Echs1, Acads/Acadsb, Acaa2) (Fig. 7B). Pathway analysis confirmed enrichment of glycolysis, glutamine metabolism, and TCA cycle-related pathways, and suppression of fatty acid $\beta$-oxidation (Fig. 7C). Pdha1, a key regulator linking glycolysis to oxidative phosphorylation, was reduced in lamin b1 negative macrophages after hyperoxia compared to the air control (Fig. 7B,D).

Fig. 7.

Metabolic gene expression and pathway enrichment in macrophage subclusters. (A,B) Dot plots show metabolic gene expression in macrophage subclusters (A) and across AirD7, O2D7, and SD7 groups (B). The size of the dot represents the percentage of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. (C) Metabolic pathway enrichment was assessed using the compare_pathways() function in SCPA, with Wilcoxon tests for significance. ^****p 0.0001, while ns indicates no significance. (D) Dual immunofluorescence of Pdha1 and lamin b1 in pnd7 lungs; lamin b1^–/Pdha1⁺ intensity was compared between air and hyperoxia groups. White arrows indicate representative macrophages showing reduced Pdha1 and lamin b1 immunoreactivity. Inset depicts a magnified view of the boxed region to highlight the immunofluorescent signal. Scale bar: 100 µm for the main panel, and 10 µm for the inset. N = 3. ^*p 0.05.

Because clusters 0, 1, and 3 constituted 62.8% of SD7 macrophages, we assessed metabolic changes specifically within these clusters. In all three clusters, glycolysis and glutamine metabolism were increased, while $\beta$-oxidation genes were consistently reduced in the SD7 group compared to the AirD7 and O2D7 groups (Figs. 8,9,10). The TCA genes were unaltered in clusters 1, 2 or 3 among AirD7, O2D7 and SD7 groups (Fig. 10). However, mitochondrial complexes I, II, IV, and V genes were upregulated, whereas the complex III gene Bcs1l was decreased in the SD7 group compared to the AirD7 and O2D7 groups (Fig. 11). Together, these findings suggest that senescent alveolar macrophages undergo a metabolic shift toward glycolysis and glutamine utilization.

Fig. 8.

Metabolic profiles in alveolar macrophage subcluster 0. (A) Dot plot shows metabolic gene expression in subcluster 0 across AirD7, O2D7, and SD7 groups. The size of the dot represents the percentage of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. (B) Pathway enrichment was evaluated using compare_pathways() in SCPA; significance was assessed by Wilcoxon test. ^***p 0.001, ^****p 0.0001, while ns indicates no significance.

Fig. 9.

Metabolic profiles in alveolar macrophage subcluster 1. (A) Dot plot displays metabolic gene expression in subcluster 1 across AirD7, O2D7, and SD7 groups. The size of the dot represents the percentage of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. (B) Pathway enrichment analysis was performed using SCPA; statistical significance is indicated as ^*p 0.05, ^**p 0.01, ^***p 0.001, ^****p 0.0001, respectively, while ns indicates no significance.

Fig. 10.

Metabolic profiles in alveolar macrophage subcluster 3. (A) Dot plot shows metabolic gene expression in subcluster 3 among the three groups. The size of the dot represents the percentage of cells within a cluster expressing the indicated gene, while the intensity of expression is indicated by the color legend. (B) Pathway enrichment in subcluster 3 assessed with SCPA. Wilcoxon test was used to evaluate statistical significance among these groups. ^*p 0.05, ^***p 0.001, ^****p 0.0001, respectively, while ns indicates no significance.

Fig. 11.

Mitochondrial electron transport chain genes in alveolar macrophages. Dot plots show expression of mitochondrial electron transport chain genes in macrophage subclusters 1, 2, and 3 across AirD7, O2D7, and SD7 groups. Dot size represents the proportion of cells expressing each gene, and color intensity indicates expression level.

Using SCPA to compare GO biological processes across groups, we identified pathway alterations characteristic of senescent macrophages. Volcano plots show broad upregulation and downregulation of signal pathways in SD7 relative to AirD7 and O2D7 (Fig. 12A; Tables 10,11,12,13). Upregulated pathways in SD7 group included metal metabolism and homeostasis, amino acid metabolism, interleukin signals, immune effector process, robo receptor signaling, apoptotic signal pathway, heme oxygenase 1 regulation, and RNA process and translation compared to AirD7 (Table 10) and O2D7 (Table 11). Transferrin receptor is a membrane protein which binds to transferrin and mediates cellular iron uptake. This protein is also regulated by intracellular iron levels. Transferrin receptor expression was increased in lamin b1-negative macrophages after neonatal hyperoxia (Fig. 12B), supporting enhanced iron uptake and metabolism in senescent cells.

Fig. 12.

Pathway enrichment in macrophages. Pathway analysis was performed using compare_pathways() in SCPA. (A) Volcano plots show up- and down-regulated pathways in SD7 vs. AirD7 (top) and O2D7 (bottom). (B) Dual immunofluorescence of transferrin receptor (TR) and lamin b1 in pnd7 lungs; and its intensity was compared between air and hyperoxia groups. White arrows indicate representative lamin b1⁺/TR⁺ macrophages in lung sections. Inset depicts a magnified view of the boxed region to highlight the immunofluorescent signal. Scale bar: 50 µm for the main panel, and 10 µm for the inset. N = 3. (C) Violin plots show differences in p38 MAPK, ATM, mTOR (increased) and JAK/STAT (decreased) pathways in SD7 vs. AirD7 and O2D7; SCPA used for statistics. ^*p 0.05, ^**p 0.01, ^***p 0.001, ^****p 0.0001, respectively, while ns indicates no significance.

Upregulated pathways in SD7 group included lung morphogenesis, ROS/H₂O₂ detoxification, and Rho GTPase signals were downregulated compared with both AirD7 and O2D7 groups (Tables 12,13). Pathways including morphogenesis, development, branching, vasculogenesis, response to retinoic acid, insulin-like growth factor, and Runx2 were downregulated in the SD7 group compared with the AirD7 group (Table 12). Compared to O2D7 group, filament, cytoskeleton and cell junction organization, adhesion, phagocytosis, Eph/ephrin signaling, and DNA repair pathway were also downregulated in SD7 group (Table 13).

In addition, p38 mitogen-activated kinase (p38 MAPK), ataxia-telangiectasia mutated (ATM), and mechanistic target of rapamycin (mTOR) pathways were strongly enriched in SD7, whereas JAK/STAT signaling was reduced (Fig. 12C). Because p38 MAPK, ATM, and mTOR are central regulators of the SASP [21, 22], these findings suggest that senescent macrophages activate pathways that promote SASP production.

PDK phosphorylates the $\alpha$ subunit of E1, e.g., Ser293, leading to Pdh inactivation. To determine whether activating Pdh represses hyperoxia-induced macrophage senescence and lung injury, mice received a prototypical PDK inhibitor DCA (15 and 30 mg/kg, i.p.) daily between pnd4 and pnd6 after hyperoxia. As shown in Fig. 13A, injection of DCA increased lung Pdh activity in mice exposed to hyperoxia at pnd7. The numbers of lamin b1^–/CD68⁺ cells and p21⁺/CD68⁺ were increased by hyperoxia at pnd7, and this was decreased by DCA injection at a dose of 30 mg/kg (Fig. 13B,C). At pnd14, injection of DCA significantly decreased mean linear intercept (Lm) and increased RAC in hyperoxia-exposed mice (Fig. 13D). Furthermore, neonatal hyperoxia-induced reduction of vWF-positive vessels was attenuated by the treatment at pnd14 (Fig. 13E). Altogether, activating Pdh decreases neonatal hyperoxia-induced macrophage senescence and lung injury. The animal experiments presented in Fig. 13 were approved under Protocol 2024-003, while all remaining animal work described in this study was conducted under Protocol 21-08-0003.

Fig. 13.

Activating Pdh decreases neonatal hyperoxia-induced macrophage senescence and lung injury. C57BL/6J mice (12 h old) were exposed to air or hyperoxia for 3 days and allowed to recover in room air until pnd7 (A–C) or pnd14 (D,E). DCA (15 and 30 mg/kg, i.p.) was injected daily between pnd4 and pnd6. (A) Lung Pdh activity was measured. (B) Dual immunofluorescence of p21 and CD68 was performed in mice injected with a dose of DCA at 30 mg/kg. White arrows indicate representative CD68⁺ macrophages showing p21 induction. Inset depicts a magnified view of the boxed region to highlight the immunofluorescent signal. Scale bar: 20 µm for the main panel, and 10 µm for the inset. (C) Dual immunofluorescence of lamin b1 and CD68 was performed in mice injected with a dose of DCA at 30 mg/kg. White arrows indicate representative CD68⁺ macrophages showing altered lamin b1 expression. Inset depicts a magnified view of the boxed region to highlight the immunofluorescent signal. Scale bar: 20 µm for the main panel, and 10 µm for the inset. At pnd14, lung Lm and RAC via H&E staining (D) and vessel numbers via CD31 immunostaining (E) were evaluated. Scale bar: 50 µm. N = 4–6. ^**p 0.01, ^*⁣**p 0.001 vs. Veh/Air; ^†p 0.05, ^††p 0.01, ^†⁣††p 0.001 vs. Veh/O₂.

We previously demonstrated that dasatinib plus quercetin decreases early senescence in hyperoxia-exposed neonatal lungs [6]. To test whether this treatment also mitigates persistent injury, mice received the senolytic cocktail quercetin/dasatinib (2.5 and 5 mg/kg, i.p.) on pnd4 and pnd14 after neonatal hyperoxia. The numbers of lamin b1^–/CD68⁺ cells and p21⁺/CD68⁺ were increased by hyperoxia at pnd7, and this was decreased by injection of quercetin/dasatinib (5 mg/kg) (Fig. 14A,B). At pnd60, hyperoxia increased Lm and decreased RAC, consistent with impaired alveolarization, but both changes were significantly improved by senolytic treatment in a dose-dependent manner (Fig. 14C). Neonatal hyperoxia-induced vascular simplification, measured by reduced vWF-positive vessels, was also attenuated by the treatment (Fig. 14D). These findings indicate that reducing senescent lung cells ameliorates long-term structural lung injury after neonatal hyperoxia.

Fig. 14.

Senolytic cocktail mitigates hyperoxia-induced macrophage senescence and persistent lung injury. C57BL/6J mice (12 h old) were exposed to air or hyperoxia ($>$95% O₂) for 3 days, followed by air recovery until pnd7 (A,B) or pnd60 (C,D). Quercetin (Que, 2.5 and 5 mg/kg) and dasatinib (Das, 2.5 and 5 mg/kg) were administered intraperitoneally every other day from pnd4 to pnd14. (A,B) Dual immunofluorescence of lamin b1 and CD68 (A) and p21 and CD68 (B) was performed in mice injected with a dose of Que/Das at 5 mg/kg. White arrows indicate representative CD68⁺ macrophages exhibiting altered lamin b1 expression in lung sections. Scale bar: 20 µm. (C,D) Lung alveolar size (Lm) and RAC were assessed by H&E staining, and microvessel (100 µm) counts were determined by vWF immunofluorescence. N = 4–5. ^*⁣**p 0.001 vs. corresponding Air. ^†p 0.05, ^††p 0.01, ^†⁣††p 0.001 vs. corresponding Veh/O₂.