Circular RNAs expression profiles and bioinformatics analysis in bronchopulmonary dysplasia

Abstract Background Bronchopulmonary dysplasia (BPD) has long been considered the most challenging chronic lung disease for neonatologists and researchers due to its complex pathological mechanisms and difficulty in prediction. Growing evidence indicates that BPD is associated with the dysregulation of circular RNAs (circRNAs). Therefore, we aimed to explore the expression profiles of circRNAs and investigate the underlying molecular network associated with BPD. Methods Peripheral blood was collected from very‐low‐birth‐weight (VLBW) infants at 5–8 days of life to extract PBMCs. Microarray analysis and qRT‐PCR tests were performed to determine the differentially expressed circRNAs (DEcircRNAs) between BPD and non‐BPD VLBW infants. Simultaneous analysis of GSE32472 was conducted to obtain differentially expressed mRNAs (DEmRNA) from BPD infants. The miRNAs were predicted by DEcircRNAs and DEmRNAs of upregulated, respectively, and then screened for overlapping ones. GO and KEGG analysis was performed following construction of the competing endogenous RNA regulatory network (ceRNA) for further investigation. Results A total of 65 circRNAs (52 upregulated and 13 downregulated) were identified as DEcircRNAs between the two groups (FC >2.0 and p.adj <0.05). As a result, the ceRNA network was constructed based on three upregulated DEcircRNAs validated by qRT‐PCR (hsa_circ_0007054, hsa_circ_0057950, and hsa_circ_0120151). Bioinformatics analysis indicated these DEcircRNAs participated in response to stimulus, IL‐1 receptor activation, neutrophil activation, and metabolic pathways. Conclusions In VLBW infants with a high risk for developing BPD, the circRNA expression profiles in PBMCs were significantly altered in the early post‐birth period, suggesting immune dysregulation caused by infection and inflammatory response already existed.


| INTRODUC TI ON
Bronchopulmonary dysplasia (BPD) is the most common chronic lung disease in neonates, especially with a high incidence rate of approximately 45% in preterm infants with a gestational age of less than 29 weeks, 1 leading to severe complications and poor long-term outcomes, including impaired pulmonary function and neurodevelopmental disorders. 2,3 The main pathological changes of BPD include simplified alveolar structures and pulmonary capillary dysplasia. In addition to barotrauma and volutrauma, oxygen toxicity, fluid overload, and genetic factors, the immature lung may also suffer damage due to infection or inflammatory response. 4,5 A complex etiology and perplexing pathogenesis have led to an absence of specific biomarkers for prediction and precision therapy of BPD. 6 The rapid advancement of high-throughput microarray technologies and genome-wide sequencing techniques offers the potential for prediction, diagnosis, therapy, and pathologic insights into human diseases. [7][8][9] Mounting evidence has revealed that lncRNAs and circRNAs are critical regulators of gene expression. In addition to encoding protein, circRNAs could act as sponges for miRNAs to regulate the expression of downstream mRNAs. 10-12 A significant advantage of circRNAs over other RNAs is their more stable structure 13 and a longer half-life in various tissues and body fluids, making them more amenable to being specific biomarkers and therapeutic targets. 14,15 Considering this, mounting studies have focused initially on the epigenetic perspective and, more recently, on circRNAs to explore the molecular mechanisms of BPD. 16 Several studies have discovered that aberrant expression of circRNAs is associated with the pathogenesis of BPD. [17][18][19] First, Zou Z et al constructed a mice model of LPS-induced lung injury to identify the different expression profiles of circRNAs. This study observed numerous significantly upregulated circRNAs accompanied by abnormally increasing TNFα and IL-1β in the serum. Astonishingly, after applying P2X7 receptor antagonists, the expression of circRNAs (circ_0001679 and circ_0001212) and the cytokine mRNA levels were restrained, alleviating lung damage caused by sepsis. 20 Another study demonstrated that circRNAs could act as miRNA sponges and subsequently regulate target gene expression correlated with lung injury. In neodymium oxide-treated human bronchial epithelial (16HBE) cells, hsa_circ_0000638 inhibited the activation of NF-κB by competitively binding to miR-498-5p and further downregulated the expression of IL-8 and IL-1β. 18 Additionally, in the pulmonary hypertension mice models, Zhang J et al identified the ceRNA network and further verified that the circ-calm4/miR-337-3p/Myo10 signal transduction axis could modulate the proliferation of pulmonary artery smooth muscle cells. 21 Although partially circRNAs have been implicated in the pathogenesis of BPD, further research is necessary to confirm these findings. 19 Due to the difficulty of obtaining human lung tissues, peripheral blood is one of the most common substitution sources. 22 A growing number of studies have recently employed gene expression profiles obtained from peripheral blood mononuclear cells (PBMCs) as both a method for pathogenesis research and a diagnostic tool for several lung diseases. 23 As part of this study, we used microarray technology to detect the circRNA expression profiles in PBMCs between BPD and normal VLBW infants. Based on the circRNA expression profiles, DEcircRNAs were identified, and a portion with significant upregulation was further validated with qRT-PCR. Simultaneously, we obtained the DEmRNA expression profiles between BPD and non-BPD VLBW infants from GEO databases. Potential miRNAs binding to circRNAs and mRNAs were predicted for further construction of circRNA-miRNA pairs and miRNA-mRNA pairs, respectively. Finally, we constructed the circRNA-miRNA-mRNA regulatory network based on screening circRNA-miRNA pairs and miRNA-mRNA pairs with overlapped miRNAs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to investigate the molecular function and signaling pathways affecting BPD progression. These identified DEcircRNAs and DEmRNAs obtained from moderate to severe BPD infants may provide novel insights for early disease prediction and management.
The flowchart of this study is shown in Figure 1.

| Specimen collection
PBMCs were derived from peripheral blood samples obtained from VLBW infants (birth weight <1500 g and gestational age < 32 weeks) born in Suzhou Municipal Hospital from May to December 2020. For each VLBW infant, 0.5 ml peripheral blood specimen was collected 5-8 days after birth. The PMBCs were separated using lymphocyte isolation solution and transferred into enzyme-free 1.5 ml EP tubes.
Samples were immediately frozen at −80°C after adding the appropriate amount of TRlzol reagent (Invitrogen, United States). The entire manipulation process was completed within 2 h.
The categorized grade of BPD was defined based on the infant's requirement for oxygen and different respiratory support modalities at 36 weeks corrected gestational age (moderate BPD: oxygen requirement <30% FiO 2 ; severe BPD: oxygen requirement >30% FiO 2 or continuous positive airway pressure (CPAP)/mechanical ventilation). [24][25][26] Exclusion criteria were (1) infants with chromosomal abnormalities, (2) significant congenital deformities, and (3) treated with corticosteroids in the first week of birth. Finally, three VLBW infants with moderate to severe BPD were selected as the BPD group, and three VLBW infants with a matched gestational age of non-BPD were chosen as the control group. Basic information about these infants is shown in Table 1.

| CircRNA microarray analysis
Total RNA was extracted using the TRLzol reagent (Invitrogen, United States). The quality and quantity of the RNA samples were detected with a NanoDrop ND-1000 spectrophotometer and agarose gel electrophoresis. Only samples with OD260/280 ratios be-

| Quantitative real-time RT-PCR validation
We selected five significantly upregulated DEcircRNAs (hsa_ All measurements were performed in triplicate, and the gene expression levels were standardized to β-actin. The primers were synthesized by Genscript Corporation, and the primer sequences are listed in Table 2.

| Microarray Data acquisition from publicly available database
We screened the GEO database for microarray sequencing analysis data that met the following criteria: (1) studies containing preterm infants of less than 32 weeks gestational age and birth weight ≤1500 g, (2) studies with blood samples that were drawn from infants on the 4th to 10th day of life for the assessment of gene expression in peripheral blood leukocytes, and (3)

| Identification of DEmRNAs and DEmiRNAs
First, data downloaded from the GEO were preprocessed and normalized using the preprocess Core package in R version 3.6.2 ( Figure S2). Then, between the moderate to severe BPD group and non-BPD group, mRNAs with a fold change >2.0 and adjusted p-value <0.05 were identified as differentially expressed (DEmRNAs).
We integrated the Circular RNA Interactome database to identify the downstream target miRNAs for each DEcircRNA and  and the mirDIP database. 30 Following that, we used Venny 2.1 (https://bioin fogp.cnb.csic.es/tools/ venny/) to make a Venn diagram and screened for overlapping miRNAs common to these data sets.

| Construction of ceRNA network
Based on the DEcircRNAs identified after qRT-PCR validation, the overlapped prediction miRNAs, and the DEmRNAs, we constructed the circRNA-miRNA pairs and the miRNA-mRNA pairs.
Subsequently, a novel circRNA-miRNA-mRNA regulatory network based on the ceRNA hypothesis was constructed using Cytoscape 3.7.1.

| GO and KEGG functional enrichment and statistical analysis
To investigate the molecular interaction and mechanism underlying the ceRNA regulatory network, we performed GO and KEGG functional enrichment analysis based on g: Profiler. 31

| General clinical data
Among the six VLBW infants, there was no significant difference in birth weight or gestational age between the BPD and control groups. For the BPD group, the median gestational age at birth Additional clinical character details of the six VLBW infants are shown in Table 1.
From the GEO database of GSE32472, we screened a total of 28 VLBW infants with moderate to severe BPD and 35 VLBW infants with non-BPD based on the time of blood collection and the diagnostic grading criteria for BPD. However, the BPD group's gestational age and birth weight were significantly lower than the control group, except for gender. BPD group, the median gestational age at birth was 25.5 (24.5-27) weeks, and the non-BPD group was 29 [28][29][30] weeks; similarly, the median birth weight of the BPD group was 705 (600-910) g, and the control group was 1200 (1150-1400) g, reveal- ing the uneven distribution of both gestational age and birth weight (p-Value = 0).

| CircRNA expression profiles
Among all the expression dysregulation circRNAs of the BPD group, 65 circRNAs were confirmed to be significantly differentially downregulated. We compared the samples' post-normalization expression value distributions using a box plot. Each sample's log2 ratio followed the same pattern. We noticed that most dysregulated circRNAs were generated by exons; the others were derived from introns and intergenic regions. In the stacked bar plot, the distribution of circRNAs on different human chromosomes was revealed ( Figure 2). Volcano plots, heatmap of hierarchical clustering analysis, and scatter plots ( Figure S1B) were used to illustrate the differences in circRNA expression between the two groups. The top 10 up or downregulated circRNAs are summarized in Table 3A,B.

| Quantitative Real-time RT-PCR validation
Results indicated that more upregulated circRNAs were expressed in the BPD group than downregulated ones. Thus, to confirm the findings of the circRNA microarray analysis, qRT-PCR was utilized to validate the expression of five significantly upregulated circR-NAs. Finally, four circRNAs (hsa_circ_0007054, hsa_circ_0050386, hsa_circ_0057950, and hsa_circ_0120151) were consistent with the microarray data ( Figure 3); nevertheless, hsa_circ_0057953 did not amplify by qPCR.

| GO and KEGG functional enrichment analysis
GO and KEGG pathway analyses were used to annotate the function of DEmRNAs in the ceRNA network. A total of 408 GO terms and 9 KEGG pathways were significantly enriched (p.adj <0.05), and the top 20 most significantly enriched terms and pathways are listed in Figure 6. GO biological process (BP) was mainly enriched in "immune system process," "response to external stimulus," "myeloid leukocyte activation," "response to biotic stimulus," and "neutrophil activation."

F I G U R E 3
The expression levels of the circRNAs validation by qRT-PCR Gene expression was assessed using ΔCt values, normalized to β-actin levels. Means ± SEM was presented in the table. *p < 0.05.

| DISCUSS ION
Over the past several decades, neonatal health care workers and and fungi) were demonstrated with a significantly higher risk of BPD compared with those without sepsis, with the relative risk (RR) ranging from 2.6 (95% CI: 1.5-4.6) to 9.40 (95% CI: 3.83-23.08) among different studies. 41,42 All these suggest that the inflammation response initially triggered by diverse pathogens could contribute to BPD progression and lead to long-term adverse outcomes.
Alternatively, receptor binding may activate many signaling Additionally, we employed KEGG pathway analysis to explore notable signaling pathways in the ceRNA network. Aside from F I G U R E 6 GO and KEGG functional enrichment analysis For GO enrichment analysis, three components were employed for demonstrating different levels of biological functions: the molecular function (MF), the biological process (BP), and the cellular component (CC). In addition, the KEGG enrichment analysis was applied to assess the level of gene enrichment in various pathways. Based on the most significant adjusted p-Value, the top 5 GO terms and KEGG pathways were selected (p.adj <0.05).
the fact that the renin-angiotensin system was involved in pulmonary fibrosis and the development of BPD under hyperoxia exposure, 47 another interesting finding was the enriched metabolic pathway. We inferred these were linked to lung microbes due to more recent data revealing that the abnormal distribution of the microbiota in the lung is inextricably linked to the progression of BPD. 40 53 In contrast to earlier results, however, no evidence of these signaling pathways was detected in our study.
We inferred two main factors might be responsible for this deviation. First, the development of BPD is multidimensional, and the progression of the disease is also dynamic. 39 In this study, we col- day compared with those on the 5th day postnatally in peripheral blood. 27 In our study, the significantly enriched pathways during the early stage were most inflammation-related, consistent with the disease progression pattern of BPD. Another reason for this discrepancy was the difficulty obtaining human pulmonary tissue other than peripheral blood samples. The gene expression profiles of peripheral blood leucocytes cannot entirely replace those of lung tissue. However, some studies demonstrated that proinflammatory cytokines in peripheral blood were associated with the phenotype of BPD. 54 Numerous recent studies have already used this approach as an alternative to lung tissues to identify novel biomarkers for BPD. 27 We acknowledge the limitations of our research. First, in this study, differentially expressed circRNAs were identified in a limited number of infants and should be confirmed in a larger cohort.
Additionally, many biological processes' crucial changes are not only manifested by changes in RNA levels but also at DNA and protein levels. 55 These molecular alterations need to be validated further.
Even so, the exact role and molecular mechanism of circRNAs affecting the progression of BPD are yet to be determined.

| CON CLUS ION
We performed microarray technology to identify different expression profiles of circRNAs in peripheral blood from moderate to severe BPD infants 5-8 days after birth. Hsa_circ_0007054, hsa_circ_0057950, hsa_circ_0050386, and hsa_circ_0120151 were verified by qRT-PCR to be differentially expressed between the two groups. Based on mining data from public databases, we constructed a ceRNA regulatory network, including 3 circRNAs, 22 miRNAs, and 16 mRNAs. Gene ontology enrichment and KEGG pathway analysis revealed that these genes were associated with immune cell activation and cytokine receptor activation, which promote the inflammatory response and lung injury. 56

CO N FLI C T O F I NTE R E S T
For this study, Yang Zuming was fortunate to get financial assistance from the Jiangsu Commission of Health and Yu Lun got fund from the Suzhou Commission of health. Yu Lun, Junlong Hu, and Yang Zuming have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
Due to the sensitivity of human data, the data that support the findings of this study are not publicly available, but are available from the corresponding author upon reasonable request.