Next-generation sequencing to investigate circular RNA profiles in the peripheral blood of preterm neonates with bronchopulmonary dysplasia.

BACKGROUND: Circular RNAs (circRNAs) are emerging noncoding RNAs that are involved in many biological processes and diseases. The expression profile of circRNAs in preterm neonates with bronchopulmonary dysplasia (BPD) remains unresolved.
METHODS: In BPD infants, peripheral venous blood was drawn and circRNAs were extracted and sequenced by next-generation sequencing. The levels of the selected circRNAs were measured by real-time quantitative reverse transcription PCR.
RESULTS: Among thousands of circRNAs, 491 circRNAs were significantly changed. Among the top 10 changed circRNAs, hsa_circ_0003122, hsa_circ_0003357, hsa_circ_0009983, hsa_circ_0003037, and hsa_circ_0009256 were significantly increased, while hsa_circ_0014932, hsa_circ_0015109, hsa_circ_0017811, hsa_circ_0020588, and hsa_circ_0015066 were significantly decreased. These altered circRNAs are involved in complicated biological functions and signaling pathways. Additionally, hsa_circ_0005577 (hsa_circ_FANCL), which was significantly increased in the moderate-to-severe BPD subjects, was correlated with oxygenation therapy.
CONCLUSION: These results suggest that an aberrant circRNA profile in the peripheral blood of BPD infants might be important in BPD pathogenesis.

INTRODUCTION

Bronchopulmonary dysplasia (BPD) is one of the most common chronic lung diseases in infants.1 With arrested pulmonary development, the alveolarization and vascularization of the lung are severely reduced. Although hyper‐oxygen exposure, mechanical ventilation, and congenital infection have been postulated as initiators of the disease,2, 3 the pathogenesis of BPD is still largely unknown.

Noncoding RNAs have emerged as potential key players in the pathogenesis of BPD. MicroRNAs (miRNAs) are 22‐nucleotide long noncoding RNAs that interact with the 3’ untranslated region to promote mRNA degradation. In animal models of BPD, several miRNAs in lung tissues showed differential expression during the process of BPD.4 Long noncoding RNAs (lncRNAs) exceed 200 nucleotides in length and exert various effects on lung development. In the animal model of BPD, lncRNA profiles in lung tissues have been documented.5, 6 Different from miRNAs and lncRNAs, circular RNAs (circRNAs) have emerged as a novel type of noncoding RNAs that form covalently closed continuous loops without 5′ to 3′ polarities or poly (A) tails.7 With advances in high‐throughput technology, thousands of circRNAs have been documented in lung diseases, especially in lung cancer.8, 9 The roles of circRNAs in noncancerous pulmonary diseases remain largely unclear.

In the present study, we explored the aberrant circRNA profile in the blood of BPD infants. The identification of circRNAs closely associated with BPD may shed light on the pathogenesis of BPD.

MATERIALS AND METHODS

BPD infants and blood samples

This study was carried out in accordance with the recommendations of the Committee on Human Rights Related to Research Involving Human Subjects, Faculty of Affiliated Children Hospital, Nanjing Medical University (Number: NJCH2016003). Written informed consent was obtained from the parents of the infants in this study. Preterm infants with or without BPD born at gestational age ≤32 weeks were included in the study. To screen the circRNA profile, 3 mL peripheral venous blood samples were drawn from 3 BPD infants and 3 non‐BPD controls. To validate the expression of the selected hsa_circ_0005577 (hsa_circ_FANCL), peripheral blood samples were collected from 16 BPD infants and 14 non‐BPD controls. Each sample was mixed with TRIzol and frozen in liquid nitrogen until it was used for circRNA sequencing or quantification. The clinical characteristics of BPD patients and controls are presented in Table 1 and Table 2. Infants suffering from severe infections, shock, inherited metabolic diseases, or severe congenital malformations were excluded from the study. The infants were diagnosed with BPD according to the National Institutes of Child Health and Human Development/National Heart, Lung, and Blood Institute and Office of Rare Diseases workshop definitions.10

Table 1

Clinical characteristics of the BPD and non‐BPD infants in the RNA‐seq analysis

	Non‐BPD (n = 3)			BPD (n = 3)			P‐value
	#1	#2	#3	#1	#2	#3
Infants' characteristics
Sex gender	Male	Male	Male	Female	Male	Female	ND
Birthweight (grams)	1990	2200	3000	1170	970	1150	.0144
Gestational age (d)	251	249	246	190	186	205	.0008
Apgar 1 min	9	10	10	7	6	10	.1841
Apgar 5 min	10	10	10	8	9	10	.1583
Intraventricular hemorrhage (IVH)	Yes	Yes	Yes	Yes	Yes	Yes	ND
Necrotizing enterocolitis (NEC)	No	No	No	No	No	No	ND
Late‐onset neonatal sepsis (LOS)	No	No	No	No	No	Yes	ND
Mechanical ventilation (d)	0	0	0	19	43	29	.0121
CPAP (d)	0	0	0	10	5	6	.0102
Days with oxygen	0	0	0	14	27	11	.0242
Surfactant treatment	No	No	No	Yes	Yes	No	ND
Patent ductus arteriosus (PDA)	Yes	Yes	Yes	No	No	No	ND
Hospitalization days	15	8	10	71	76	94	.0007
Maternal characteristics
Preeclampsia	No	No	No	No	No	No	ND
Antenatal steroids	No	No	Yes	No	Yes	No	ND
Premature rupture of membranes (PROM)	No	Yes	No	Yes	No	No	ND

Abbreviation: ND, not done.

Table 2

Clinical characteristics of the BPD and non‐BPD infants in the qRT‐PCR quantification of circular RNAs

	Non‐BPD (n = 14)		BPD (n = 16)		P‐value
	Mean ± SEM	n (%)	Mean ± SEM	n (%)
Infants' characteristics
Male gender		10 (71)		8 (50)	.4112
Birthweight (grams)	1452 ± 72.81		1311 ± 166.1		.4648
Gestational age (d)	207.1 ± 2.832		199.1 ± 5.259		.2047
Apgar 1 min	8.200 ± 0.1333		7.500 ± 0.3743		.1419
Apgar 5 min	9.000 ± 0.1491		8.071 ± 0.3701		.0549
Intraventricular hemorrhage (IVH)		14 (100)		16 (100)	1.0000
Necrotizing enterocolitis (NEC)		2 (14.3)		3 (18.8)	.8700
Late‐onset neonatal sepsis (LOS)		3 (21.4)		10 (62.5)	.0580
Mechanical ventilation (MV) (d)	0.1429 ± 0.1429		28.07 ± 7.084		.0007
CPAP (d)	1.214 ± 0.8266		15.27 ± 3.984		.0024
Days with oxygen	3.350 ± 1.029		14.40 ± 2.760		.0002
Days with MV + CPAP+Oxygen	3.938 ± 1.769		57.73 ± 7.234		<.0001
Surfactant treatment		1 (7.1)		12 (75.0)	<.0001
Patent ductus arteriosus (PDA)		8 (57.1)		12 (75.0)	.5177
Hospitalization days	20.10 ± 3.497		63.13 ± 7.597		<.0001
Maternal characteristics
Preeclampsia		1 (7.1)		1 (6.2)	.5249
Antenatal steroids		9 (64.2)		9 (56.2)	.9405
Premature rupture of membranes (PROM)		4 (28.5)		4 (25.0)	.8469

RNA preparation and circular RNA sequencing

Total RNA from whole peripheral venous blood was extracted using TRIzol reagent (Invitrogen) following the manufacturer's procedure. The RNA amount and purity of each sample were quantified using a NanoDrop ND‐1000 (NanoDrop). A Ribo‐Zero™ rRNA Removal Kit (Illumina) was used to deplete ribosomal RNA in approximately 5 μg of total RNA. Then, the remaining RNA samples were treated with RNase R (Epicentre Inc) to remove linear RNAs and to enrich circRNAs. The cleaved RNA fragments from the enriched circRNAs were reverse‐transcribed to create cDNA. The average insert size for the final cDNA library was 300 ± 50 bp. Finally, we performed paired‐end sequencing on an Illumina HiSeq 4000 (LC Bio). The reads that contained adaptor contamination, low‐quality bases, and undetermined bases were removed with Cutadapt.11 The sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Bowtie212 and TopHat213 were used to map reads to the reference genome. The remaining reads (unmapped reads) were still mapped to the genome using TopHat‐Fusion.14 CIRCexplorer15, 16 was used to de novo assemble the mapped reads to circRNAs first; then, back‐spliced reads were identified among the unmapped reads by TopHat‐Fusion and CIRCexplorer. All samples generated unique circRNAs. The differentially expressed circRNAs were selected with log2 (fold change) >1 or log2 (fold change) <−1 and with statistical significance (P‐value < .05) by t test.

RNA quantification

Real‐time quantitative reverse transcription PCR was performed to confirm circRNA expression in peripheral venous blood samples. Briefly, circRNAs were extracted from peripheral venous blood samples using the miRNeasy Mini Kit (217004, QIAGEN), reverse‐transcribed into cDNA using the TUREscript 1st stand cDNA synthesis kit (Aidlab), and quantified using SYBR Premix Ex Taq™ II (638515, TAKARA). The primers are listed in Table 3. The cycling conditions were 95°C for 10 minutes, followed by 95°C for 15 seconds and 60°C for 60 seconds for up to 40 cycles. The relative expression levels of the circRNAs were normalized to that of the internal control GAPDH by using the 2−ΔΔCt cycle threshold method.

Table 3

Primers in the circular RNA quantification

Primers	Sequence (5′ to 3′)
GAPDH‐F	AAAGGGTCATCATCTCTG
GAPDH‐R	GCTGTTGTCATACTTCTC
Hsa_circ_0004033‐F	ACCACCATGGCCTCCG
Hsa_circ_0004033‐R	GCTGCTTTGTGTCCAGCTTC
Hsa_circ_0005577‐F	GCACTACCTCCTCCTCCC
Hsa_circ_0005577‐R	ATCCACTAAGTATTGTTCTCAGC
Hsa_circ_0008959‐F	GGGATACAACAGGCCTTTACA
Hsa_circ_0008959‐R	GATGACTGTGAAGATCAGGC

Statistical analysis

The results for variables that were normally distributed are presented as the mean ± SEM. ANOVA was performed to establish equal variance, and a 2‐tailed Student's t test with Bonferroni correction was applied to determine statistical significance using GraphPad Prism version 5.0. Correlation analysis between hsa_circ_0005577 in peripheral blood and patient characteristics was performed using Spearman rank correlation in MedCalc Software version 15.0. The ROC curve of hsa_circ_0005577 in the diagnosis of BPD was generated by MedCalc Software version 15.0.

RESULTS

Circular RNA profiles of BPD infants

RNAs from whole peripheral blood were extracted, and circRNAs were assayed. Of the 60 498 sequenced genes, 23 262 were annotated as circRNAs. Most of these circRNAs (97.11% ~ 99.4%) were mapped to exons. A volcano map showed that the expression levels of hundreds of circRNAs were altered. Among them, 127 circRNAs were significantly increased and 364 circRNAs were significantly decreased in the BPD subjects. Furthermore, heatmap analysis revealed that the inner‐group variation was small, suggesting that these alterations may be common in BPD pathogenesis (Figure 1).

Figure 1

Blood circRNA profile of BPD infants. A, Volcano plots presenting the differential expression of circRNAs. The vertical line corresponds to the fold change, and the horizontal line represents the p‐value. The red dots are the differentially expressed circRNAs with statistical significance (log2 (fold change) >1 and a P‐value < .05). B, In total, 364 significantly downregulated circRNAs and 127 upregulated circRNAs were observed. C, Heatmap and hierarchical clustering analysis of the altered circRNAs

Hsa_circ_0005577 was increased in the moderate‐to‐severe BPD infants

As described above, we observed that the circRNA profile from peripheral blood was altered in BPD infants. Based on the fold change of circRNAs between the BPD and non‐BPD subjects, the top 10 circRNAs with significant changes were selected. Among them, hsa_circ_0003122, hsa_circ_0003357, hsa_circ_0009983, hsa_circ_0003037, and hsa_circ_0009256 were significantly increased, while hsa_circ_0014932, hsa_circ_0015109, hsa_circ_0017811, hsa_circ_0020588, and hsa_circ_0015066 were significantly decreased. In the analysis step, we set cutoff value to avoid the low abundant circular RNAs, and these most changed 10 circRNAs, however, were not most abundant in expression. Considering the importance of expression levels in RNA sequencing, we selected the 3 circRNAs with the top fragments per kilobase of transcript per million mapped reads (FPKM) (Table 4). In the quantitative RT‐PCR analysis from the BPD and non‐BPD samples, the expression of hsa_circ_0004033, hsa_circ_0008959, and hsa_circ_0005577 was comparable. The disease severity of BPD was defined as suggested by the bronchopulmonary dysplasia workshop summary.10 The results demonstrated that hsa_circ_0005577 was significantly increased in the moderate‐to‐severe BPD preterm infants (Figure 2). Arguably, hsa_circ_0005577 may be associated with the disease severity of BPD.

Table 4

Circular RNAs for qRT‐PCR

Accession	chr	Start	End	Strand	Exon count	circType	Isoform name	geneName
Hsa_circ_0005577	chr2	58221942	58232112	−	4	circRNA	ENST00000403295	FANCL
Hsa_circ_0005989	chr1	10134988	10137205	+	2	circRNA	ENST00000253251	UBE4B
Hsa_circ_0004033	chr7	157186834	157187021	+	1	circRNA	ENST00000348165	UBE3C

Figure 2

Quantification of selected circRNAs. A‐C, The expression levels of hsa_circ_0004033, hsa_circ_0008959, and hsa_circ_0005577 were comparable between BPD patients and non‐BPD controls. D, Compared with that in the non‐BPD controls, hsa_circ_0005577 was significantly increased in the moderate‐to‐severe BPD preterm infants (P = .0324). E, Hsa_circ_0005577 was significantly increased in the moderate‐to‐severe BPD preterm infants compared with the mild BPD patients (P = .0137)

Hsa_circ_0005577 expression was associated with oxygen therapy

BPD is usually diagnosed as oxygen dependence at 36 weeks of postmenstrual age. Considering that hsa_circ_0005577 may be associated with disease severity, we explored whether hsa_circ_0005577 contributed to oxygen therapy. Preterm neonates may acquire oxygen through mechanical ventilation, CPAP, or oxygenation. A positive correlation between hsa_circ_0005577 and mechanical ventilation, continuous positive airway pressure (CPAP), or oxygenation was established. Although the P‐value in the individual correlation analysis did not reach significance, hsa_circ_0005577 was significantly correlated with the duration of oxygen therapy (Figure 3), further suggesting that hsa_circ_0005577 may be involved in BPD pathogenesis. Next, we tested whether circRNA could be used to diagnose BPD. In the ROC analysis (Figure 3), the specificity and sensitivity were 50% and 75%, respectively. Therefore, hsa_circ_0005577 alone may not be a useful marker in BPD diagnosis.

Figure 3

Correlation analysis of hsa_circ_0005577 with oxygen therapy. A, Hsa_circ_000557, was not correlated with the mechanical ventilation days (P = .3506). B, Hsa_circ_0005577 was not correlated with CPAP days (P = .1586). C, Hsa_circ_0005577 was not correlated with the oxygenation days (P = .5233). D, Hsa_circ_0005577 was significantly associated with oxygen therapy days (P = .0407). E, ROC analysis of hsa_circ_0005577 for BPD diagnosis. The associated criterion was > 3.589, the sensitivity was 75%, and the specificity was 50%

Potential targets of hsa_circ_0005577

Hsa_circ_0005577 is back‐spliced from exon 4 of mRNA FANCL (Fanconi anemia, complementation group L), which is associated with acute lung injury.17 Therefore, hsa_circ_0005577 was renamed hsa_circ_FANCL. As noncoding RNAs, circRNAs may regulate miRNA functions.18 To predict the potential targets, analyses in TargetScan19 and miRanda20 were performed; circRNAs may sponge miRNAs with complementary sequence in the seed region.21 The potential targets of hsa_circ_FANCL were let‐7, miR‐196, miR‐20a, miR‐22, and miR‐26a among other miRNAs.

DISCUSSION

Due to improved clinical management and treatment, morality in preterm infants has greatly decreased. Consequently, the incidence of BPD, which has been attributed to lung immaturity in preterm infants, has risen.22 Oxygen requirement, positive pressure ventilation, and the immaturity of the respiratory system have been suggested to contribute to the pathogenesis of BPD. Blood tests, which quantify the oxygen saturation, might be valuable in BPD diagnosis. However, BPD diagnosis is still functional and is largely based on an oxygen requirement at 36 weeks of postmenstrual age. There is no specific diagnostic marker for BPD.23 CircRNAs, which are usually stable, abundant, and conserved, have been postulated as diagnostic markers in cancer.24, 25 Herein, we determined that the circRNA profile was significantly altered in BPD infants. In total, 491 circRNAs were significantly changed.

Among these altered circRNAs, hsa_circ_FANCL (hsa_circ_0005577) was significantly increased in the moderate‐to‐severe BPD blood samples. As competing endogenous RNA,21 hsa_circ_FANCL may sponge targeted miRNAs. Dysregulated miRNA expression has been associated with BPD. In the peripheral blood samples from BPD infants, miR‐133b and miR‐7 were significantly increased, whereas miR‐152 and miR‐30a‐3p were significantly reduced.26 In the lung tissues from BPD subjects, let‐7 and other miRNAs were significantly altered.4 In the neonatal murine model of BPD, overexpression of let‐7 could alleviate the disease, which may be associated with the downregulation of the TGF‐β pathway.27 miRNA let‐7 was shown to be reduced in human lung cancer, and overexpression of let‐7 inhibited cancer cell growth.28 The potential targets of let‐7 include RAS and HMGA2,29 which promote cancer advancement. It is interesting to note that RAS and HMGA2 are also closely associated with many other human illnesses, including immunological and inflammatory disorders.30 In HMGA2 knockout mice, the proper cell proliferation and distal epithelium differentiation during embryonic lung development were compromised.31 Therefore, hsa_circ_FANCL targeting let‐7 may contribute to BPD pathogenesis via the regulation of the TGF‐β/RAS/HMGA2 pathways.

miRNAs are considered to be conserved in different species. In contrast to miRNAs with ~22 nucleotides, lncRNAs with >200 nucleotides are less conserved among species. However, hundreds of lncRNAs are still conserved. Results from animal studies5, 6 demonstrated that Smgc and other lncRNAs were significantly altered. Moreover, the lncRNA MALAT1 protected mice with BPD by inhibiting cell apoptosis.32 We could not exclude the possibility that hsa_circ_FANCL may regulate BPD progression by binding to lncRNAs.

Our study was limited by some disadvantages. First, the clinical characteristics of BPD and non‐BPD controls in the study were uncertain. Intraventricular hemorrhage (IVH) is common in preterm neonates.33 In Tables 1 and 2, all premature infants developed grade 1 or silent IVH and were diagnosed by craniocerebral ultrasound. Neonates with severe IVH were excluded from the research. In mainland China, the incidence of IVH is as high as 70% in low‐weight preterm infants,34 which may explain the high ratio of IVH in our cohort. In addition to IVH, patent ductus arteriosus (PDA) is common (~70%) in premature infants.35 As shown in Table 1, PDA was present in the non‐BPD preterm infants but absent from BPD neonates. As shown in Table 2, with more samples included, the PAD incidence between non‐BPD and BPD preterm babies was comparable. BPD is a heterogeneous disease with a complicated etiology and different clinical phenotypes.36, 37 Therefore, more samples and definitive criteria in the BPD diagnosis should be evaluated in the study. Second, only hsa_circ_FANCL was selected in the qRT‐PCR analysis. Hsa_circ_FANCL may contribute to BPD disease severity and pathogenesis. However, hsa_circ_FANCL alone could not be used to diagnose BPD. The combination of several different circRNAs may provide better specificity and sensitivity for BPD diagnosis. Moreover, targets of hsa_circ_FANCL were only predicted using two different methods. Further studies are needed to verify potential targets of hsa_circ_FANCL.

CONCLUSIONS

In summary, we identified a circRNA profile in the peripheral blood of BPD infants and found that 491 circRNAs were significantly changed. Hsa_circ_FANCL, which was significantly increased in the moderate‐to‐severe BPD preterm infants, may regulate BPD progression by acting as a sponge of miRNAs. In the future, we would like to evaluate the value of the circRNA profile in BPD diagnosis and the mechanisms of action of hsa_circ_FANCL in the pathogenesis of BPD.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.