Circular RNAs are abundant and dynamically expressed during the embryonic lung development of C57BL/6 mice.

Circular RNAs (circRNAs), a novel type of endogenous RNAs, can function as microRNA (miRNA) sponges capable of regulating gene transcription, binding to RNA-associated proteins, and even encoding proteins. CircRNAs are involved in various cell behaviors, such as proliferation and apoptosis. The mouse model has also been demonstrated to be similar to that of humans in many studies. To explore the profile of circRNAs during embryonic lung development and their potential functions in lung development-related diseases, mouse embryos at the pseudoglandular phase, canalicular phase, saccular phase, and alveolar phase were collected. High-throughput sequencing was then used to identify a total of 1,735 circRNAs (junction reads ≥5 and p < 0.05). It is well known that the functions of circRNAs are related to host genes. In our study, bioinformatics analysis indicated that the screened host genes were closely associated with lung development and included the Hippo signaling pathway, PI3K-Akt signaling pathways, and TGF-β signaling pathways. Moreover, miRNA sponges are another mechanism involved in lung development. Therefore, we predicted many miRNAs binding to circRNAs, such as miR-17 and miR-20, using the TargetScan and miRanda databases. Previously, miRNAs were proven to be necessary for lung development. The peak expression of circRNAs is distributed at different time points, suggesting their involvement in different stages of embryonic mouse lung development.

Introduction

Circular RNAs (circRNAs) are novel endogenous RNAs that have been highlighted in the field of RNA biology. A recent study reported that abundant, stable, and conserved circRNAs play crucial functions in numerous biologically active processes [1]. A previous study also reported that circRNAs have many complex biological functions, including serving as microRNA (miRNA) sponges, regulating transcription, binding RNA-binding proteins, and encoding proteins [2]. Most importantly, circRNAs were found to be promising molecular biomarkers with a bright future in clinical disease diagnosis and therapy [3].

In recent years, studies have shown that circRNAs demonstrate time specificity and tissue specificity during development [4, 5, 6]. In the porcine model, circRNA-7 (CiRS-7/CDR1as) is upregulated in the cerebellum at embryonic day 115 (E115) and is expressed at its highest level at E60 in other brain tissues. CDR1as is significantly expressed in the nervous tissues, while its content is close to zero in other tissues [7, 8]. Prepulse inhibition owing to impaired neuropsychiatric symptoms was observed in CDR1as−/- mice. MiRNA sponge transcripts, so-called competing endogenous RNA, have been confirmed to suppress miRNA activity, resulting in increased levels of miRNA targets. Previous studies have demonstrated that circRNAs can function as miRNA sponges. For example, CDR1as contains more than 70 selectively conserved miR-7 target sites and was shown to act as a miR-7 sponge regulates post-transcriptional gene expression [9]. The expression of miR-7 was downregulated, while its direct target was upregulated in CDR1as-deficient brains [10]. Another study showed that CDR1as downregulation can result in upregulated miR-7 to further attenuate the activity of ubiquitin protein ligase A, which is involved in clearing toxic amyloid peptide from the brains of patients with Alzheimer's disease [10]. In summary, these studies demonstrating the time specificity and tissue specificity of circRNAs during normal development were of great significance for clinical diseases. However, the association of circRNAs with embryonic lung development has not yet been reported. Therefore, we decided to examine the time-specific and tissue-specific circRNA expression during embryonic lung development and surmised that they play a potential role in lung-related diseases.

Studies have shown that the expression of circRNAs is time and tissue-specific. Few studies have examined the functional mechanism of circRNAs involved in lung development. Mouse embryonic development is similar to human embryonic development, and is one of the most commonly used lung development models [11]. Therefore, we performed circRNA sequencing and analysis of lung tissues in four key stages of mouse embryonic development.

In this study, circRNAs in mouse embryonic lung tissues were detected by RNA-seq and CIRI algorithm, and the functions of these differentially expressed circRNAs during mouse embryonic lung development were further analyzed. After that, a network of circRNA–miRNA was constructed to provide evidence for future research on the molecular mechanism regulating embryonic lung development. Differentially expressed circRNAs at different stages are linked to specific characteristics in the development of embryonic lungs.

Materials and methods

Experimental design

High-throughput sequencing was conducted to assess the differentially expressed circRNAs in embryonic lung tissues divided into four different groups: pseudoglandular phase (E14.5, S1), canalicular phase (E16.5, S2), saccular phase (E18.5, S3), and alveolar phase (postnatal day (PN) 7.5, S4) [n = 2]. Embryonic lung tissues were harvested from C57BL/6 mice at 14.5, 16.5, and 18.5 days of gestation and 7.5 days after birth. For each group, one pair of embryonic lung tissues was randomly selected and retained for histological observation (Figure 1). After the histological evaluation of lung development, we extracted total RNA, removed the ribosomal RNA and linear RNA, and constructed a cDNA library for high-throughput sequencing. Next, we applied two algorithms to detect circRNAs, CIRI and KNIFE. Subsequently, we performed bioinformatics analysis and the differentially expressed circRNAs were further validated. Finally, the circRNA–miRNA network was established using Cytoscape software.

Figure 1

Sample collection and analysis. Embryonic lung tissues were collected from C57BL/6 mice at 14.5, 16.5, 18.5 days of gestation and 7.5 days after birth. Some of them were retained for histological observation, and the remaining lung tissues were placed in liquid nitrogen and preserved at -80 °C. Total RNA was extracted, and a cDNA library was established after the removal of rRNA and the digestion of linear RNA with RNase R. After analyzing these differentially expressed circRNAs, a miRNA-circRNA network was developed.

Sample collection

Experimental protocols were approved by the Nanjing Medical University Animal Care and Use Committee. C57BL/6 mice (20 females and 20 males) from a specific pathogen-free environment were provided by the animal center of the Nanjing Medical University. Pregnant C57BL/6 mice were sacrificed, and whole embryonic lung tissues were immediately harvested from the fetuses on E14.5, 16.5, and 18.5. Newborn mice were sacrificed by decapitation, and the lung tissues were removed on PN7.5d. Based on the time points, the experimental groups were assigned to group S1 (E21), group S2 (E19), group S3 (E16), and group S4 (PN7.5d). Five pregnant mice in each group were sacrificed, and their lung tissues were collected and washed with phosphate-buffered saline (PBS). Two mice in each group were sacrificed for the histological observation of lung tissues. After histological observation, the remaining lung tissues were recovered and RNA was extracted for sequencing. All lung tissues were immediately placed in liquid nitrogen and preserved at -80 °C for follow-up studies.

Histological examination

The four groups of embryonic lung tissues were fixed with 4% paraformaldehyde, embedded in paraffin, and serially sliced at 3 μm. Hematoxylin and eosin (H&E)-stained sections were imaged under a light microscope at a magnification of ×400 to observe embryonic lung development at the four time points.

CircRNA enrichment

The sequencing service (Nanjing Decode Genomics, China) were submitted to the GEO website (GEO accession number: GSE131839). The extraction of total RNA from the lung tissues of neonatal and embryonic mice was performed using TRIzol reagent (Invitrogen, USA). RNA purity was tested using Nanodrop-2000 after the removal of total genomic DNA using DNase I (New England Biolabs Inc., USA). RNA was subject to ribosomal RNA depletion following the protocols of the RiboMinus kit (Life Technology, USA). Next, RNA was fragmented into 200-bp segments using the RNA fragmentation kit (Ambion, USA) and quantified with Nanodrop (Thermo Scientific, USA). The RNA library was loaded on the Illumina HiSeq 2500 for 2 × 150 bp paired-end (PE) sequencing. Initially, basecalls were performed using Solexa pipeline v1.8, followed by alignment of trimmed reads (trimmed 5′,3′-adaptor bases using cutadapt) to the reference genome using Hisat2 (2.0.5) software (genome build mm10). The transcript abundance was estimated with StringTie (1.3.1c). The SRPBM value was calculated with htseq. The circRNA reads were aligned to the reference genome using BWA detecting back-spliced junction reads with CIRI algorithm. The circRNAs in this study were newly discovered and are not easily found in public databases, so we renamed them as “circ + host gene name.”

Bioinformatics analysis

The host gene of circular RNA refers to the region of the gene from which the circular RNA is derived. Host genes include coding and noncoding genes. CircRNAs stem from their host genes and may have regulatory relationships with them. Gene ontology (GO) analyses included biological processes, cellular components, and molecular functions, and they were performed on Gene Ontology (http://geneontology.org/). The significance of the target gene in each GO stage was calculated by the enrichment scoring approach and represented by the p-value. Kyoto Encyclopedia of Genes and Genomes (KEGG) was conducted in corresponding host genes on the KOBAS website (http://kobas.cbi.pku.edu.cn/anno_iden.php). Here, log10 (p-value) reflected the importance of the pathway correlations.

qRT-PCR

The results of the circRNA-seq were quantified by quantitative real-time reverse transcription PCR (qRT-PCR) using the SYBR method. The total RNA from the lung tissues of neonatal and embryonic mice was extracted using Trizol. The reverse transcription mixture containing 500 ng/RNA concentration total-RNA, 2 μl SuperMix (Vazyme Biotech Co., China), and DEPC complemented to a total of 10 μl was prepared. The reaction was performed at 50 °C for 150 min and 85 °C for 2 min, and cDNA was extracted. Next, 2 μl cDNA was added to 10 μl master mix containing 5 μl ChamQ SYBR qPCR Master Mix (Vazyme Biotech Co., China), 2.6 μl DEPC, and 0.2 μl reverse and forward primers. The amplification was carried out using thermal cycler T100 (Bio-Rad, USA). qRT-PCR was performed at 95 °C pre-denaturation for 30 s; 40 cycles of cyclic reaction at 95 °C for 10 s and 60 °C for 30 s; followed by dissolution curve at 95 °C for 15 s, 60 °C for 60 s, and 95 °C for 15 s. The relative levels were calculated using the 2−ΔΔCT method [12]. Data were calculated from three independent records. The primer sequences are depicted in Table 1.

Table 1

Verified circular RNAs and primers in the present study.

RNA names	Primer names	Primer sequences (5′→3′)
chr2:76853605|76856366:ENSMUSG00000051747:exon	circTtn-F	TCCTCTTCCCTATACACCGAAA
	circTtn-R	TAGAAGAACCTCCCCCAACC
chr7:111597803|111599123:ENSMUSG00000038296:intron	circGalnt18-F	TGCGGCACCTCCTCATGT
	circGalnt18-R	ACCCACCTTCACCTCATTGA
chr2:166048455|166051706:ENSMUSG00000027678:exon	circNcoa3-F	CACATCAGCTTTTTGAACATCA
	circNcoa3-R	TTCCCATCCAGTCCTGAGAG
chr15:88811297|88812159:ENSMUSG00000035845:exon	circAlg12-F	CAGCAACCGTCAGCCATATAG
	circAlg12-R	GTCCTCAACATCATGGCTGC
chr16:57569650|57572425:ENSMUSG00000043336:exon	circFilip1l-F	CTGGGACGCAGTTTCTTGAT
	circFilip1l-R	ACCAACAGCCACACCTCTTC
chr17:66809466|66817977:ENSMUSG00000033278:exon	circPtprm-F	GTGCTGCTCATCGTCTCCTT
	circPtprm-R	ACTCAGATGAAGTGCGCTGA
GAPDH	GAPDH-F	GTCTTCACTACCATGGAGAAGG
	GAPDH-R	TCATGGATGACCTTGGCCAG

MiRNA sponge prediction

Currently, circRNAs are receiving attention due to their function as miRNA sponges. TargetScan and MiRanda software were applied to predict the targeting relationship of miRNA-circRNA and identify the binding sites of miRNAs on circRNAs. CircRNAs with more miRNA binding sites were screened as candidate “sponge” circRNAs. Finally, the circRNA–miRNA network was established using Cytoscape 3.6.1 software.

Statistical analyses

Data analysis and mapping were performed using GraphPad Prism7 software. The one-way ANOVA was applied to assess the significant differences between time points (S1 vs. S2, S1 vs. S3, S1 vs. S4, S2 vs. S3, S2 vs. S4, and S3 vs. S4). A p < 0.05 was regarded as significant.

Results

Histopathological assessment

In the S1 group (Figure 2A), the primitive bronchioles and their branches were distributed and extended with a dendritic shape. The epithelium was high columnar (red arrows) with the epithelial cells arranged in rings, a thick interstitium, rare capillaries, and a lack of normal alveolar structure.

Figure 2

A–D Histological morphology of fetal lung tissues by H&E staining, ×400 magnification. Representative images of group S1 (A, E14.5, red arrows indicate the primitive bronchiole epithelial cells), S2 (B, E16.5, the yellow arrow indicates epithelial cells), S3 (C, E18.5, the green arrow indicates secondary septation; the red circle indicates capillaries), and S4 (D, PN7.5, the blue arrow indicates capillaries). H&E, hematoxylin and eosin; E, embryonic day; PN, postnatal day.

In the S2 group (Figure 2B), we observed that the lumen of the bronchi and bronchioles had become larger and the alveolar epithelial cells had changed from high columnar to short columnar square (yellow arrow) with thinner intermedia and more capillaries.

In the S3 group (Figure 2C), the terminal alveoli formed rapidly and secondary septation appeared (green arrow), which divided the alveolar ducts into terminal alveoli.

In the S4 group (Figure 2D), the cavities were extremely dilated with the appearance of mature alveoli around the bronchi, a thin interstitium, a reduced number of cells, and was rich in capillaries (blue arrow).

Differential expression levels of circRNAs

We determined the expression levels of circRNAs in the four groups of embryonic lung tissues using sequencing technology. As a result, 1,735 circRNAs were identified in the four periods (S1 vs. S2, S1 vs. S3, S1 vs. S4, S2 vs. S3, S2 vs. S4 and S3 vs. S4) (Figure 3A). In particular, 520 differentially expressed circRNAs were selected (fold change ≥2, p < 0.05). The volcano plots showed the differentially expressed circRNAs between any two groups (Figure 3C). Among the differentially expressed circRNAs, there were 29 downregulated and 27 upregulated circRNAs in S1 vs. S2, 45 downregulated and 34 upregulated circRNAs in S1 vs. S3, 44 downregulated and 52 upregulated circRNAs in S1 vs. S4, 47 downregulated and 48 upregulated circRNAs in S2 vs. S3, 62 downregulated and 42 upregulated circRNAs in S2 vs. S4, and 39 downregulated and 51 upregulated circRNAs in S3 vs. S4. We selected 13 circRNAs s for further analysis (Figure 3B). The selection criteria were as follows: (a) circRNAs continuously regulate lung development, (b) circRNAs with variation trends, (c) p < 0.05, and (d) fold change ≥5.

Figure 3

A–C Detection of all circRNAs by circRNA-seq technology in four different stages of embryonic development. A. The thermal map revealed that all circRNAs were detected in the E14.5, E16.5, E18.5, and newborn PN7.5 of lung tissues. B. The thermal map revealed 13 differentially expressed circRNAs at four time points. The selection criteria were as follows: (a) circRNAs continuously regulate lung development, (b) circRNAs with variation trends, (c) p < 0.05, and (d) fold change ≥5. circRNAs, circular RNAs; E, embryonic day; PN, postnatal day. C (a–f). Volcano plots of differentially expressed circRNAs between any two groups. The red (downregulated) and green (upregulated) plots represent significantly differentially expressed circRNAs between the two compared samples (fold change ≥2.0, P-value < 0.05).

Functional prediction of circRNAs

Previous studies have shown that most functions of circRNAs are related to host genes [13]. To further clarify the potential functions of the differentially expressed circRNAs, 520 circRNAs (fold change ≥2; p < 0.05) in lung tissues were analyzed. We performed GO analysis and KEGG pathway analysis on the host genes of the circRNAs. The target genes were classified and analyzed according to biological processes, cellular components, molecular function, and KEGG pathways. Figures 4A–D show the KEGG pathway analysis and GO analysis on group S1 vs. S2, S1 vs. S3, S1 vs. S4, S2 vs. S3, S2 vs. S4, and S3 vs. S4, respectively. In the KEGG pathway analysis, circRNAs were identified to be related to the Hippo, PI3K-Akt, Wnt, MAPK, and TGF-β signaling pathways and small cell lung cancer. As a result, bioinformatics analysis indicated that these circRNAs were closely related to lung development.

Figure 4

KEGG analysis and GO analysis of 520 differentially expressed circRNAs. The potential function of circRNA is indirectly inferred by functional enrichment analysis of the circRNAs host gene. The histogram shows the potential biological functions of the differentially expressed genes. The vertical axis and the horizontal axis represent the functional category and correlation, respectively (-log10 P-value). (A) Predicted target genes between the two compared groups identified by KEGG analysis using DAVID online analysis tools. (B) Top 10 enriched GO terms in the biological process (BP) categories of differentially expressed circRNAs between the two compared groups. (C) Top 10 enriched GO terms in the cellular component (CC) categories of differentially expressed circRNAs between the two compared groups. (D) Top 10 enriched GO terms in the molecular function (MF) categories of differentially expressed circRNAs between the two compared groups. KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Verification of the differentially expressed circRNAs

To confirm the reliability of our sequencing results, six out of the 13 differentially expressed circRNAs (circGalnt18, circAlg12, circFilip1, circTtn, circPtprm, and circNcoa3) were verified by qRT-PCR. GAPDH was considered as the internal reference.

The qRT-PCR results verified that the differential expression of circGalnt18, circAlg12, circFilip1l, circPtprm, circNcoa3, and circTtn were consistent with the sequencing results (Figures 5A–F), suggesting the reliability of the circRNA-seq results.

Figure 5

Validation of the circRNA-Seq data. A-F, relative gene expression levels of six circRNAs' (circAlg12, circPtprm, circFilip1l, circTtn, circGalnt18 and circNcoa3). Samples were run in triplicate. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001. The data are presented as mean ± standard deviation. circRNA: circular RNA; qRT- PCR: quantitative real-time PCR.

Prediction of circRNA–miRNA interactions

Many studies have reported that circRNAs can function as miRNA sponges to regulate gene expression [14, 15]. Therefore, this study predicted the downstream miRNA binding sites of circRNAs in the TargetScan and miRanda databases. We observed that 520 differentially expressed circRNAs theoretically have miRNA binding sites. Here, the interactions among the six differentially expressed circRNAs and miRNAs were organized into a network map (Figure 6) and these relationships require further study.

Figure 6

CircRNA–miRNA interaction network showing the interaction between six circRNAs and their potentially combined miRNAs. Red represents differentially expressed circRNAs and yellow represents predictive miRNAs that bind to them.

Discussion

With the recent advances in research methods, increasing numbers of studies have shown that circRNAs can regulate tissue and organ development. A screening study found that 957 circRNAs were expressed abnormally in the lungs of rats with lipopolysaccharide-induced acute respiratory distress syndrome [16]. A recent study described the expression of circRNAs in rat lung development. Alterations in the expression of circRNAs were examined using sequencing technology at three time points during rat fetal development [E19, E21, and N3]. They identified three circRNAs that were characterized by having consistent fold changes (≥1.5) at three time points [17]. It is of significance to uncover circRNA functions in embryonic lung development. It contributes to the prevention, diagnosis, and treatment of abnormal embryos in pregnancy or neonatal lung development-related diseases. Meanwhile, it can provide data for future studies on lung diseases.

The trends of circRNA expression in C57BL/6 mice were examined by circRNA-seq at four developmental stages in our study. The experimental results showed that many circRNAs demonstrated stage-specific differential expression patterns across the four developmental stages. Regarding the validation by qRT-PCR of the results from high-throughput sequencing, we analyzed 13 differentially expressed circRNAs and are only presenting the results that displayed statistically significant differences. There are two possible reasons for the failure of 7 circRNAs to be detected. On the one hand, the design of primers probably be defective; on the other hand, we couldn't detect these circRNAs at specific time points because the expression of some circRNA is low. Of these differentially expressed circRNAs, circGalnt18 was hardly expressed in the pseudoglandular phase (S1) but then gradually increased and reached a peak in the saccular phase (S3), and circNoca3 and circPtprm were almost undetectable in the pseudoglandular phase (S1) but then gradually increased in the canalicular phase (S2), saccular phase (S3), and alveolar phase (S4). Therefore, these altered circRNAs during embryonic lung development were of significance in regulating embryonic lung development. However, there is little research in this field, and further functional investigation into the stage-specific differential expression of these circRNAs is needed.

Previous studies have found that the functions of circRNAs are related to their host genes [18]. Therefore, the target genes were subjected to GO analysis and KEGG pathway analysis. Among the KEGG and GO terms found in this study, many signaling pathways were associated with lung maturation, including the Hippo, Wnt, TGF-β, and PI3K-Akt signaling pathways. Knockout mice for Taz, a Hippo signaling effector, exhibited renal cysts and defects in alveolarization [19, 20]. Furthermore, we found the host genes of circRNAs enriched for the “Hippo signaling pathway” in KEGG pathway analysis are well correlated with the biological functions of embryonic lung development. For example, Bmpr1b exerts two vital roles in early-stage respiratory development. Firstly, it accelerates tracheal formation by suppressing Sox2. Secondly, it restricts the site where the lung bud is initiated [21]. In addition, the “Wnt signaling pathway” was demonstrated to participate in embryonic lung development. In the pseudoglandular stage, the Wnt signaling pathway is involved in the bud formation by the induction of apical constriction. Deficiency of the Wnt signaling pathway is required for air sac formation in the canalicular and saccular stages [22]. Thus, these differentially expressed circRNAs might be closely associated with embryonic lung development.

CircRNAs have a clear function of acting as miRNA sponges [23]. The cCDR1as/miR-7 signal axis has been investigated in many diseases [24, 25, 26]. CDR1as/ciRS-7, a circRNA molecule generated by the transcription of the CDR1 antisense strand, is widely present in the brains of humans and mice, and in zebrafish embryos, both the overexpression of CDR1as and the knockdown of miR-7 resulted in reduced midbrain and impaired midbrain development [1]. In our study, circGalnt18 and its host gene Galnt18 were highly expressed in the lung tissue and circGalnt18 was predicted to have many binding sites of miRNAs. miR-17 is one of the predicted downstream targets of circGalnt18. The miR-17 family has been reported to negatively regulate the fibroblast growth factor (FGF) signaling pathway by acting on Mapk14 (p38-α) and Stat3, which affect the expression of E-Cadherin, and exert a specific function in the branching of mouse epithelial tubules and the maintenance of structural homeostasis. This study found that the downregulation of miR-17 expression led to changes in E-Cadherin distribution and thus interfered with epithelial tubule formation, suggesting that miR-17 is critical for lung development during the pseudoglandular phase [27]. Therefore, circGalnt18 may regulate lung development by absorbing miR-17. In previous studies, miRNAs were proven necessary in lung development. These findings are the basis for further exploration of the function of circRNAs as miRNA sponges, which comprise molecules containing copies of a miRNA binding site that effectively sequester the miRNAs, leading to a loss-of-function phenotype. For example, Zhang et al. studied the role of miR-26 in the formation of pulmonary surfactants and observed that miR-26 overexpression downregulated its target gene SMAD1 and inhibited the production of surfactants in rat type II alveolar epithelial cells [28]. Studies have also reported that miR-126 can enhance the angiogenesis of vascular endothelial growth factor and FGF by inhibiting the expression of Spred-1, a signaling molecule that inhibits angiogenesis [29]. In our study, many binding sites between miRNAs and circRNAs were discovered by software analysis and require further research.

In summary, our study investigated the dynamic changes of circRNAs in mice from the pseudoglandular stage to the alveolar stage. CircRNAs were differentially expressed between the different stages. In addition, the functions of target genes were predicted by bioinformatics analysis, and a candidate circRNA–miRNA network was identified for further study. This research examined the pattern of dynamic circRNA expression during mouse embryonic lung development, and will help with the further exploration of human embryonic lung development. It also provides a novel idea for the diagnosis and prevention of lung development-related diseases during pregnancy and even the neonatal period.

Declarations

Author contribution statement

Meng-Jie Zhang, Jiang-Wen Yin, Jin-Huan Wu: Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Cai-Yun Yuan, Jie Gu: Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Hong-Jun Miao, Zhang-Bin Yu: Conceived and designed the experiments.

Funding statement

This work was supported by the National Natural Science Foundation of China (Grant No. 81870240) and Natural Science Foundation of Jiangsu Province (BK20170147).

Competing interest statement

The authors declare no conflict of interest.

Additional information

Data associated with this study has been deposited at Gene Expression Omnibus under the accession number GSE131839.