Differential expression of circRNAs during rat lung development.

At present, thousands of circular RNAs (circRNAs) have been found in cancer and various tissues from different species. However, the expression of circRNAs during rat lung development remains largely unknown. In the present study, circRNA expression profiles were screened in three mixed rat lung tissues at 3 time‑points [embryonic day (E) 19, E21 and post‑natal (P) day 3] during fetal rat development with circRNA high‑throughput sequencing. Preliminary results were verified by reverse transcription‑PCR (RT‑PCR) at 4 time‑points (E16, E19, E21 and P3). A total of 375 circRNAs were differently expressed in E19 vs. E21 (fold change ≥1.5; P<0.05). At the same time, a total of 358 circRNAs were differently expressed in E21 vs. P3 (fold change ≥1.5; P<0.05). A total of 3 circRNAs (rno_circ:chr7:24777879‑24784993, rno_circ:chr14:14620910‑14624933 and rno_circ:chr3:1988750‑​1998592) were characterized by having consistent fold changes (≥1.5) between 3 time‑points (E19, E21 and P3) and were selected for RT‑PCR at 4 time‑points (E16, E19, E21 and P3). Subsequently, Gene Ontology and Kyoto Encyclopedia of Genes and Genomes analysis of parent genes of the differentially expressed circRNAs revealed that these circRNAs may serve important roles in lung development. The present results support that these new found circRNAs participate in lung development. Furthermore, these findings may help to clarify the physiopathological mechanism of normal rat lung development, and may further provide a physiopathological basis of lung developmental diseases.

Introduction

Lung development is the physiological basis of breathing. Although, at present, the mortality rate of premature infants has decreased, the incidence of lung developmental diseases, such as bronchopulmonary dysplasia (BPD) and respiratory distress syndrome (RDS) is increasing (1,2). Lung development is a complex process that extends from the embryonic period to the fetal period to birth, and involves cell proliferation, differentiation and morphogenesis (3). The classic molecular pathways include the Wnt, fibroblast growth factor (FGF) and bone morphogenetic protein pathways, and the microRNAs (miRNAs/miRs) involved in regulating lung development include miR-7, let-7, miR-221, miR-30 and miR-93 (4,5).

A number of circular RNAs (circRNAs) have been discovered in different kinds of cell lines and species (6). This small RNA is a novel class of non-coding RNAs characterized by a covalently closed-loop structure without a 5′-terminal cap or 3′-terminal poly tail, which is generated through a special type of alternative splicing, termed backsplicing (7). Generally, there are three types of circRNA, including exonic, intronic and retained-intronic circRNA (8). Accumulating evidence has demonstrated that circRNAs can act as both an miRNA sponge and RNA blinding protein sponge, which regulates alternative splicing and gene transcription (8-10). This suggested that the sponge activity of circRNA affects the pathway regulation by targeting miRNAs. Furthermore, previous studies have shown that circRNAs have a spatiotemporal specificity during species development. It was reported that circRNAs are differentially expressed during nervous system development in drosophila (11), porcine embryos (12) and human fetuses (13). However, the role of circRNAs in lung development remains unknown.

The 5 stages of lung development consist of the embryonic (first 13 days), pseudoglandular (13-18 days), cana-licular (18-20 days), saccular (20 days-full term) and alveolar (following birth) phases (14). During the canalicular phase, lamellar bodies begin to emerge, marking a significant point in lung maturation. During the saccular phase, alveolar ducts and air sacs form, an event that is critical for postnatal gas exchange (14). During the alveolar phase there is complete alveoli formation (14,15). Therefore, in the present study, 3 representative time points [embryonic day (E) 19, E21 and postnatal (P) day 3] were selected for circRNA high-throughput sequencing. In addition, 4 time points (E16, E19, E21 and P3) were selected for reverse transcription-PCR (RT-PCR) verification to provide further insight for rat lung development. To the best of our knowledge, circRNA high-throughput sequencing was used for the first time to perform differential expression profiling of circRNAs during rat lung development.

Materials and methods

Rat fetal lung tissues

A total of 12 10-week-old healthy pregnant female Sprague-Dawley rats (270-360 g) obtained from The Animal Center of Nanjing Medical University, were raised in a specific-pathogen-free animal facility at 21-25°C, with a humidity 45-65%, a 12 h light/dark cycle, and ordinary feed and drinking water at The Animal Center of Nanjing Medical University. The duration of the gestation period for the rat strain used in the present study was 23 days. The animal use protocol has been reviewed and approved by The Nanjing Medical University Animal Ethical and Welfare Committee (approval no. IACUC-1809020). All pregnant SD rats were sacrificed by cervical dislocation after anesthesia with 2% chloral hydrate (400 mg/kg), and then all the rat fetuses at E16, E19 and E21 were sacrificed with direct decapitation, while the rat fetuses at P3 were sacrificed with direct cervical dislocation. The whole fetal lungs were isolated from the rat fetuses on E16, E19, E21 and P3. These 4 time points were named as the E16 group, E19 group, E21 group and P3 group. There were 3 pregnant rats in each group following a random contrast rule. The lung tissues were washed with normal saline after isolation. In total, 1 sample of fetal lungs was randomly selected for morphological observation in each group, E19, E21 and P3. In the E19, E21 and P3 groups, 3 samples of fetal lungs in each group were mixed to produce 1 sample to be used in high-throughput sequencing. In the E16, E19, E21 and P3 groups, 10 samples of fetal lungs were used for RT-PCR.

Hematoxylin and eosin staining

The lung tissues were fixed with 4% paraformaldehyde for 24 h at 4°C, embedded in paraffin, and cut into longitudinal sections of 4 µm. Then, staining and mounting were performed, according to the method of hematoxylin and eosin staining (16). The lung tissues were observed under a light microscope (DM2500; Leica Microsystems GmbH; magnifications, ×20 and 40).

Library construction and sequencing

Total RNA was isolated from tissues using TRIzol® (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. RNA purity was assessed using ND-1000 Nanodrop (Thermo Fisher Scientific, Inc.) requiring A260/280 ≥1.8, A260/A230 ≥2.0. RNA integrity (RIN) was evaluated using Agilent 2200 TapeStation (Agilent Technologies, Inc.) requiring RIN ≥7.0. Briefly, ribosomal RNAs were removed from total RNA using the Epicentre Ribo-Zero rRNA Removal kit (Illumina, Inc.). Then, RNA was treated with RNase R (Epicentre; Illumina, Inc.) and fragmented to ~200 bp. Subsequently, the purified RNA fragments were subjected to first strand and second strand cDNA synthesis. The RT temperature protocol consisted of the first strand cDNA synthesis (10 min at 25°C, 15 min at 42°C, 15 min at 70°C and hold at 4°C) and the second strand cDNA synthesis (incubate in a thermal cycler for 1 h at 16°C, with heated lid set at 40°C), following adaptor ligation and enrichment with a low-cycle using the NEBNext® Ultra™ RNA Library Prep kit for Illumina (NEB; Illumina, Inc.) according to the manufacturer's protocol. The purified library products were evaluated using the Agilent 2200 TapeStation and Qubit® 2.0 (Thermo Fisher Scientific, Inc.) and then sequenced on HiSeq 3000 with 2×150 bp mode.

Pre-processing of sequencing reads/quality control

Raw reads were treated with Trimmomatic tools (V0.36) to remove adapters (17). Following reads quality control: Scan the read with a 4-base wide sliding window, cutting when the average quality per base drops below 15, drop reads which are less than 35% of initiation read length. Then, the reads quality was inspected using FastQC software (Version no. 0.11.8; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) then output statistical result.

Identification and quantification of circRNAs

A total of two algorithms, CIRI2 (v2.0.6; https://sourceforge.net/projects/ciri/files/CIRI2/) and CIRCexplorer2 [V2; (18)], were used to detect circRNAs. Reads were mapped to human reference genome GRCh37/hg19 (http://genome.ucsc.edu/) by Burrows-Wheeler-Aligner (BWA)-MEM [BWA-0.6; (19)] or Tophat (Bowtie 2; http://ccb.jhu.edu/software/tophat/tutorial.shtml), respectively. CIRI2 detects the paired chiastic clipping signals from the mapping information of reads by local alignment with BWA-MEM and combines with systematic filtering steps to remove potential false positives. CIRCexplorer2 uses TopHat and TopHat-Fusion alignment output to detect circRNAs. If a circRNA can be detected by both methods, it will be considered as an identified cirRNA. Back-spliced junction reads identified in CIRI2 were combined and scaled to Reads Per Million mapped reads (bwamem mapping) to quantify every circRNA. Some previous studies used fold change >2, P<0.05 (20,21) and a previous study used fold change >1.5, P<0.05 (22) to screen the differential expression of circRNAs between two samples. In the present study, the criteria were defined as fold change ≥1.5 and P<0.05 to explore as many differentially expressed circRNAs as possible among three groups.

Bioinformatics analysis

Gene Ontology (GO; http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; release 88.0; http://www.genome.jp/) were used to analyze the parent genes to predict circRNA functions. Each differentially expressed circRNAs-targeted miRNAs was predicted with miRanda (August 2010 Release; http://www.microrna.org/microrna/home.do), RNAhybrid (RNAhybrid.2.1; http://bibiserv.techfak.uni-bielefeld.de/download/tools/rnahybrid.html) and TargetScan (Release 7.2; http://www.targetscan.org) software. A total of four software packages [TargetScan (Release 7.2; http://www.targetscan.org), miRDB (5.0; http://www.mirdb.org/), miRTarBase (Release 7.0; http://mirtarbase.mbc.nctu.edu.tw/php/index.php) and miRWalk (miRWalk.2.0; http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/index.html)] were used to predict the target genes of miRNAs in a circRNA-miRNA network. Then, these target genes were subjected to GO and KEGG analysis.

RT-PCR

Total RNA was isolated from 10 samples (fetal lungs) in each group using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, Inc.). The RNA quantity control and concentration were detected using a NanoDrop2000 Spectrophotometer (Thermo Fisher Scientific, Inc.). Subsequently, 1 µg total RNA was isolated as aforementioned and was converted into cDNA using the HiScript® II Q Select RT SuperMix for qPCR (R232-01; Vazyme) according to the manufacturer's protocol [4 µl 5X HiScript II Select qRT SuperMix, 1 µl Random hexamers (50 ng/µl) and 1 µg RNA]. The RT reaction was conducted at 37°C for 15 min and 85°C for 2 min. Next, the RT-PCR reaction was performed using AceQ® qPCR (Q131-01; Vazyme). For RT-PCR, 1 µl cDNA was added to 9 µl master mix, including 5 µl SYBR® Green Master Mix (Low Rox Premixed; Q131-01; Vazyme), 0.2 µl reverse and forward primers, and 3.6 µl diethypyro-carbonate water. Then, the PCR was performed with an ABI 7500 thermal cycler (Applied Biosystems; Thermo Fisher Scientific, Inc.) and the reaction conditions included an initial step at 95°C for 5 min, and 40 cycles at 90°C for 15 sec and at 60°C for 15 sec, 72°C for 1 min and final extension at 72°C for 10 min. The primer sequences used are listed in Table I. PCR was performed in each plate with GAPDH as an endogenous control. All primer sequences were designed and synthesized by Guangzhou RiboBio Co., Ltd. The relative level of each circRNA expression was calculated using the 2−∆∆Cq method (23).

Table I

Primers for reverse transcription-PCR.

Gene name	Primers	Temperature, °C
rno_circ:chr7:24777879-24784993	F: 5′-TTCCAACGCTGAGGACGCT-3′R: 5′-CGTCTTCAATGTCATAGCCGCT-3′	60
rno_circ:chr14:14620910-14624933	F: 5′-GTCGGGTATTGTGCTGCTTG-3′R: 5′-GGTTAAAGTGGGTCTCTGGACA-3′	60
rno_circ:chr3:1988750-1998592	F: 5′-GACAATGCGGGTGCCAATA-3′R: 5′-GCGTCCAAGTGGTTGTTCTCT-3′	60
GAPDH	F: 5′-GAACGGGAAGCTCACTGG-3′R: 5′-GCCTGCTTCACCACCTTCT-3′	60

F, forward; R, reverse.

Statistical analysis

All data are presented as the mean ± SD. All experiments were repeated independently at least three times. The data were analyzed using SPSS 17.0 (SPSS, Inc.) and GraphPad Prism 5.0 (GraphPad Software, Inc.) statistical packages. One-way ANOVA followed by Newman-Keuls was applied to analyze statistical significance among four groups. P<0.05 was considered to indicated a statistically significant difference.

Results

Morphological results of rat lung tissue at E19, E21 and P3

In the E19 group, the respiratory bronchiolitis and the alveolar sacs began to take shape. The interstitial tissue became thinner than before (Fig. 1A). The airway tubes were lined with cuboidal epithelium (Fig. 1D). In the E21 group, the airspaces expanded rapidly and the connective tissue between the airspaces diminished (Fig. 1B), and the epithelium became flat (Fig. 1E). In the P3 group, the terminal alveoli formed rapidly (Fig. 1C) and the secondary septation appeared, which divides alveolar ducts into terminal alveoli. The double capillary network turned into a single capillary system (Fig. 1F).

Figure 1

Morphological observation of fetal lung tissues in three groups. Representative images of lung tissues stained with H&E at (A) E19 (black arrows indicate respiratory bronchiolitis), (B) E21 (black arrows indicate the expanded airspaces) and (C) P3 (black arrows indicate the secondary septation) under an optical microscope. Magnification, ×20. Representative images of lung tissues stained with H&E at (D) E19 (black arrows indicate the airway tubes lined with cuboidal epithelium), (E) E21 (black arrows indicate the flat epithelium) and (F) P3 lung tissues (black arrows indicate the secondary septation) under optical microscope at a higher magnification. Magnification, ×40. H&E, hematoxylin and eosin; E, embryonic day; P, post-natal day.

circRNA expression profile during lung development

The present study examined differentially expressed circRNAs during rat lung development by high-throughput sequencing in the E19, E21 and P3 groups. With CIRI2 and CIRC explorer packages, 9,734, 8,429 and 9,169 overlapped circRNAs were predicted in the E19, E21 and P3 groups, respectively. As illustrated in Fig. 2A and C, 375 of these exhibited differentially expressed circRNAs between E19 and E21 groups (fold change ≥1.5; P<0.05; Table SI), among which 200 were upregulated and the remaining 175 were downregulated in group E21. Meanwhile, as shown in Fig. 2B and D, 358 of these exhibited differentially expressed circRNAs between E21 and P3 (fold change ≥1.5; P<0.05; Table SI), among which 164 were upregulated and 194 were downregulated in group P3.

Figure 2

Differentially expressed circRNAs in lung development. (A) Clustered heatmap exhibited the differentially expressed circRNAs in group E19 compared with group E21. (B) Heatmap exhibited the differentially expressed circRNAs in group E21 compared with group P3. (C) Volcano plots showed the differentially expressed circRNAs in group E19 compared with group E21. (D) Volcano plots showed the differentially expressed circRNAs in group E21 compared with group P3. The red represents the circRNAs with significant difference (fold change ≥1.5; P<0.05) and the blue represents the circRNAs with no significant difference. circRNAs, circular RNAs; E, embryonic day; P, post-natal day.

From all these differentially expressed circRNAs (Table SI), a specific criterion was set for screening circRNAs, which continuously regulate lung development. The circRNAs must meet the criteria that the fold-changes of E19 vs. E21 ≥1.5 and E21 vs. P3 ≥1.5. Then, 1 consistently upregulated and 1 consistently downregulated circRNAs were screened (Table II). In addition, 4 circRNAs presented first upregulation and then downregulation (Table III). In contrast, 1 circRNA was first downregulated and then upregulated (Table III).

Table II

Specific fold-changes of the consistently up- and downregulated circular RNAs in three groups (E19, E21 and P3).

Name	E19 vs. E21		E21 vs. P3		Gene symbol
	log2FoldChange	Regulation	log2FoldChange	Regulation
rno_circ:chr7:24777879-24784993	−2.11547721741994	Down	−1.54432051622381	Down	Polr3b
rno_circ:chr19:24859866-24860320	2.92599941855622	Up	2.05324202127698	Up	Adgre5

E, embryonic day; P, post-natal day.

Table III

Specific fold-changes of first increased and then decreased and first decreased and then increased circular RNAs.

Name	E19 vs. E21			E21 vs. P3			Gene symbol
	P-value	log2FoldChange	Regulation	P-value	log2FoldChange	Regulation
rno_circ:chr14:14620910-14624933	0.003912883	2.3310112	Up	0.032264398	−1.722379793	Down	Fras1
rno_circ:chr1:255636891-255637533	0.018024097	1.685443894	Up	0.009200469	−2.028644787	Down	Btaf1
rno_circ:chr3:1988750-1998592	0.000002489	Positive infinite	Up	0.0485656156	−2.096861539	Down	Ehmt1
rno_circ:chr4:65761402-65770025	0.000015442	Positive infinite	Up	0.014539291	−2.955108543	Down	Atp6v0a4
rno_circ:chr10:16724062-16730140	0.027029049	−2.64385619	Down	0.001996808	3.313660479	Up	Crebrf

E, embryonic day; P, post-natal day.

RT-PCR of the significantly differentially expressed circRNAs

When PCR verification was performed, the specified cycle threshold value of rno_circ:chr19:24859866-24860320, rno_circ:chr4:65761402-65770025, rno_circ:chr1:255636891-255637533 and rno_circ:chr10:16724062-16730140 could not be detected successfully at different time points (data not shown). Therefore, 3 out of the 7 significantly differentially expressed circRNAs were confirmed successfully by RT-PCR in groups E16, E19, E21 and P3. Among them, 2 circRNAs first presented upregulation and then downregulation (rno_circ:chr14:14620910-14624933 and rno_circ:chr3:1988750-1998592; Table III). Another differentially expressed circRNA is rno_circ:chr7:24777879-24784993 (Table II), which showed consistent downregulation.

Analysis by RT-PCR revealed that these three circRNAs exhibited changes in their expression levels, which were consistent with the results of circRNA high-throughput sequencing. The relative expressions of these 3 circRNAs are shown in Fig. 3. It is worth noting that rno_circ:chr7:24777879-24784993 showed a continuously downregulated tendency at the four continuous time points (E16, E19, E21 and P3; Fig. 3C). Meanwhile, the relative expression level of rno_circ:chr3:1988750-1998592 in the E16 group was significantly lower than E19 (Fig. 3B). Similarly, the relative expression level of rno_circ:chr14:14620910-14624933 in the E16 group was significantly lower than E21 as well (Fig. 3A).

Figure 3

Validation of circular RNA high-throughput sequencings results by reverse transcription-PCR. (A) Relative expression of rno_circ:chr14:14620910-14624933 in four groups (E16, E19, E21 and P3) showed a tendency of first upregulation and then downregulation. (B) Relative expression of rno_circ:chr3:1988750-1998592 in four groups (E16, E19, E21 and P3) showed a tendency of first upregulation and then downregulation. (C) Relative expression of rno_circ:chr7:24777879-24784993 in four groups (E16, E19, E21 and P3) showed consistent downregulation. Triplicate assays were performed from each RNA sample. Data are normalized using GAPDH as an endogenous control for RNA input. Error bars represent the mean ± SD. All exhibited statistical significance of reverse transcription-PCR was tested by one-way ANOVA followed by Newman-Keuls. *P<0.05, **P<0.001, ***P<0.0005. E, embryonic day; P, post-natal day.

Bioinformatic analysis of differentially expressed circRNAs

The functions of circRNAs may be associated with their parent genes, as circRNAs share the same parent genes with mRNA and compete with them by linear splicing (24,25). Therefore, parent genes of differentially expressed circRNAs were subjected to KEGG and GO analysis based on their mRNAs to explore the potential functions of circRNAs. GO and KEGG analysis suggested that these differentially expressed circRNAs are involved in biological process, cellular component, molecular function and several biological pathways, such as the cGMP-PKG, Hippo, Wnt, regulating pluripotency of stem cells, TGF-β and PI3K-Akt signaling pathways. Notably, many of these signaling pathways were associated with lung development closely. The top 30 KEGG pathways of differentially expressed circRNAs between E19 and E21 are shown in Fig. 4A. Similarly, the top 30 KEGG pathways which were identified for the parent genes of the upregulated and downregulated circRNAs between P3 and E21 are illustrated in Fig. 4B. In addition, the role of these differentially expressed circRNAs were further investigated by GO analysis. Fig. 4C and D showed the top 10 enriched GO terms in molecular function, cellular component and biological process of differentially expressed circRNAs among three groups.

Figure 4

KEGG pathway analysis and GO analysis of the parent genes of circRNAs. circRNAs with fold change ≥1.5 and P<0.05 were selected from the differentially expressed circRNAs in three groups, and the parent genes of these circRNAs were identified using bioinformatics tools. (A) Top 30 KEGG pathways of differentially expressed circRNAs in E19 compared with E21 using KEGG analysis, according to the number of enriched genes. (B) Top 30 KEGG pathways of differentially expressed circRNAs in E21 compared with P3 using KEGG analysis, according to the number of enriched genes. (C) Top 10 enriched GO terms in the molecular function, cellular component and biological process categories of differentially expressed circRNAs in E19 compared with E21. (D) Top 10 enriched GO terms in the molecular function, cellular component and biological process categories of differentially expressed circRNAs in E21 compared with P3. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; circRNAs, circular RNAs; E, embryonic day; P, post-natal day.

circRNA-miRNA network

The 3 circRNAs verified by RT-PCR were selected to construct a representative circRNA-miRNA network (Fig. 5). The downregulated rno_circ:chr7:24777879-24784993 was predicted to increase the expression level of 25 miRNAs. Rno_circ:chr14:14620910-14624933 and rno_circ:chr3:1988750-1998592 were predicted to downregulate 25 and 23 miRNAs, respectively. In particular, it was identified that rno_circ:chr14:14620910-14624933 can act as a let-7 family sponge. To explore the respective potential functions of the three circRNAs, the target genes of these miRNAs were predicted. In addition, KEGG and GO analysis of miRNAs were conducted to gain insight into each of the three candidate circRNAs. KEGG analysis suggested that the 3 circRNAs involved in some signaling pathways, like Wnt, Hippo, PI3K-Akt, NF-κB and vascular endothelial growth factor (VEGF) signaling pathways. Fig. 6A, C and E showed that the top 10 enriched GO terms of these miRNAs are similar to the analysis of differentially expressed circRNAs in the three categories. Additionally, the top 30 KEGG pathways for the three miRNAs are demonstrated in Fig. 6B, D and F.

Figure 5

Representative circRNA-miRNA network. This network was based on the expression profile results and the related software. The three candidate circRNAs, rno_circ:chr7:24777879-24784993, rno_circ:chr14:14620910-14624933 and rno_circ:chr3:1988750-1998592 (purple nodes) having been validated, were predicted to be functionally connected with their targeted miRNAs (blue nodes) in the network. In particular, rno_circ:chr14:14620910-14624933 can act as a let-7 family sponge, as 'rno-let-7g-3p' is a target of rno_circ:chr14:14620910-14624933. Blue lines indicate targeted miRNAs for rno_circ:chr3:1988750-1998592, green lines indicate targeted miRNAs for rno_circ:chr14:14620910-14624933 and grey lines indicate targeted miRNAs for rno_circ:chr7:24777879-24784993. circRNAs, circular RNAs; miRNA/miR, microRNA.

Figure 6

GO analysis and KEGG pathway analysis of the target genes of miRNAs. miRNAs were presented in the circular RNA-miRNA network, and the target genes of these miRNAs were identified by bioinformatics tools. (A) Top 10 enriched GO terms in the molecular function, cellular component and biological process categories of miRNAs connected to rno_circ:chr3:1988750-1998592. (B) Top 30 KEGG pathways of miRNAs connected to rno_circ:chr3:1988750-1998592 were identified via KEGG analysis according to the number of enriched genes. (C) Top 10 enriched GO terms in the molecular function, cellular component and biological process categories of miRNAs connected to rno_circ:chr7:24777879-24784993. (D) Top 30 KEGG pathways of miRNAs connected to rno_circ:chr7:24777879-24784993 were identified via KEGG analysis according to the number of enriched genes. (E) Top 10 enriched GO terms in the molecular function, cellular component and biological process categories of miRNAs connected to rno_circ:chr14:14620910-14624933. (F) Top 30 KEGG pathways of miRNAs connected to rno_circ:chr14:14620910-14624933 were identified via KEGG analysis according to the number of enriched genes. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; miRNA, microRNA.

Discussion

In the present study, differentially expressed circRNAs were identified in fetal lung development for the first time, to the best of our knowledge. Of these differentially expressed circRNAs, rno_circ:chr7:24777879-24784993 was consistently downregulated among the three groups. Although when PCR verification was performed, the specified cycle threshold value of 4 circRNAs could not be detected successfully at different time points. It was hypothesized that there are two possible reasons for this result. One possibility is that the designed primer might be defective, and the other possibility might be that the content of some circRNA was too low to be detected at specific time points. Following the study by Li et al (26), which examined the function of circRNAs on the regulation of their parental genes during transcription, several previous studies reported that the splicing efficiency of certain genes was negatively correlated with the level of circRNA, suggesting that there are competitive splicing mechanisms between circular transcripts and linear processing (25-27). These findings suggested that specific circRNAs could regulate their parent genes in both a negative and a positive manner at different levels (25-27). In the present study, a significant amount of KEGG and GO terms associated with pulmonary angiogenesis [NF-κB and VEGF signaling pathways (28,29)] and vascular development [Hippo and TGF-β signaling pathways (30-32)] were identified. Lung maturation requires the involvement of various signaling pathways, including the Hippo, Wnt, TGF-β and PI3K-Akt signaling pathways (30-37).

For example, the Hippo signaling pathway plays a critical role in lung development (30). Evidence has demonstrated that the Hippo pathway effector Yap is a key regulator of the formation and branching of epithelial tubes, and mammalian STE20-like protein kinase 1/2 during the differentiation of lung epithelial cells (30,31). At the same time, the TGF-β signaling pathway interacts with the Hippo pathway to regulate epithelial-mesenchymal interactions, distal airway formation and alveolarization (32). In addition, the Wnt signaling pathway, essential during lung development, was also identified. Mice with wnt2 knockout were characterized by decreased cell proliferation and dilated endothelial vasculature (33). Additionally, mice with loss of wnt7 succumbed to respiratory failure due to hypoplastic lungs and abnormal vascular development (34). The Wnt signaling pathway has also been shown to cooperate with the Hippo and PI3K-Akt signaling pathways to regulate epithelial-mesenchymal interactions that are the basis for fetal lung morphogenesis (35-37). In conclusion, the signaling pathways identified in the present study play a vital role in lung development. However, how these circRNAs influence the downstream signaling pathways remains unknown.

It is widely accepted that circRNAs can act as miRNA sponges; however, the potential functions of circRNAs are unclear. As miRNAs regulate a large set of biological processes, the circRNA sponge activity also affects these processes (10). In addition, several previous studies have reported that circRNAs could affect the expression of miRNAs and regulate the signaling pathways in several diseases, including lung cancer (38-40). Therefore, the miRNAs of these differentially expressed circRNAs were predicted. Among them were some specific miRNAs such as cystic fibrosis- [miR-126 (41)], asthma- [miR-221 (42,43)], BPD- [miR-29 (44) and miR-152 (45)] and RDS- linked miRNAs [miR-26a (46)]. It was hypothesized that these differentially expressed circRNAs were closely associated with lung development and diseases. In the circRNA-miRNA network of three candidate circRNAs, it was observed that rno_circ:chr7:24777879-24784993 was the sponge of miR-93. miR-93 was differentially expressed in developing mouse embryos, and was found to promote lung development by regulating VEGF expression and stem cell differentiation (47). Similarly, miR-221, closely linked to rno_circ:chr3:1988750-1998592, was reported to target two Hox genes known to have important functions in embryonic lung branching morphogenesis and epithelial cell fate (43). In particular, rno_circ:chr14:14620910-14624933 can act as a let-7 family sponge. let-7 and its family members, which were highly conserved across species in functions and sequence, were originally discovered in the nematode Caenorhabditis elegans, and regulate cell proliferation and differentiation (48,49). Numerous previous studies have reported that let-7 was expressed in lung tissue and directly regulated RAS expression (49-53). The RAS/mitogen-activated protein kinase signaling pathway is one of the major downstream targets of the FGF signaling pathway (54). The FGF signaling pathway plays a key role in lung development, including lung bud formation, pulmonary branching morphogenesis and lung epithelial cell proliferation (55-57). Therefore, it was hypothesized that these circRNAs were associated with lung development.

To further explore the functions of 3 candidate circRNAs (rno_circ:chr7:24777879-24784993, rno_circ:chr14:14620910-14624933 and rno_circ:chr3:1988750-1998592), circRNA-miRNA network analysis was performed, in addition to GO and KEGG analysis of the miRNA-targeting genes. The Wnt, Hippo, TGF-β and PI3K-Akt signaling pathways were all identified in the KEGG analysis. In addition to these pathways, the VEGF and NF-κB signaling pathways were also included, although some pathways were not included in the top 30 KEGG terms.

In the clinical setting, with the development of neonatology, the survival rate of premature infants is increasing. However, the incidence of lung development-related diseases, such as BPD, RDS and cystic fibrosis, have been gradually increasing as well, severely influencing the prognosis of preterm infants (58). The lungs of premature infants with BPD and RDS are characterized by a reduced alveolar number, thickened septa, malformed pulmonary circulation and lack of pulmonary surfactant, which, in combination with different risk factors, ultimately gives rise to lung injury and lung-related diseases (1,2). Therefore, exploring the physiopathological mechanisms of lung development is crucial. The present study identified certain circRNAs possibly associated with lung development. However, it is necessary for future studies to carry out homologous database analysis in humans and rats to validate these predicted miRNAs, target genes and signaling pathways. Future studies may aim to demonstrate the importance of selected circRNAs in regulating lung development using knocked out rat models. Finally, more experimental and clinical data are required to demonstrate that the selected circRNAs may play a unique, beneficial and essential role in lung development. A challenge will be to further delineate the mechanisms of circRNAs, to further progress the development of therapeutic strategies to selectively block or enhance these mechanisms to effectively prevent or treat lung development-related diseases. We hope to clarify the potential physiopathological mechanism of human lung development, and then these certain circRNAs could become the therapeutic target in the future.

In conclusion, the present study is the first study, to the best of our knowledge, to profile differentially expressed circRNAs at 3 key time points during rat lung development. In total, 7 consistently differentially expressed circRNAs were identified. The circRNA-miRNA interactions were also predicted and circRNA-miRNA networks were constructed for 3 candidate circRNAs (rno_circ:chr14:14620910-14624933, rno_circ:ch r3:1988750 -1998592 a nd rno_circ:chr7:24777879-24784993). These results supported that these novel circRNAs participate in lung development. These findings may also help clarify the physiopathological mechanisms of normal rat lung development, and may further provide a physiopathological basis for lung development-related diseases.