Introduction

BPD, a potentially fatal illness, is a significant public health concern because of its association with an aberrant immune response to infections and severe dyspnea.¹ Despite major efforts to minimize injurious but often life-saving postnatal interventions (such as oxygen, mechanical ventilation and corticosteroids), BPD remains the most common frequent complication of extreme preterm birth.^2,3 Nearly one-third to one-half of all infants born before 28th week of gestation develop BPD.⁴ However, the precise mechanisms of BPD remain unclear.

An increasing body of researches suggests that inflammatory mediators, including cytokines and growth factors, play a role in BPD pathogenesis, and their activity can be upregulated by inflammatory processes.^5,6 Autophagy has been implicated in the regulation of immune cell infiltration and is associated with lung injury.⁷ According to an investigation, BPD is positively correlated with the occurrence of intrauterine infection.⁸ Besides, the umbilical cord tissue not only reflects but also can influence the intrauterine environment. Moreover, umbilical cord blood contains not only hematopoietic stem cells but also NK cells, lymphocytes, monocytes, mesenchymal stem cells, and other kinds of cells.⁹ The wide origins of these many types of cells make single-cell sequencing possible, allowing us to identify more extensively the links between genes and immune cells, as well as understand the probable mechanisms of action of BPD-related genes and immune cells. Therefore, we could evaluate the intrauterine infection by detecting the umbilical cord tissue to identify possible factors in the pathogenesis of BPD.

Microarray expression profiling provides a comprehensive portrait of the transcriptional world, enabling us to view the organism as a “system” that is more than the sum of its parts.¹⁰ Machine-learning algorithms have shown great promise in assessing transcriptome data and locating genes that are of biological significance.¹¹ Machine learning and bioinformatics is an innovative research method. This approach helps improve the precision of our research on the relationship between gene targets and diseases, and enables us to more efficiently discover additional gene targets that we have not previously identified. However, there are few investigations of biomarkers based on machine learning about BPD. What remains unclear is the extent to which biologic markers obtained from the umbilical cord might predict the occurrence of BPD and have an influence on preterm infants. This study provides new insights into the potential therapeutic interventions of BPD.

Method

Source of Microarray Data

In this study, we downloaded a microarray dataset (GSE8586) from the Gene Expression Omnibus (GEO) database¹² (https://www.ncbi.nlm.nih.gov/geo/). The GSE8586 dataset was applied to identify DEGs, which contained umbilical cord tissue samples from infants born before 28th week of gestation, including 20 infants with BPD and 34 neonates with unaffected controls, respectively. BPD is defined as the requirement for supplemental oxygen at 36 weeks postmenstrual age. To ensure the reliability and accuracy of our subsequent analyses, we used the Affy package of R to perform normalization and background correction. The “sva” package was used to perform batch normalization. Besides, we downloaded 4409 autophagy-related genes (ARGs) in GeneCards. A schematic diagram of our study is depicted in Figure 1.

Figure 1 The overall research design. Screening DE-ARGs through GSE8586 and autophagy-related genes firstly. Using cell-type identification with CIBERSORT to analyze the immune status of BPD. The hub genes were identified by the use of Lasso and Cytoscape with four machine-learning algorithms (MCC, Degree, MCODE and LASSO). To validate hub genes, we use ROC, single cell sequence and IHC in hyperoxia rats. At last, we searched the drug targets of these hub genes, and established a nomogram model for predicting the risk of BPD.

Sample Detection and Differential Expression Genes Analysis

The DEGs between normal and BPD groups in the GSE8586 dataset were identified using the limma package. |log2 FC| > 0.5 and p value < 0.05 were regarded as significant. Besides, volcano plot was produced to highlight the differential expression of DEGs. A heatmap was generated based on the screened DEGs.

Protein–Protein Interaction and Module Analysis

Protein–protein interactions (PPIs) were analyzed using STRING for interacting Genes/Proteins dataset (http://cn.string-db.org,version10.0), and an interaction with a combined score >0.4 was considered statistically significant.¹³ We used the Cytoscape, an open-source bioinformatics software platform, to visualizing molecular interaction networks.¹⁴ Besides, the most significant modules in the PPI networks were identified using MCODE, MCC and Degree. The criteria for selection were as follows: MCODE scores > 5, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100 and κ-score=2.

Least Absolute Shrinkage and Selection Operator (LASSO)

To reduce data dimensionality and identify genetic biomarkers for BPD, we used LASSO with R (4.2.1) version and glmnet [4.1.7] to analyze the cleaned data and obtain variable coefficient values, logarithm of lambda values, likelihood values or classification error rates.¹⁵ This allows variable selection and complexity adjustment by controlling the lambda parameter. Finally, optimum gene biomarkers for BPD were discovered using biomarkers produced by Lasso algorithm and 3 distinct algorithms in Cytoscape.

Annotation and Enrichment Analysis

Functional enrichment was applied to the data to verify the likely functions of potential targets. Gene ontology (GO) contained biological pathways (BP), cellular components (CC) and molecular functions (MF). Kyoto Encyclopedia of Genomes (KEGG) pathway enrichment was used to investigate the activities of genes and the high-level genomic information to connect to those functions.¹⁶ Both GO and KEGG analyses were carried out by cluster Profiler [4.4.4], GOplot [1.0.2], ggplopt2 [3.3.6], igraph [1.4.1] and ggraph [2.1.0].

Immune Infiltration Analysis

CIBERSORTx, an analytical tool from the Alizadeh Lab and Newman Lab, was utilized to impute gene expression profiles and estimate the abundances of 22 different member immune cell types in a mixed cell population,¹⁷ using gene expression data in BPD samples from the GSE8586 dataset. All immune cell-type fractions were summed up to one for each sample. The Spearman association between unique diagnostic markers and immune invading cells was analyzed using the “ggstatsplot” and “ggplot2” packages to illustrate the results.

Predicting Drug Targets for Hub Genes

We accessed the DSigDB drug database on Enrichr (https://maayanlab.cloud/Enrichr/) and inputted the five hub genes we studied. Using P < 0.05 as a condition, we screened for drug targets to predict the drug targets of the genes.

Receiver Operating Characteristic (ROC) Curve

The ROC curve was analyzed using R (4.2.1) version, pROC [1.18.0] package, and the results were visualized using ggplot [3.3.6] package.¹⁸ The Area Under Curve (AUC) of the ROC curve is used to evaluate the accuracy of the diagnostic model for autophagy-related genes in BPD. The AUC value ranges from 0.5 to 1, with values closer to 1 indicating better diagnostic performance of the variable in predicting the outcome. We obtained the final AUC results with the function of pROC and evaluation of the AUC.

Single-Cell Sequencing Analysis

The single-cell sequencing data about the hub genes were used from Single Cell Portal. We used single-cell sequencing data from human lung tissues to explore the expression of hub genes.

Establishment of Nomogram

A diagnostic nomogram based on the selected candidate DE-ARGs was established using the “rms” package in the R software. The consistency between our predicted values and reality was assessed using a calibration curve. ROC curve analysis and decision curve analysis (DCA) were carried out to determine whether or not the decisions made based on the model were favorable to the patients.

Hematoxylin and Eosin (HE) Staining

After the left lung tissues were fixed with 4% paraformaldehyde, dehydrated with gradient ethanol, 75%, 85%, 90%, 95% alcohol and anhydrous ethanol I, II for 2 h, respectively. Xylene I, II, III for 40 min, respectively. 65°C molten paraffin I for 0.5 h, 65°C molten paraffin II for 1 h, 65°C molten paraffin III for 2 h and 45 min, and paraffin embedded, slides with 5μm paraffin sections were dewaxed and hydrated. Then, the sections incubated with hematoxylin for 5 min, stained with eosin for 2 min. Finally, using the microscope to observe the morphology of lung tissues (×200 magnification).

Immunohistochemical (IHC)

The methodology employed utilizes rats left lung tissue, initially fixed and embedded, subsequently sectioned and subjected to a process of fixing and antigen repair. The sections are then immunostained and sealed for observation. The expression of hub genes were evaluated according to the area of positive cells. The primary antibodies used for IHC staining overnight were BRIX1 (Cat No. 17295-1-AP; Proteintech; 1:100), JUN (Cat No. 24909-1-AP; Proteintech; 1:100), PES1 (Cat No. 13553-1-AP; Proteintech; 1:200), NR4A1 (Cat No. 25851-1-AP; Proteintech; 1:200), RRP9 (Cat No. 13553-1-AP; Proteintech; 1:100). This study was approved by the Institutional Animal Experiment Committee of Southern Medical University Nanfang Hospital.

Results

Screening of DEGs and Immune Landscape Analysis in BPD

After normalizing each sample in GSE8586 dataset (Figure 2A), the expression patterns of the DEGs were visualized in a clustering heatmap (Figure 2B). Then, we found 647 DEGs between BPD and control groups based on the criterion with |log2 FC| > 0.5 and p value < 0.05 (Figure 2C and D). In addition, we used CIBERSORT to estimate the immune infiltration between BPD and control group. The results showed that T cells follicular helper (Tfh cells), Monocytes, plasma cells, B cells memory, Eosinophils, NK cells resting, T cells CD8 and Macrophages M1 were the main immune cells that infiltrated the plaque (Figure 2E). The abundance of B cells memory, T cells gamma delta, Eosinophils (p value < 0.05), Macrophages M1 and Dendritic cells resting (p value < 0.01) were significantly increased in patients with BPD, while Tfh cells, NK cells resting and Monocytes (p value < 0.05) were significantly decreased in preterm infants with BPD (Figure 2F). The correlation between immune cells in BPD was further investigated by Pearson’s correlation analysis (Figure 2G).

Figure 2 Selection of DEGs and immune infiltration analysis from dataset GSE8586. (A) Normalization of each sample data. (B) Heatmap displaying the expression levels of various genes. (C) Selection of DEGs illustrated with a volcano plot. (D) Differential analysis chart presenting the DEGs. (E) Cibersort immune infiltration analysis of all samples containing DEGs from GSE8586. (F) Expression profiles of 22 immune cell types in BPD and control groups. (G) Heatmap of the correlation among immune cells. ****P<0.0001, ***P<0.001, **P<0.01, *P<0.05.

Functional Enrichment of 73 DE-ARGs

To further explore whether autophagy played a role in BPD, we downloaded 4409 autophagy-related genes (ARGs) from GeneCards. There were 73 the differentially expressed and autophagy-related genes (DE-ARGs) by overlapping the DEGs in GSE8586 and ARGs, including 43 upregulated genes and 30 downregulated genes (Figure 3A). To gain a deeper understanding of the biological functions of the DE-ARGs, we conducted functional analysis. The results of GO enrichment analysis revealed that these DE-ARGs were linked to response to mechanical stimulus, release of cytochrome c from mitochondria, apoptotic mitochondrial changes, connexin complex, gap junction, preribosome, R-SMAD binding, gap junction channel activity and wide pore channel activity. KEGG analysis was associated with IL-17 signaling pathway and colorectal cancer (Figure 3B–D).

Figure 3 Enrichment analysis of differentially expressed genes and autophagy-related genes. (A) Intersection of DEGs from GSE8586 with autophagy-related genes (ARGs) from the Genecards dataset. (B) Chord diagram illustrating combined GO and KEGG analysis with LogFC. (C) GO and KEGG analysis of the 73 intersecting genes. (D) Network diagram of the GO and KEGG enrichment analysis results.

Identification and Immune Landscape Analysis of 5 Hub Genes as Diagnostic Genes for BPD

Then, the protein–protein interaction network of 73 DE-ARGs was conducted, with 30 downregulated genes and 43 upregulated genes (Figure 4A). And the Cytoscape was used to visualize these more interacting genes (Figure 4B). To evaluate the diagnostic potential of autophagy-related genes between BPD and control samples, we conducted 3 distinct machine learning algorithms (MCC, Degree and MCODE) on the GSE8586 dataset (Figure 4C–G). By utilizing a Venn diagram to compare the overlapping regions of DEGs and key module genes, we were able to identify 10 genes as module genes (Figure 4H). To further minimize the potential diagnostic genes for BPD, we used the LASSO algorithm to identify TBL3, PES1, BRIX1, RRP9, JUN, DUSP1, NR4A1 (Figure 5A and B). Then, we drew a one-dimensional bar chart of 7 genes based on their lambda min, and selected 5 most influential genes, PES1, BRIX1, RRP9, JUN, and NR4A1 (Figure 5C). Finally, we used ROC to validate the 5 hub genes, each AUC of hub genes was higher than 0.67 (Figure 5D). Subsequent GO enrichment analysis showed that 5 hub genes were correlated with ribosome biogenesis, rRNA metabolic process, rRNA processing, transcription regulator complex, preribosome, preribosome large subunit precursor, MAP kinase tyrosine phosphatase activity, U3 snoRNA binding, protein tyrosine phosphatase activity, while KEGG enrichment analysis presented these genes were linked to MAPK signaling pathway, fluid shear stress and atherosclerosis and cocaine addiction (Figure 5E–G). In addition, PES1 was correlated with RRP9, while Jun was correlated with NR4A1 (Figure 5H). And the expression of 5 hub genes were significantly different between the BPD and control groups (Figure 5I). Therefore, BRIX1, JUN, PES1, NR4A1, and RRP9 may be potential therapeutic targets for patients with BPD.

Figure 4 PPI interaction and module analysis of autophagy-related genes. (A) PPI network analysis of autophagy-related differentially expressed genes (DE-ARGs). (B) Visualization of more interacting DE-ARGs using Cytoscape software, with 23 downregulated genes (blue) and 11 upregulated genes (red). (C and D) Selection of the top 10 genes using the MCC (Maximal Clique Centrality) and Degree algorithms. (E–G) Identification of 3 clusters using the MCODE (Molecular Complex Detection) algorithm. (H) Intersection of the three algorithms in a Venn diagram to identify 10 common genes.

Figure 5 Screening of hub genes and analysis of their functional enrichment and immune infiltration. (A) Diagnostic Lasso coefficient screening. (B) Diagnostic Lasso variable trajectory. (C) One-dimensional bar chart screening of the top 5 significant impact factors. (D) ROC curve confirming the reliability of the 5 hub genes. (E) Heatmap presentation of hub gene expression. (F) Cluster analysis of hub genes using GO and KEGG. (G) Combined LogFC enrichment analysis of GO and KEGG. (H) Correlation among hub genes. (I) Expression of five hub genes between BPD and control groups; a Mann–Whitney U-test (Wilcoxon rank sum test) was employed for statistical analysis, with bean plots illustrating the comparative expression of the five hub genes across multiple dimensions between the two groups. (J) Expression of the five hub genes across 22 types of immune cells. ***P<0.001, **P<0.01, *P<0.05.

To study immune cell infiltration in BPD and control, the corrected expression matrix was subjected to CIBERSORT to estimate the abundances of infiltrating immune cells in a mixed-cell population.¹⁷ We analyzed immune cell infiltration in the 5 hub genes. We found that JUN was positively correlated with plasma cell, T cells CD8, Tfh cells, NK cells resting, and monocytes, while negative correlated with B cells, Macrophages M2 (p value < 0.01), Macrophages M1, T cells CD4, and Dendritic cells resting (p value < 0.05). NR4A1 was positively correlated with T cells CD4, Tfh cells, monocytes (p value < 0.05). BRIX1 was positive correlated with plasma cells, while negatively correlated with Eosinophils (p value < 0.05). RRP9 was positive correlated with Dendritic cell, while negatively correlated with Neutrophils (p value < 0.05). However, PES1 was not correlated with these immune cells (Figure 5J).

Structure and Function of Hub Genes

PES1 (Pescadillo ribosomal biogenesis factor 1) located on chromosome 22 is a protein-coding gene, which is related to rRNA processing in the nucleus and cytosol, processing of capped intron-containing pre-mRNA and processing cell cycle regulation.^19,20 BRIX1, located in chromosome 5 and nucleolus, is involved in ribosomal large subunit assembly.²¹ RRP9 (Ribosomal RNA processing 9) encodes a member of the WD-repeat protein family, which participate in the processing and modification of pre-rRNA. Ectopic expression of RRP9 promotes tumor cell proliferation, colony formation, and cell migration.²² JUN is the putative transforming gene of avian sarcoma virus 17. It encodes a protein, which is highly similar to the viral protein, and which interacts directly with specific target DNA sequences to regulate gene expression. Together with FOSB, plays a role in activation-induced cell death of T cells by binding to the AP-1 promoter site of FASLG/CD95L, and inducing its transcription in response to activation of the TCR/CD3 signaling pathway.²³ NR4A1 (Nuclear Receptor Subfamily 4 Group A Member 1) encodes a member of the steroid-thyroid hormone-retinoid receptor superfamily. The encoded protein acts as a nuclear transcription factor.²⁴ Translocation of the protein from the nucleus to mitochondria induces apoptosis (Supplementary Table 1 and Supplementary Figure 1).

Drug Targets in Hub Genes

Using data from the DSigDB database, we identified Quinpirole HL60 DOWN, Pergolide HL60 DOWN, Atrazine CTD 00005450, Cylindrospermopsin CTD 00003136, Alprostadil HL60 DOWN, 5-Fluorouracil CTD 00005987, Genistein CTD 00007324, Progesterone CTD 00006624, Calcitriol CTD 00005558, Estradiol CTD 00005920, Hydrogen peroxide CTD 00006118, and Trichostatin A CTD 00000660 were identified as candidate drugs for most gene interactions (Table 1). Since these iconic drugs have been detected in the core genes, they may target PES1, BRIX1, RRP9, JUN and NR4A1.

Table 1 Drug Targets in Hub Genes

Clinical Implications

To further analyze the relevance of these drugs in treating BPD and how they interact with the identified hub genes, we conducted a detailed research on these hub genes-targeted drugs’ potential therapeutic effects in BPD. Clinically, Alprostadil, prostaglandin E1, can be used to dilate pulmonary blood vessels, alleviating the respiratory distress symptoms caused by BPD-related pulmonary vascular narrowing.²⁵ Estradiol may improve lung function and outcome in preterm infants.²⁶ Calcitriol is a receptor agonist of vitamin D. Vitamin D protects against hyperoxia-induced BPD in newborn rats by regulating the vitamin D-VDR signaling pathway.²⁷ Trichostatin interferes with the development of alveoli in newborn rats²⁸ (Table 2). These gene-targeting drugs provide practical and promising clinical value for the treatment of BPD diseases.

Table 2 Hub Genes-Targeted Drugs’ Potential Therapeutic Effects in BPD

Single-Cell Sequencing Analysis

We used the Single-Cell Portal to investigate single-cell sequencing of 5 hub genes expression in the lung, which suggested that the hub genes were highly expressed in the human lung. The most significant expression of hub genes is JUN. These data were classified for 13 subgroups: apoptotic cell, B cell, endothelial cell, epithelial cell, fibroblast, lymph node lymphatic vessel endothelial cell, lymphocyte, mast cell, mitotic cell cycle, myeloid cell, smooth muscle cell, typeI pneumocyte and typeII pneumocyte. Interestingly, all hub genes mainly located in B cell, epithelial cell, fibroblast, endothelial cell, smooth muscle cell, typeI pneumocyte and typeII pneumocyte (Figure 6).

Figure 6 Expression of hub genes in single-cell sequencing dataset. (A) Using single-cell sequencing data from human lung tissue via the Single Cell Portal, cells are categorized into 13 subgroups. (B–F) In the single-cell sequencing data, the expression of hub genes ranges from strong to weak: JUN, NR4A1, BRIX1, PES1, RRP9.

Establishment of Nomogram

We developed a diagnostic nomogram based on these five biomarkers (Figure 7A). ROC curves revealed that the model distinguished between BPD and unaffected control samples (Figure 7B). The predictive accuracy of the nomogram was highlighted using the calibration curves (Figure 7C). The line in the DCA curve remained above the gray and black lines from 0 to 1, suggesting that decisions based on the nomogram model may favor patients with BPD (Figure 7D).

Figure 7 Establishment of the nomogram. (A) Establishment of a nomogram based on the five biomarkers. (B) ROC curve for the model. (C) Calibration curve revealed the predictive ability of the nomogram. (D) The nomogram may be used to make decisions that favor BPD patients.

Validation of Hub Genes in Animal Experiment

In order to show the expression levels of ARGs in rats lung, animals’ studies were divided into two groups, one exposure of male and female newborn SD rats background to hyperoxia (85% O₂) for 14 days, while other group exposure to normoxia (21% O₂), as shown in Figure 8A. Alveolarization was greater in 85% O₂ rats than in 21% O₂ rats (Figure 8B). The IHC results showed that expression of PES1, BRX1, RRP9, JUN, NR4A1 was lower in the hyperoxia group compared to normoxia group (Figure 8C and D). There were 3 animals analyzed in each group, P < 0.05.

Figure 8 Validation of hub genes in rats experiment. (A) Schematic representation of hyperoxia model in BPD rats. (B) Representative images of hematoxylin and eosin-stained tissue sections showing alveolarization in rats at postnatal day 14 (P14). (C and D) Representative immunohistochemical plots of ARGs in rats' lung tissue. The analysis of the ARGs positive area in IHC. **P<0.01, *P<0.05.

Discussion

Immune Disorder of BPD Diseases

BPD remains the most common complication associated with prematurity and is increasing in prevalence, most likely due to the increased survival of preterm infants.³⁵ Inflammatory activity in the lungs leads to the dysregulation of pro-inflammatory cytokines and growth factors, including the migration of inflammatory cells to the lungs.³⁶ In our results, we found higher abundance of Macrophages M1, B cells memory, T cells gamma delta, and Eosinophils, while lower abundance of Tfh cells, NK cells resting and Monocytes in BPD.

Autophagy-Dependent Mechanism Influence BPD

Autophagy is considered a fundamental regulator of the immune response and may also influence the pathogenicity and inflammatory responses in diseases.³⁷ The progression of BPD is primarily driven by the apoptosis of alveolar epithelial cells, and autophagy plays a role in apoptosis. Our research showed that most autophagy-related genes were correlated with immune infiltration in BPD samples. Our findings were accordance with other researches. For example, an experiment demonstrated that autophagy-promoting peptide promoted autophagy and prevent apoptosis caused by exposure to hyperoxia-induced BPD model.³⁸ Another study represented that AMP-activated protein kinase-induced macrophage M2 with autophagy-dependent mechanism and played a protective role in a murine model of BPD.³⁹ A report showed that autophagy acted as an important regulator of lung development and morphogenesis, contributing to the BPD phenotype when impaired.⁴⁰ These results support the general understanding that inflammatory factors with autophagy-dependent mechanism influence the progress of BPD.

Hub Genes Involve in the Progress of Autophagy

Between the BPD and control groups, our research revealed a total 73 DE-ARGs, with 43 genes being upregulated and 30 genes being downregulated. On the basis of these DE-ARGs, we incorporated Lasso and three machine-learning algorithms to screen and identify BPD markers. We identified 5 hub genes (BRIX1, JUN, PES1, NR4A1 and RRP9) that could play a pivotal role in the development of BPD, which were correlated with inflammation and autophagy. For example, the lung macrophages expansion correlated with IL4 production and was affected by the deficiency of NR4A1.⁴¹ Translocation of the NR4A1 protein from the nucleus to mitochondria induces apoptosis. NR4A1 can regulate the inflammatory response in macrophages by regulating metabolic adaptations during inflammation, including repressing the transcription of genes involved in the citric acid cycle (TCA).^42,43 NR4A1 also inhibits NF-kappa-B signaling by binding to low-affinity NF-kappa-B binding sites, such as at the IL2 promoter.⁴⁴ Take another example, PES1 is related to rRNA processing in the nucleus and cytosol and processing cell cycle regulation. There is a study revealed that overexpression of PES1 may induce TG dysregulation by inhibiting autophagy.⁴⁵ Our results demonstrated that the expression of NR4A1 and PES1 were lower in patients with BPD compared to the control group. These results support the general understanding that BPD worsen with immune dysfunction, inhibition of autophagy and increased apoptosis.

The Relationship Between Hub Genes and Single-Cell Sequencing

Based on single-cell sequencing results, we found all hub genes mainly located in B cell, epithelial cell, fibroblast, endothelial cell, smooth muscle cell, typeIpneumocyte and typeII pneumocyte. On the one hand, we speculate that these genes promote the expression of lung fibroblasts, leading to the formation of alveolar epithelial-like tissue fibrosis in the lung, which contributes to BPD. JUN promotes pulmonary fibrosis by regulating CD36 to initiate fibrosis, increasing the number of reticular fibroblasts, and reducing the number of adipofibroblasts.⁴⁶ PES1 induces cell proliferation by altering pre-RNA.⁴⁷ We speculate that this gene promotes the proliferation of fibroblasts, leading to the production of excess fibrous tissue, which in turn promotes pulmonary fibrosis. Research related to BRIX1 and RRP9 is still limited. BRIX1 promotes ribosome synthesis through the selective translation of GLUT1.⁴⁸ We believe that excessive ribosome synthesis will lead to increased cellular physiological activity, resulting in enhanced gene transcription and protein synthesis. This, in turn, may cause immune cells or fibroblasts to excessively secrete proteins or proliferate, leading to pulmonary fibrosis. The protein–protein interaction network of RRP9 in the SSU processome, where a reduced expression of RRP9 leads to unstable interactions between U3/pre-rRNA base pairs or a decrease in the length of their connecting segments, suggesting may cause abnormal cell proliferation and promote the proliferation of fibrosis-related cells.⁴⁹ NR4A1 in human myofibroblasts reduces TGFβ1-induced collagen deposition and the expression of fibrosis-related genes,⁵⁰ which slow down the progression of BPD lung fibrosis. On the other hand, in our cord blood array analysis, five hub genes were lower in BPD patients, indicating genes were protective indicators of BPD. Autophagy of AT2 cells and lamellar body development enhanced the alveolarization and reduced the incidence of BPD in premature infants.⁵¹ Therefore, we supposed that our screening autophagy-related hub genes reduced the incidence of BPD by promoting autophagy of AT2 cells. However, the exact relationship between the relevant cells and the hub genes still need more researches to validate.

Signaling Pathways Mediate BPD

Our results suggested that these genes were correlated with ribosome processing. The ribosome is the primary site of protein translation, responsible for synthesizing nearly half of the cellular proteins. Autophagy is an important pathway for the regulation of ribosomal homeostasis.⁵² Our findings also showed that in the umbilical cord tissue homogenates, the activation of MAPK and IL 17 signaling pathway was higher in neonates with BPD compared to healthy controls. These findings corroborate a previous study showing that IncRNA H19 implicated in the progression of BPD via enhanced MAPK signaling pathway.⁵³ Another study showed that in immature rats chronic airway exposure to LPS induced pulmonary inflammation, increased IL 17a expression and hypoalveolarization.⁵⁴ These results, to some extent, provide additional information regarding the pathophysiology of BPD.

Limitation of Our Research

However, our research still has some limitations. First, we discovered one of the largest BPD datasets for analysis, but the number of included samples was still small. Therefore, a large cohort is required to corroborate our results. There were not enough in experiments to verify the level of expression of the genes identified, thus reducing the precision of our study. Consequently, further validations are required to make the findings more convincing.

Conclusions

In summary, our study results suggest that autophagy-related genes may serve as biological markers for BPD, potentially regulating immune infiltration, providing new insights and approaches for the early diagnosis and potential clinical treatment of BPD.