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

Abstract

Background: Bronchopulmonary dysplasia (BPD) refers to a chronic lung disease which is commonly observed in preterm infants. It can usually be caused by several pathological processes that endanger the long-term lung development, such as inflammation and immune dysfunction.

Methods: In this study, a bioinformatics approach was applied to identify the differentially expressed immune-related genes (DEIRGs). We downloaded the transcriptional profiles (GSE32472 dataset) from the Gene Expression Omnibus (GEO) database and performed gene set enrichment analysis (GSEA). Cell type Identification By Estimating Relative Subsets of RNA Transcripts (CIBERSORT), microenvironment cell populations counter (MCPcounter), and Estimation of STromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) were used for the analysis of the immune cell infiltration landscape of BPD. A weighted co-expression network was subsequently constructed using weighted gene co-expression network analysis (WGCNA) to screen candidate differentially expressed immune related genes (DEIRGs).

Results: GSEA results indicated that immune-related pathways were mainly involved in BPD. Ten significantly different immune cell types were observed between BPD and normal groups. A total of 228 DEGs in the turquoise module were identified, and 31 DEIRGs were further identified. Cluster of the differentiation 8 alpha (CD8A) expression was down-regulated in BPD, and its expression was validated by the GSE25286, GSE25293, GSE99633 datasets and qRT-PCR. In addition, CD8A expression was closely associated with immune cells infiltration, especially T cells CD8 and neutrophil.

Conclusion: A distinct immune cell infiltration landscape was found between BPD and normal group. CD8A can be a novel candidate biomarker for BPD, which plays an essential role in the onset and progress of hyperoxia-related BPD via the disruption of immune cell functions.

Keywords: CD8A; bronchopulmonary dysplasia; hyperoxia; immunocytes infiltrate; preterm infants.

Figure 1

The schematic graph of the bioinformatics analysis.

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.

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.

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.

Figure 6

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

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.

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.

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.

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.