Investigation of the miRNA-mRNA Regulatory Circuits and Immune Signatures Associated with Bronchopulmonary Dysplasia

Introduction

Bronchopulmonary dysplasia (BPD), characterized by impaired alveolarization and dysregulated vascularization in the lung, is an important cause of morbidity and mortality in preterm infants worldwide.¹ The younger the gestational age and the lower the birth weight, the higher the incidence and severity of BPD, with infants born at less than 1, 250 g accounting for 97% of BPD cases.² Advances in neonatal intensive care have improved the survival of preterm infants, but not the incidence of BPD. BPD survivors often have cardiopulmonary and neurological sequelae, resulting in heavy medical care and financial burdens on families and society.³ Unfortunately, the etiology and pathogenesis of BPD are still ambiguous so far.^4,5

In the past decades, immune reactions and inflammation have been topics of interest to investigators on the mechanisms of BPD.⁶ Besides hyperoxia exposure, antenatal and postnatal inflammation induced by injection of LPS could also cause BPD in animal models.^7,8 Immune cell-derived inflammatory cytokines, including IL-1, IL-6, TNF-α, were reported to be positively correlated with BPD severity.^9–11 Nevertheless, although excessive inflammatory responses were recognized to be key risk factors for BPD, which immune cell types played important roles over the course of BPD progression remained obscure. Numerous studies have reported that microRNAs (miRNAs) could affect BPD through multiple mechanisms, involving regulation of angiogenesis, extracellular matrix, and inflammation.^12–14 However, whether miRNA-mRNA regulatory circuits could regulate specific immune cells and thus affect BPD processes is not yet clear.

In this study, we compared miRNA and mRNA expression between BPD patients and controls through datasets obtained from the Gene Expression Omnibus (GEO). Then, miRNA-mRNA regulatory networks were constructed, and enrichment analysis of target genes was performed. Furthermore, we screened hub genes by protein-protein interaction (PPI) network analysis and inferred the proportion of 22 immune cell types in BPD using the CIBERSORT algorithm. Finally, we conducted real-time PCR (RT-PCR), immunohistochemistry (IHC), and flow cytometry to validate our results of bioinformatic analyses in BPD mice. Together, using bioinformatic techniques followed by confirmation in an experimental mouse model, this study aimed to investigate miRNA-mRNA regulatory circuits and immune signatures associated with BPD, which might improve current understanding of the mechanisms underlying BPD and provide clues to potential therapeutic targets for BPD.

Materials and Methods

Gene Expression Data Acquisition

The miRNA expression dataset (GSE108755) and mRNA expression dataset (GSE108754) were acquired from the GEO database (http://www.ncbi.nlm.nih.gov/geo). Each of the two datasets included blood samples from 5 BPD patients born at 28 weeks’ gestational age and 6 non-BPD age-matched controls. A flow diagram demonstrating the overall design of the study was shown in Figure 1.

Figure 1 Flow chart of the whole study.

Identification of Differentially Expressed miRNAs (DEMs)

The gene expression matrix from GSE108755 and its corresponding platform annotation file were downloaded from the GEO database. The matrix was log₂ transformed and quantile normalized using the “Normalize Between Arrays” function of “limma” package under R (version 4.2.1).¹⁵ Next, differential expression analysis between the BPD group and controls was performed using the limma algorithm. The DEMs were screened at a threshold of p-value of 0.05 and |log₂ fold change (logFC)| of 1.5. The “ggplot2” package was used to draw the heatmap.¹⁶

Construction of DEMs-Target Genes Interaction Network

The potential target genes of DEMs were predicted through the miRWalk, miRNet, miRDB, and TargetScan databases.^17–20 Concordant target genes of the four databases were identified by drawing a Venn diagram using an online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). The DEMs-target genes interaction network was visualized using Cytoscape software (version 3.8.2).²¹

Functional Enrichment Analysis

Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) analysis were performed using the DAVID functional annotation tool (version 6.8). GO terms include 3 categories: biological process (BP), cellular component (CC), and molecular function (MF). GO terms and KEGG pathways with a p-value < 0.05 and q value < 0.2 were considered to be significantly enriched. Visualization of the results of enrichment analysis was conducted using the “ggplot2” package in R.

PPI Network and Identification of Hub Genes

The PPI network of target genes was established using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (http://string-db.org/) with default parameters.²² Then, the interaction file containing the source and target nodes was downloaded and imported into Cytoscape software. CytoHubba, a plugin of Cytoscape, was used to rank the nodes. The top 30 nodes in the PPI network were selected as hub genes based on the maximal clique centrality (MCC) method, which was recognized to have a better performance on precise predictions of essential proteins.²³

Prediction of Upstream Transcription Factors of Hub Genes

The ChEA3 database, which contains 6 primary reference gene set libraries created from multiple resources, was used to predict the potential upstream transcription factors of hub genes.²⁴ In this study, the ENCODE library was selected.

Autophagy-Related Genes (ARGs)

To investigate the involvement of autophagy in the pathogenesis of BPD, a total of 222 genes related to autophagy were obtained from the Human Autophagy Database (http://www.autophagy.lu/index.html).

Validation of Hub Genes and ARGs

GSE108754, available in the GEO database, was used to validate the expression of hub genes. Differential gene expression analysis was performed using the “limma” package. Genes with p-value < 0.05 and |logFC| >1 were considered as differentially expressed genes (DEGs). Valid hub genes need to satisfy both of the following criteria:1) the direction of hub genes should be opposite to their corresponding DEMs; 2) DEGs. The intersections of DEGs and ARGs were also selected.

Immune Infiltration Analysis

Immune infiltration analysis of 22 types of immune cells was conducted using the CIBERSORT algorithm in R, and the results were visualized using the “ggplot2” package.²⁵

BPD Model and Tissue Collection

Neonatal C57BL/6J mice in this study were purchased from Guangzhou Zhiyuan Biotechnology. Mice were randomly assigned to the room air (RA) group (exposure to room air, n = 15) and hyperoxia group (exposure to 85% oxygen, n = 16) within the first 24 hours after birth. Dams were exchanged between hyperoxia and RA groups every 24 hours. Lung tissues of mice on day 7 after birth were collected for flow cytometry (n = 11). Other mice were sacrificed 14 days after birth for experiments except flow cytometry (n = 20).

RNA Extraction and RT-PCR

Total RNA was isolated from mice lungs with the RNA-easy Isolation Reagent (Vazyme, China). The total RNA was reverse transcribed using HiScript III RT SuperMix for qPCR (+ gDNA wiper) (Vazyme, China) analysis, ChamQ Universal SYBR qPCR Master Mix (Vazyme, China). RT-PCR was performed using QuantStudio 6 Flex real-time PCR system (Applied Biosystems, USA). Primer sequences are displayed in Table 1.

Table 1 Primers Used for RT-PCR of Mouse Genes

IHC

The lung paraffin sections were dewaxed to water, and then antigen retrieval was performed with pH 9.0 Antigen Retrieval Solution. After blocking endogenous peroxidase with 3% H₂O₂, the sections were blocked with 5% fetal bovine serum and incubated with the primary rabbit monoclonal anti-amyloid precursor protein antibody (ab32136, Abcam) at 4°C overnight. Sections were washed with PBS, then incubated with HRP-conjugated goat anti-rabbit IgG (ab6721, Abcam) for 1 hour at room temperature. After washing sections, the DAB substrate was added to the sections until color developed. Finally, the sections were examined using a Leica microscope and photographed.

Flow Cytometry

On day 7, the lungs were isolated from mice, grounded, and subsequently passed through a 70μm cell strainer. Ammoniumchloride-potassium lysing buffer was used for erythrocyte lysis. Multiparameter assessments were performed using BD FACSAria II (BD Biosciences) and data was analyzed with Flowjo software (version 10.6). After staining Zombie Aqua™ (Biolegend) to exclude dead cells, single-cell suspensions were incubated with the fluorochrome-conjugated antibodies in PBS containing 2% fetal bovine serum and then washed twice before detection. Antibodies specific to mouse used included CD16/32 (BD, 553,141), CD45-APC-Cy7 (Biolegend, 103,116), CD11b-AF700 (Biolegend, 101,222), Ly6G-PE (Biolegend, 127,608), CD62L-BV650 (Biolegend, 104,453), CD3-BV711 (Biolegend, 100,241), CD19-BV711 (Biolegend, 115,555), NK1.1-Pacific Blue (Biolegend, 109,816).

Statistical Analysis

Statistical analysis was performed using R (version 4.2.1) or GraphPad Prism (version 9). Data were expressed as mean ± SD. Comparisons of means between two groups were done with Student’s t-test. For correlation analysis, Spearman’s method was used. p-value < 0.05 was considered statistically significant.

Results

Identification of DEMs

The GSE108755 dataset was downloaded from the GEO database. After log₂ transformation and normalization, differential expression analysis was performed using the “limma” package in R. The results were shown as a heatmap in Figure 2A. Two miRNAs (miR-15b-5p and miR-20a-5p) were finally identified to be DEMs.

Figure 2 Identification of DEMs and their target genes. (A) Heatmap of miRNAs in GSE108755. (B) Venn diagram of potential target genes of DEMs predicted by miRNet, miRWalk, miRDB and TargetScan. (C) DEMs-target genes network. Abbreviation: BPD, bronchopulmonary dysplasia; miRNAs, microRNAs; DEMs, differentially expressed microRNAs.

Construction of DEMs-Target Genes Network

The potential target genes of DEMs were predicted using the miRWalk, miRNet, miRDB, and TargetScan databases, and 157 target genes were eventually screened (Figure 2B). The DEMs-target genes network was constructed using Cytoscape software (Figure 2C).

Functional Enrichment Analysis of Target Genes

To gain further insight into functions of the 157 target genes of DEMs, GO terms enrichment and KEGG pathway analysis were performed. BP analysis revealed that DEMs target genes were mainly enriched in myeloid cell differentiation, regulation of cell cycle phase transition, and macroautophagy (Figure 3A). CC analysis demonstrated that DEMs target genes were enriched in RNA polymerase II transcription regulator complex (Figure 3A). As for MF, DEMs target genes were enriched in protein serine/threonine kinase activity, protein-macromolecule adaptor activity, and SH3 domain binding (Figure 3A). KEGG pathway analysis indicated that target genes were involved in cellular senescence, cell cycle, MAPK signaling pathway, TGF-beta signaling pathway, VEGF signaling pathway, etc (Figure 3B).

Figure 3 GO and KEGG pathway analysis of the target genes of DEMs. (A) Enriched GO terms. (B) KEGG pathway analysis. Terms with p.adjust < 0.05 and q value < 0.2 were identified to be significant. The x-axis showed the gene count (A) and gene ratio (B) of each GO or KEGG term, and the y-axis showed names of enriched terms. Dot color and size represented p.adjust and gene counts, respectively. Abbreviations: DEMs, differentially expressed microRNAs; GO, gene ontology; KEGG, kyoto encyclopedia of genes and genomes; BP, biological processes; CC, cell component; MF, molecular function; p.adjust, adjusted p-value.

Construction of PPI Network and Screening of Hub Genes

A PPI network was established with the 157 target genes using the STRING database. Then, the interaction file was downloaded and imported into Cytoscape and processed with the cytoHubba plugin to identify hub genes. The top 30 hub genes are shown in Figure 4A. For hsa-miR-20a-5p, the hub genes were WEE1, E2F2, E2F3, MAPK1, CREB1, RBL2, AGO1, VEGFA, TNRC6B, USP28, APP, E2F5, CRK, SMAD5, YOD1, BCL2, OTUD4, CLOCK, REST and RGMB. For hsa-miR-15b-5p, the hub genes were CCND3, WEE1, E2F3, CHEK1, MAPK1, VEGFA, TNRC6B, BTRC, AKT3, APP, SMURF1, CRK, CLSPN, ACVR2A, PSMD7, BCL2, PHC3 and CDC27.

Figure 4 Identification of hub genes and potential upstream transcription factors. (A) PPI network of the top 30 hub genes for DEMs. The redder the node, the larger the MCC score. (B) Predicted upstream transcription factors of hub genes using ChEA3. Abbreviations: DEMs, differentially expressed microRNAs; PPI, protein-protein interaction; MCC, maximal clique centrality.

Prediction of Upstream Transcription Factors of Hub Genes

The upstream transcription factors of the top 30 hub genes, including ZNF384, BHLHE40, E2F4, GATA2, CTCFL, and CEBPB, were predicted using ChEA3 (Figure 4B).

Validation of Hub Genes and ARGs

The expression of hub genes was further verified using the GSE108754 dataset. However, only amyloid beta precursor protein (APP) was significantly downregulated in the BPD group compared with controls (Figure 5). The expression of CHEK1, BCL2, CCND3, E2F3, CLSPN, PSMD7, VEGFA, SMAD5, CLOCK, and AKT3 were consistently decreased but did not reach statistical significance (p-value > 0.05 or |logFC| < 1). ARGs expression levels were also assessed. Using the cutoff value of p-value < 0.05 and |logFC| > 1, two ARGs were found to be upregulated (DLC1 and PARP1) and two downregulated (NLRC4 and NRG1) (Table 2).

Table 2 Differentially Expressed ARGs

Figure 5 Validation of hub genes in the GSE108754 dataset. Student’s t-test. *p < 0.05. ns not significant.

Immune Cell Infiltration Analysis

The relative proportions of 22 types of immune cells in BPD and normal samples were inferred using the CIBERSORT algorithm in R. Figure 6A showed the composition of immune cells in each sample of the GSE108754 dataset. The correlation heatmap of immune cells demonstrated that there was a strong positive correlation between resting dendritic cells and Tregs (r > 0.8). Resting NK cells and activated CD4⁺ memory T cells also had a strong positive correlation, while monocytes and CD8⁺ T cells were strongly negatively correlated (r < −0.8). Neutrophils were positively correlated with resting CD4⁺ memory T cells and negatively correlated with naïve B cells (Figure 6B). Box plots showed that naïve B cells were significantly increased in the BPD group, while neutrophils were significantly reduced (Figure 6C). Interestingly, the expression of APP was found to be negatively correlated with neutrophils (Figure 6B).

Figure 6 Relative proportions of 22 immune cell types in BPD and normal groups. (A) Compositions of immune cells in each sample. (B) Correlation heatmap of immune cells. Spearman’s method. *p < 0.05, **p < 0.01. (C) Differences in immune cell compositions between BPD and normal groups. Student’s t-test. **p < 0.01. Abbreviations: APP, amyloid beta precursor protein.

Validation of Hub Genes in Mouse Models of BPD

Consistent with the phenotypes of human BPD, hyperoxia-induced BPD mice exhibited disruption of alveolar growth and enlargement of alveolar spaces (Figure 7A), with decreased radial alveolar count (RAC) and increased mean linear intercept (MLI) (Figure 7B). The expression of genes previously identified was verified using RT-PCR in BPD mice. The expression level of App was significantly reduced in the lungs of BPD mice compared with the RA group. However, contrary to our expectation, Nrg1 expression was elevated (p > 0.05). Parp1, Nlrc4, and Dlc1 showed trends consistent with previous bioinformatic analyses, but the differences did not reach statistical significance (p > 0.05) (Figure 7C). IHC staining further confirmed the downregulation of APP levels in BPD mice (Figure 7D).

Figure 7 Experimental verification of bioinformatics findings in BPD mice. (A) The morphology of lung tissues in the RA and hyperoxia groups under the microscope at magnifications of 100 x and 200 x. (B) Quantitation of lung histopathological changes via RAC and MLI. Student’s t-test. **p < 0.01. (C) RT-PCR analysis of App, Nlrc4, Parp1 and Dlc1. Student’s t-test. *p < 0.05. ns, not significant. (D) Representative IHC staining images of APP under the microscope at magnifications of 200 x. The black arrows point to APP-positive areas. (E) Proportion of total neutrophils in the RA and BPD groups. Student’s t-test. ns, not significant. (F) Proportions of CD62L+ and CD62L− neutrophil subsets in the RA and BPD groups. Student’s t-test. **p < 0.01. (G) Mean fluorescence intensity of CD62L on neutrophils in the RA and BPD groups. Student’s t-test. **p < 0.01. Abbreviations: RA, room air; APP, amyloid beta precursor protein; RAC, radial alveolar count; MLI, mean linear intercept; MFI, mean fluorescence intensity.

CD62L⁺ Neutrophils Significantly Decreased in BPD Mice

Through flow cytometry, we found that the proportion of total neutrophils did not differ between the RA and BPD groups (Figure 7E). However, the population of CD62L⁺ neutrophil subset (p < 0.01) was dramatically decreased, whereas the population of CD62L⁻ neutrophil subset (p < 0.01) was significantly increased in BPD mice compared with controls (Figure 7F and G). These findings suggest that in response to hyperoxia, the activation status in neutrophils altered in BPD mice, suggesting a potential role of neutrophils in hyperoxia-induced injury.

Discussion

Despite rapid advances in medical technology, the incidence of BPD remains high, and effective therapies are lacking, largely due to limited understanding of its pathogenesis. The mechanisms underlying BPD are complex, involving ferroptosis, IL-6, Wnt signaling, etc.^9,26,27

This study aimed to identify the miRNA-mRNA regulatory circuits that are most likely related to BPD and explore the roles of immune cells in BPD pathogenesis. Two miRNAs (miR-15b-5p and miR-20a-5p), both upregulated, were identified as DEMs. Among them, miR-15b-5p was reported to regulate the proliferation and apoptosis of human vascular smooth muscle cells.²⁸ miR-20a-5p was found to inhibit autophagy in ovarian cancer.²⁹ Notably, the role that these two miRNAs play in BPD is still unclear.

GO enrichment analysis showed that DEMs target genes were mainly involved in myeloid cell differentiation, regulation of cell cycle phase transition, macroautophagy, ubiquitin ligase-substrate adaptor activity, and protein serine/threonine kinase activity. KEGG analysis showed that target genes were mainly enriched in the cell cycle, MAPK signaling pathway, TGF-β signaling pathway, and VEGF signaling pathway. Through PPI network analysis, 30 hub genes were screened. However, among these hub genes, only APP was significantly decreased in the mRNA expression dataset (p-value < 0.05 and |logFC| > 1). Significant downregulation of APP was also confirmed by RT-PCR and IHC in the mouse model of BPD. APP is a highly pleiotropic protein involved in numerous cellular functions.³⁰ APP misprocessing and deposition of pathogenic fragments could cause neuron and synapse loss, which is the most recognized hypothesis in the pathogenesis of Alzheimer’s disease.³¹ However, according to the Bio Gene Portal System (BioGPS, http://biogps.org/) database, apart from the brain, APP is widely expressed in many tissues, including the bronchus and lung. APP was reported to regulate global protein synthesis in a variety of human dividing cells.³² Depletion of APP could cause cell size abnormalities and death through G0 arrest.³³ To date, the exact function of APP in the lung has not been reported. Previous studies reported that APP played an important protective function in cerebral vascular endothelium by maintaining endothelial nitric oxide synthase (eNOS) expression and mediating reactivity to vascular endothelial growth factor (VEGF).^34,35 Aberrant pulmonary alveolarization and dysregulated vascularization are key features of BPD, and impaired VEGF and eNOS signaling have been proved to be related to BPD progression.^36,37 Herein, we speculated that miR-15b-5p/miR-20a-5p-APP might promote BPD via impairing cell cycle and VEGF/eNOS signaling. Dysregulation of autophagy has been reported to be associated with various lung diseases.^38,39 However, the role that autophagy plays in BPD pathogenesis is not clear. We identified 4 different ARGs (DLC1, PARP1, NLRC4 and NRG1) in BPD and verified them via RT-PCR, but none of the differences reached statistical significance.

Exaggerated immune activation and inflammation are key risk factors for BPD.⁴⁰ Continued inflammation recruits immune cells, which in turn secrete more pro-inflammatory cytokines, ultimately aggravating the inflammation and causing lung injury. Neutrophils are the earliest immune cells to infiltrate the damaged tissue after injury and have been proven to be associated with chronic lung diseases.⁴¹ However, little research has been conducted regarding changes of circulating or lung neutrophils in BPD. In this study, we estimated the composition of immune cells in BPD patients and controls using the CIBERSORT algorithm. We found that neutrophils were reduced in the peripheral blood of BPD patients, which was contrary to some previous reports.^42–44 Reasons for this inconsistency are not clear, partly due to different gestational ages and postnatal days as the number of neutrophils in the peripheral blood change drastically within the first week after birth. Besides, it was reported that neutrophils remained in the pulmonary vascular bed during acute lung injury, leading to a reduction of peripheral blood neutrophil counts.⁴⁵ Interestingly, APP levels were strongly negatively correlated with the proportion of neutrophils.

In the present study, we found that lung neutrophils in newborn mice could be divided into two subgroups: CD62L⁺ neutrophils and CD62L⁻ neutrophils, and the proportion of CD62L⁻ neutrophil subset was dramatically increased in the lungs of BPD mice. As CD62L shedding is a sign of neutrophil activation, the decreased expression of surface CD62L suggests that pulmonary neutrophils are more activated in BPD mice.^46,47 It was reported that alveolar neutrophils with reduced CD62L expression display delayed apoptosis and enhanced respiratory burst in acute respiratory distress syndrome.⁴⁸ The Wnt pathway is critical both during embryonic lung development and in BPD.²⁶ Recently, it was reported that neutrophil extracellular traps (NETs) could promote BPD progression via the Wnt/β-catenin pathway.⁴⁹ As NETs are released by activated neutrophils, it is possible that CD62L⁻ neutrophil subset promotes BPD progression via NETs. Further experiments are required to test this hypothesis.

In all, based on bioinformatic analyses, we constructed miRNA-mRNA regulatory networks, which enabled a more comprehensive understanding of genes and interactions related to BPD. In addition, we estimated changes of immune cell proportions in BPD patients applying CIBERSORT. Finally, we performed RT-PCR, IHC and flow cytometry in BPD mice to confirm our findings. However, our study also has some limitations. First, our study relied on public datasets with insufficient sample size and incomplete clinical information, which might bias the results and limit the generalizability of the study findings. Second, the specific role and mechanism of APP in BPD pathogenesis remains unclear, further experiments, particularly single-cell RNA sequencing, are warranted. Besides, the function of CD62L⁻ neutrophil subset and its role in BPD pathogenesis deserve particular attention in future research.

Conclusion

The results of this study indicated that the downregulation of APP mediated by miR-15b-5p and miR-20a-5p might be associated with BPD progression. Besides, increased CD62L⁻ neutrophil subset may promote BPD progression via excessive inflammatory responses. Our study might help to understand the pathogenesis of BPD and provide potential therapeutic targets for BPD.