CD8A is a Promising Biomarker Associated with Immunocytes Infiltration in Hyperoxia-Induced Bronchopulmonary Dysplasia

Introduction

Bronchopulmonary dysplasia (BPD) is known as a devastating chronic lung disease observed commonly in premature infants with frequent and severe respiratory infections and a lifelong impairment in pulmonary functions.^1,2 Approximately 10–30% premature infants who were born with <1000 g birthweight before 30 weeks’ gestation would have the potential to develop BPD.³ Mounting evidences have demonstrated that BPD is associated with considerable morbidity and mortality in the newborn period, childhood, and adulthood.^4,5 In addition to the respiratory system, BPD also shows adverse impacts on the neurodevelopment and other organs throughout the body.⁶ Therefore, it is of great importance to develop novel biomarkers for tracking the onset and progress of BPD and provide powerful and reliable interventions for BPD.

The pathogenesis of BPD is complex and remains largely unknown. Recent studies revealed that immune response plays an important role in the pathological process of BPD.⁷ A genomic study conducted by Revhaug et al reported that the immune system-related genes were impacted in a murine experimental model of BPD.⁸ Hyperoxia could evoke the innate immune response by promoting the generation of chemokines and cytokines and recruiting the immune cells to the pulmonary tissue in BPD model subsequently.⁹ Moreover, neonatal hyperoxic exposure can lead to the proinflammatory M1 phenotype in the alveolar macrophages.¹⁰ A recent study has further identified a set of differentially expressed genes and the dysregulated pathways associated with T cell activation in BPD subjects.¹¹ In recent years, bioinformatics approaches have been attempted to unveil the underlying pathways involved in BPD development. Angiogenesis and inflammation were considered as important biological changes that contributed to the BPD development, and the expression of interleukin 6 mRNA was observed pronouncedly increased in early BPD.¹² Some potential candidate genes (NOS2, MMP1, CRP, and LBP) in relation to BPD have been discovered in preterm newborns by exome sequencing.¹³ The down-regulation of EGR1, JUN, BTG2, and FOS has also been reported to be associated with BPD, and may serve as biomarkers for early BPD diagnosis.¹⁴ However, although multiple biomarkers have been identified for the early diagnosis of BPD, the underlying mechanism of BPD is still not fully elucidated. Many studies reported that immune cell infiltration is closely related to the etiology of different clinical outcomes, such as tumors, ulcerative colitis, psoriasis, arthritis, myocarditis, and atherosclerosis.^15–20 The role of immune cell infiltration in the BPD is, nevertheless, yet to be determined. Thus, it would be of great interest to investigate how the immune cell infiltration results in the pathological processes that contribute to lung injury in BPD and to examine the differences in the composition of various immune cells.

In the present study, we aimed to identify immune-related genes (IRGs) for BPD and explore the role of immune infiltration in BPD. First, we comprehensively applied the IBERSORT, MCPcounter, and ESTIMATE algorithms to characterize the immune infiltration landscape of BPD. Then, differential expression analysis and weighted gene co-expression network analysis (WGCNA) were used to screen immune-related modules and core IRGs. Additionally, the relationship between immune cell infiltrations and core IRGs was examined to further understand the related immune mechanism in the etiology of BPD. This study may offer a better understanding of the immune response in hyperoxia-related BPD and provide important insights into further research work or clinical practice.

Methods

Acquisition and Preprocessing of Gene Expression Omnibus (GEO) Data

The schematic diagram of this study is presented in Figure 1. The gene expression datasets of BPD (GSE32472, GSE25286, GSE25293, and GSE99633) were downloaded from the GEO database. GSE32472 dataset was applied to identify IRGs, which contained blood samples from 25 infants without BPD (no oxygen therapy) and 50 infants with severe BPD (oxygen therapy > 28 days), respectively. The GSE25286, GSE25293, and GSE99633 datasets were used to validate the results of bioinformatics analysis, and the BPD was induced by hyperoxia. Detailed information of the datasets is presented in Table 1. We used the Affy package of R to perform normalization and background correction.²¹ The “sva” package was used to perform batch normalization.²²

Table 1 The Detailed Descriptions of the Datasets in the Present Study

Figure 1 The schematic graph of the bioinformatics analysis.

Gene Set Enrichment Analysis (GSEA)

To identify differential enrichment pathways and biological processes involved in BPD development, we carried out GSEA based on the gene expression profile of the GSE32472 dataset. The DOSE (v.3.22.1), ggplot2 (v.3.3.3), enrichplot (v.1.16.2), and ClusterProfiler (v.4.4.4) packages were applied to visualize GSEA.^23,24 The reference gene sets of “h.all.v7.3.symbols” was downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). P-value <0.05 and false discovery rate <0.25 were set as the filter parameters.

Immune Cell Infiltrates Analysis

The gene expression profiles of BPD patients were analyzed to comprehensively assess the biological functions and characterize the landscape of immune cell infiltration in BPD. We applied CIBERSORT to analyze the gene expression profile of the GSE32472 dataset and compare the proportion of 22 immune cell types between normal and BPD groups. MCPcounter was applied to assess the proportion of 10 immune cell types between normal and BPD groups. ESTIMATE software was applied to calculate ESTIMATE Sore, Immune Sore, and Stromal Score. We used the ggplot2 package of R to visualize the difference of immune infiltrate levels between normal and BPD groups.

Identification of Differentially Expressed Genes (DEGs) and WGCNA

DEGs between normal and BPD groups of the GSE32472 dataset were identified using the limma package, and they were screened based on the cutoff criteria with p.adj < 0.05 and∣logFC∣ ≥ 0.3. Then, the DEGs expression profile was obtained, and we performed WGCNA to identify co-expression modules that are associated with immune-related pathways. We calculated the distance between each gene by Pearson’s correlation coefficient, WGCNA package was used to construct a weighted gene co-expression network.²⁵ The hierarchical clustering analysis was carried out using “hclust” function, and the “pickSoftThreshold” function was applied to calculate the soft thresholding power value. Furthermore, the co-expression network between normal and BPD groups was constructed using the “blockwiseModules” function. The hub genes in each valuable module were screened based on Module Membership (MM) > 0.8 and Gene Significance (GS) > 0.5.

Functional Annotation of the Valuable Modules

To identify the potential biological function of the valuable modules, we applied the clusterProfiler package of R to carry out the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology-Biological process (GOBP) enrichment analysis.²⁶ P-value <0.05 was considered significantly enriched.

Protein–Protein Interaction (PPI) Network and Core Differentially Expressed Immune Related Genes (DEIRGs) Identification

The immune-related genes were downloaded from the ImmPort database (https://www.immport.org/). The DEIRGs among WGCNA modules and IRGs were identified using the Venn tool. The STRING database (http://string-db.org) was used to construct the PPI network of DEIRGs, and the result was visualized using Cytoscape software (version 3.7.2). The core DEIRGs were screened using the CytoHubba plug-in of Cytoscape software. The samples of the GSE32472 dataset were divided into CD8A-low and CD8A-high subgroups based on the median expression of CD8A. The R heatmap package was used to visualize gene expression levels between CD8A-low and CD8A-high subgroups. The ggplot2 package was used to visualize the disparities of immune infiltrate levels between CD8A-low and CD8A-high subgroups.

Correlation Analysis Between Immune Cells and CD8A Expression

We assessed the correlation between CD8A expression and infiltrating immune cells using Spearman correlation analysis, and the results were visualized by the ggplot2 package of R.

Establishment of the BPD Model and CD8A Expression Validation

We purchased 30 neonatal Sprague-Dawley rats (3.5–4.5 g) from the laboratory animal center of Sichuan University. The animal experiment was approved by the Animal Ethics Committee of West China Second University Hospital (NO. 2022KX016). All animals were randomly divided into the normoxia group (n = 15, exposure to air) and hyperoxia group (n = 15, exposure to 90% oxygen).²⁷ After 14 days of hyperoxia induction, the rats were anesthetized with pentobarbital sodium (50 mg/kg), and blood samples were collected for a quantitative real-time polymerase-chain reaction (qRT-PCR) experiment. The TRIzol reagent (Invitrogen, USA) was used to extract the total RNA from blood samples. Two micrograms of purified RNA was reverse-transcribed into cDNA by Takara PrimeScript RT kit (Takara, Japan) following the manufacturer’s protocols. Then, the qRT-PCR was carried out using Brilliant SYBR Green QPCR Master Mix (Stratagene) on the ABI 7500HT PCR system (Applied Biosystems, MA, USA) based on the manufacturer’s instructions. The 2^–ΔΔCt method was used to quantify the relative expression and β-actin was used for normalization. The primers used were listed in Table S1.

Statistical Analysis

Normal distributed data were presented as mean ± standard deviation (SD). Statistical analysis was performed by using one-way analysis of variance (ANOVA). Dunnett’s test was conducted to correct for the multiple comparisons. P-value <0.05 was considered as statistically significant.

Results

GSEA of the GSE32472 Dataset

We performed GSEA on the gene expression profile of the GSE32472 dataset between normal and BPD groups to explore the potential signaling pathways related to the pathologies of BPD. As presented in Figure 2, our findings revealed that inflammation- and immune-related pathways were found mainly involved in BPD development, such as the B cell receptor signaling pathway, cytokine–cytokine receptor interaction, chemokine signaling pathway, an intestinal immune network for IGA production, primary immunodeficiency, T cell receptor signaling pathway, and CD8 T-cell receptor (TCR) pathway, etc. This result implied that immune response plays an important role in BPD development.

Figure 2 Gene set enrichment analysis (GSEA). Potential pathways, such as B cell receptor signaling pathway (A), cytokine receptor interaction (B), chemokine signaling pathway (C), the immune network for IGA production (D), primary immunodeficiency (E), senescence and autophagy in cancer (F), T cell receptor signaling pathway (G), CD8 TCR pathway (H) involved in BPD.

Composition of Immune Cells in BPD

As shown in Figure 3A, we used CIBERSORT to quantify the proportion of 22 immune cell types, and results indicated the proportion of B cells naïve (p < 0.05), NK cells resting (p < 0.01), macrophages M0 (p < 0.001), and neutrophils (p < 0.001) were significantly higher in the BPD group than those in the normal group. However, the proportion of T cells CD8 (p < 0.001), T cells CD4 naïve (p < 0.001), T cells CD4 memory resting (p < 0.01), and monocytes (p < 0.01) were significantly lower in the BPD group than those in the normal group. The MCPcounter was used to analyze the proportion of 10 immune cell types (Figure 3B), and the results revealed that the proportion of T cells (p < 0.001), CD8 T cells (p < 0.001), and myeloid dendritic cells (p < 0.01) were significantly lower in the BPD group than in the normal group, while neutrophils (p < 0.001) and endothelial cells (p < 0.05) were significantly higher in the BPD group than in the normal group. Furthermore, the ESTIMATE scores, Immune scores, and Stromal scores generated by the ESTIMATE software showed that Immune scores (p < 0.05) and Stromal scores (p < 0.05) were significantly different between the BPD group and the normal group (Figure 3C). These results further strengthened our hypothesis that immunocyte infiltration was highly involved in the BPD etiology.

Figure 3 Landscape of immunocytes infiltrates in BPD. (A) CIBERSORT scores of 22 immune cell types between BPD and normal groups. (B) MCPCounter scores of stromal cells or immune cells between BPD and normal groups. (C) Stromalscore, Immunescore, and ESTIMATEScore between BPD and normal groups. *p < 0.05, **p < 0.01, ***p < 0.001.

Identification of DEGs and WGCNA

We used the “limma” package to identify the DEGs between the BPD and normal groups. As shown in Figure 4A and Table S2, a total of 2000 DEGs were obtained based on the cutoff criteria with p.adj < 0.05 and∣logFC∣ ≥ 0.3, of which 834 genes were down-regulated and 1166 genes were up-regulated. And the top 50 downregulated and upregulated genes were presented in a heat map (Figure 4B). Next, we performed WGCNA based on the gene expression profile of 2000 DEGs. As shown in Figures 5A and B, we have chosen β = 5 to ensure a scale-free network and produce a hierarchical clustering tree with different modules. Then, three gene modules were identified by constructing the co-expression matrix (Figure 5C). As shown in Figure 5D, the brown module (p = 6.6e-7, r = −0.54), blue module (p = 1.1e-6, r = 0.53), and turquoise module (p = 1.1e-7, r = 0.57) showed significant correlation with BPD. Furthermore, we performed the Pearson’s correlation analysis between MM and GS in the blue module (Figure S1A, p = 3.6e-22, r = 0.47), brown module (Figure S1B, p = 7.6e-7, r = 0.45), and turquoise module (Figure S1C, p = 3.1e-34, r = 0.51).

Figure 4 Identification of differentially expressed genes (DEGs) in BPD. Volcano plot (A) and heatmap (B) of DEGs between BPD and normal groups.

Figure 5 Construction of co-expression modules. (A) The y-axis represents the scale-free topology model fit, the x-axis represents the soft threshold. The red star indicates that β (5, 0.92) was chosen to construct co-expression modules. (B) The y-axis represents the mean connectivity, the x-axis represents the soft threshold. The red star indicates that β (5, 21.11) was chosen to construct co-expression modules. (C) Dendrogram of all genes clustered based on a dissimilarity measure. The different colours represent the different co-expression modules. (D) Heatmap of the module–trait relationships. The brown module (p = 6.6e-7, r = −0.54), blue module (p = 1.1e-6, r = 0.53), and turquoise module (p = 1.1e-7, r = 0.57) showed significant correlation with BPD.

Functional Enrichment Analysis of Key Modules

To screen out immune-related modules, functional enrichment analysis was performed for each key module. Our results indicated that blue module genes were mainly enriched in cytoplasmic translation, ribonucleoprotein complex biogenesis, ribosome biogenesis, ribosome, EGFR tyrosine kinase inhibitor resistance, etc. (Figure 6A and Table S3). The brown module genes were mainly involved in cytoplasmic translation, aerobic respiration, oxidative phosphorylation, amyotrophic lateral sclerosis, Parkinson disease, etc. (Figure 6B and Table S4). The turquoise module genes were mainly enriched in lymphocyte differentiation, T cell differentiation, immune response-activating signal transduction, Th17 cell differentiation, T cell receptor signaling pathway, primary immunodeficiency, etc. (Figure 6C and Table S5). These findings revealed that immune-related pathways were mainly involved in the turquoise module. This module was, therefore, chosen for the subsequent analysis.

Figure 6 Enrichment analysis of blue (A), brown (B), and turquoise (C) module.

Protein–Protein Interaction Network and Validation of Core DEIRGs

We intersected IRGs and the genes in the three immune-related module (turquoise), and 31 DEIRGs were identified (Figure 7A). Then, the PPI network of DEIRGs was constructed, and core DEIRGs (CD8A) with the highest degree were identified (Figure 7B). We found the expression of CD8A was significantly down-regulated in the BPD group compared to the normal group (Figure 7C). Additionally, the expression of CD8A was further verified in multiple datasets and in-house experiment. As presentation in Figures 7D–G, the CD8A was significantly down-regulated in BPD group in GSE32472 (p = 1.3e-06), GSE25286 (p = 4.4e-04), GSE25293 (p = 4.7e-03), GSE99633 (p = 0.03), and in-house experiment (p = 2.4e-11), these results of which were consistent with the bioinformatics analysis results in GSE32472 dataset.

Figure 7 Protein–protein interaction network and core DEIRGs identification. (A) The intersection genes of turquoise module and IRGs. (B) Protein–protein interaction network of 31 DEIRGs. The darker the colour of the node, the more important it is in the network. Therefore, we chose the CD8A as the core gene. (C) CD8A expression level in GSE32472 dataset. (D–F) Validation of CD8A by external datasets (GSE25286, GSE25293 and GSE99633). (G) Validation of CD8A by in-house experiment.

In addition, the diagnostic ability of CD8A in distinguishing BPD from the normal samples exhibited a good diagnostic value, with an AUC of 0.81 in GSE32472, AUC of 1 in GSE25286, AUC of 1 in GSE25293, AUC of 0.8 in GSE99633, and AUC of 1 in our in vivo experiment (Figure S2).

GSEA and Gene Set Variation Analysis (GSVA) Analyses Between CD8A-Related Subgroups

We divided the GSE32472 dataset into two subgroups based on the median expression level of CD8A, and the top 50 DEGs associated with CD8A expression level were presented in Figure S3. GSEA and GSVA analyses were performed to further investigate the potential functions regulated by CD8A in BPD. Based on GSEA results, regulation of exocytosis (ES = −0.4654, P = 0.002), regulation of myeloid leukocyte mediated immunity (ES = −0.5991, P = 0.006), interleukin 6 mediated signaling pathway (ES = −0.5901, P < 0.001), leukotriene metabolic process (ES = −0.6913, P < 0.001), leukocyte mediated immunity (ES = −0.5346, P = 0.0078), toll-like receptor 4 signaling pathway (ES = −0.6471, P = 0.0081), interleukin 6 production (ES = −0.5325, P = 0.0059), and exocytosis (ES = −0.5979, P < 0.001) were enriched in the CD8A-low subgroup (Figure 8A). Furthermore, the GSVA results showed that immune- and inflammation-related pathways were activated in CD8A-low subgroup, including the immune effector process (p < 0.05), T cell antigen processing and presentation (p < 0.01), negative regulation of myeloid cell differentiation (p < 0.001), regulation of macrophage-derived form cell differentiation (p < 0.01), myeloid cell homeostasis (p < 0.05), interleukin 6 mediated signaling pathway (p < 0.001), positive regulation of myeloid leukocyte mediated immunity (p < 0.05), etc. (Figure 8B). These findings suggested that CD8A may show significant impact on the progression of BPD via regulating immune functions.

Figure 8 GSEA and GSVA analyses of CD8A in the BPD. (A) Single gene GSEA of CD8A. The gene sets with p-value < 0.05 were considered significantly enriched. (B) The histogram represented the results of GSVA between the two subgroups. *p < 0.05, **p < 0.01, and ***p < 0.001.

Immune Infiltration Analysis of CD8A

Subsequently, CIBERSORT was used to compare the proportion of 22 immune cell types, and we found that B cells memory (p < 0.01), T cells CD8 (p < 0.001), T cells CD4 naïve (p < 0.001), and T cells CD4 memory resting (p < 0.01) were decreased, whereas macrophages M0 (p < 0.001) and neutrophils (p < 0.001) were increased in the CD8A-low subgroup compared to the CD8A-high subgroup (Figure 9A). MCPcounter was applied to assess the proportion of 10 immune cell types, and results showed that T cells (p < 0.001), CD8 T cells (p < 0.001), cytotoxic lymphocytes (p < 0.001), B lineage (p < 0.05), NK cells (p < 0.01), and myeloid dendritic cells (p < 0.05) were decreased, whereas neutrophils (p < 0.001) and endothelial cells (p < 0.01) were increased in the CD8A-low subgroup compared to the CD8A-high subgroup (Figure 9B). ESTIMATE analysis showed that the Stromal score (p < 0.05) was increased, whereas the Immune score (p < 0.001) was decreased in the CD8A-low subgroup compared to the CD8A-high subgroup (Figure 9C).

Figure 9 Immune cell infiltration analysis. (A) CIBERSORT scores of 22 immune cell types between CD8A low and high groups. (B) MCPCounter scores of stromal cells or immune cells between CD8A low and high groups. (C) Stromalscore, Immunescore, and ESTIMATEScore between CD8A low and high groups. *p < 0.05, **p < 0.01, and ***p < 0.001.

Furthermore, we carried out a Spearman correlation analysis of the CD8A expression in immune cells. CD8A was positively correlated with T cells CD8 (p = 0, r = 0.84), T cells CD4 memory resting (p = 0, r = 0.52), T cells CD4 naïve (p = 0.0005, r = 0.46), B cells memory (p < 0.001, r = 0.38), eosinophils (p = 0.02, r = 0.26), T cells CD4 memory activated (p = 0.03, r = 0.24), T cells (p = 0, r = 0.83), cytotoxic lymphocytes (p = 0, r = 0.58), NK cells (p = 0.0003, r = 0.4), B lineage (p = 0.004, r = 0.32), and myeloid dendritic cells (p = 0.004, r = 0.32). CD8A was negatively correlated with macrophages M0 (p = 0, r = −0.58), and neutrophils (p =, r = −0.73) (Figure 10). Therefore, CD8A was associated with immune cell infiltrations in BPD.

Figure 10 Association of CD8A with immune infiltrates. Correlation analysis between immune cells and CD8A gene expression by using CIBERSORT (A) and MCPCounter (B) algorithms. Bold p values represent statistically significant differences.

Discussion

BPD has been recognized as one of the most common respiratory problems due to the prematurity in infants,²⁸ which remains to be a huge challenge in neonatal medicine. This longstanding complication has increased respiratory morbidity and also the incidence of neurological diseases, aggravating the financial burden of the family across the world.²⁹ Despite mounting effort has been devoted to improving the survival rates, the prevalence of BPD has remained the same or even increased.³⁰ To date, the pathogenesis of BPD is still not fully understood; novel and effective therapeutic strategies were still not available. Therefore, investigations about the underlying pathogenesis of BPD may provide the potential insights for early diagnosis and developing treatment strategies in the future. In the present study, we applied the bioinformatics approaches to investigate the differentially expressed immune-related genes and its immune infiltration in BPD. Three main findings were obtained from this study. Firstly, the results of GSEA and immune infiltration analysis indicated that immune-related pathways and immunocytes infiltrate were involved in the etiology of BPD. Second, a core biomarker (CD8A) for BPD was identified via a combination of DEG and WGCNA. Furthermore, CD8A expression was associated with immune cell infiltration in BPD, which may shed some light on the BPD immunotherapy.

CD8A is the gene encoding the CD8a chain of the dimeric CD8 protein, which is primarily responsible for cell-mediated immune defense and T-cell development.^31,32 CD8A has been widely reported as a potential prognosis and diagnostic marker for several diseases, such as inflammatory disorders and tumors.^33,34 The CD8A gene can also predict severity factors in chronic rhinosinusitis³³ and can act as a diagnostic biomarker for rheumatoid arthritis.¹⁷ The expression of CD8A was significantly up-regulated in the early N stage and related to favorable overall survival in lung adenocarcinoma.³⁴ A recent study has revealed that CD8A expression was down-regulated in a clinically relevant animal model of BPD.³⁵ Consistently, we also observed that CD8A expression was decreased in BPD in our study.

Epigenetic and genomic results reported the significant immune dysfunction in a mouse model of BPD.⁸ Increasing evidence also indicated that the imbalance between pro- and anti-inflammatory factors would induce injuries to the pulmonary interstitium in BPD patients.³⁶ Macrophages are important orchestrators of the pulmonary inflammatory and repair responses process, involving both the initiation and development of inflammation.³⁷ However, the specific pattern of immune cell infiltration in BPD is still unclear. As a result, we performed a comprehensive assessment of BPD immune cell infiltration to further investigate the potential pathogenesis related to the immune dysfunction. Our findings showed that T cells, CD8 T cells, and myeloid dendritic cells were significantly lower in the BPD group than in the normal group, while neutrophils and macrophages M0 were significantly higher in the BPD group than in the normal group. This finding revealed that immune infiltration is highly associated with the development of BPD. Furthermore, by the immunocytes infiltrate analysis of groups with low and high expression of CD8A, we also found T cells, CD8 T cells, NK cells, and myeloid dendritic cells levels were decreased, while neutrophils and macrophages M0 levels were increased in CD8A-low subgroup. It has been demonstrated that T lymphocytes plays a crucial role in the pathogenesis and progression of chronic lung disorder in babies born prematurely.³⁸ Emerging studies have suggested the CD8+ T cells may contribute to the immunopathology in BPD and may increase the risk of respiratory morbidity.³⁹ Furthermore, the high neutrophil-to-lymphocyte ratio is an early predictor of severe BPD and also a biomarker for the preterm infants with intrauterine infections.⁴⁰ The level of macrophage is significantly higher in the BPD patients.⁴¹ Macrophage activation is associated with IL-6/STAT3 axis in clinical and experimental BPD. Suppression of macrophage-associated IL-6 trans-signalling contributes to lung growth in experimental BPD.⁴² Dendritic cells are an integral component of lung development, and BPD pathology is significantly associated with dendritic cell recruitment in the lung tissue.⁴³ Neutrophils count may be a risk factor for the development of new BPD in neonates.⁴⁴ Thus, our findings are consistent with previous studies and reveal the importance of immune cell infiltration in the pathogenesis of BPD based on bioinformatics analysis. Considering all the results discussed above, we argued that CD8A is accountable for the development of BPD via the regulation of T cells, CD8 T cells, macrophages M0, neutrophils, and myeloid dendritic cells.

To our knowledge, this is the first study to characterize the immune infiltration landscape of BPD and investigate the role of CD8A in BPD; however, several limitations should be noted. First, since the validation set is not a human-related dataset and the sample size is limited, the expression of the core gene needs to be further validated in clinical samples with more participants engaged. Second, more experiments are necessary to further examine the interactions between immune cells and CD8A expression in BPD.

Conclusion

Our study not only provides a novel insight into the landscape of immune cells infiltrate associated with BPD etiology but also screened potential diagnostic marker for BPD. CD8A was down-regulated in hyperoxia-induced BPD, and it was associated with specific immune cell infiltrates in BPD. CD8A may involve in the development of BPD via regulation of immune functions. Our findings provided researchers and clinicians more evidence with immunotherapeutic strategies in the treatment of BPD.