During embryonic development, the process of neural tube closure is very sensitive to the external environment, and failure in neural tube closure leads to neural tube malformations^[1]. Incomplete neural tube closure, including brain and spinal cord abnormalities, can result in birth defects called neural tube defects (NTDs)^[2]. NTDs in craniofacial region result in exencephaly, which could lead to embryonic lethality^[3]. Genetics and the environment (e.g., lack of folate or vitamin A) are risk factors of NTDs; however, the mechanisms underlying NTDs in human remain largely unknown^[4]. Some studies stated that fetal neural tube malformed because of the lack of vitamin B12 and folic acid. Vitamin B12 is important for nerve cell metabolism and works with folic acid to promote nervous system development and cell replication. Deficiency of vitamin B12 or folic acid in pregnant women would lead to abnormalities in the development of nervous system, which increases the risk of neural tube abnormalities^[5]. Shanxi is a region with high incidence of birth defects in China. Our research group conducted a large number of investigations and studies in this region several years ago, and confirmed that the high rate of birth defects in this region is related to the lack of one-carbon nutrients locally such as folic acid, and the supplementation of folic acid during pregnancy significantly reduces the rate of birth defects^[6].

A large number of noncoding RNAs are relevant to nervous system development have been detected by large-scale genome sequencing, including microRNAs (miRNAs) and long noncoding RNAs (lncRNAs)^[7]. Unlike linear RNAs have 5′ and 3′ ends, circular RNAs (circRNAs), are endogenous noncoding RNAs with a closed circular structure without 50→30 polarities and poly (A) tails, and resistant to RNA exonucleases. Most circRNAs are highly stable and cytoplasmic, and many of them are tissue-specific, disease-specific, and timing-specificity^[8]. CircRNAs are encoded primarily in exons or introns. There are three main types of cyclization: exonic circRNA, intronic circRNAs, and exon–intron circRNAs, which are formed from back-splicing of exons and/or introns. For cyclization to occur, small nuclear ribonucleoproteins are continuously assembled on the mRNA precursors to catalyze the connection between the downstream 5′ donor sites and the upstream 3′ acceptor sites. Then, the cord tail is inserted to form a circle, and the circRNA is formed by shearing. Cis-acting regulatory elements promote circRNA formation. CircRNAs have abundant miRNA binding sites and act as miRNA sponges to prevent miRNAs interacting with their target genes and regulating gene expression, in what is called a competing endogenous RNA (ceRNA) network^[9]. CircRNA_100290 was shown to function in a ceRNA network to inhibit members of the miR-29 family and regulate the expression of CDK6, which prolonged the cell cycle and proliferation^[10]. In other studies, circRNA MYLK was found to bind competitively to miRNA-29a-3p, leading to increased expression of its target genes DNMT3B, VEGFA, and ITGB1 and the progression of bladder cancer^[11]. CircHIPK3 acted as an miR-558 sponge and suppressed heparanase (HPSE) expression^[12]. In Hirschsprung’s disease, circRNA-ZNF609 was shown to regulate AKT3 by sponging miR-150-5p in bladder cancer cells.

CircRNAs are known to be involved in many of pathological processes and biological functions, including differentiation, proliferation, apoptosis, development, stress response, disease occurrence, and development^[13]. Importantly, some circRNAs are tissue-specific and have stable structures because they do not have accessible terminals and are highly resistant to exonucleases^[14]. CircRNAs have been found to have conserved sequences and expression patterns, especially in brain^[15]. In mouse and human brain tissues, a circRNA, Cdr1as, that binds specifically to miR-7 and mir-671 was detected. The CRISPR-cas9 genome editing technique was used to selectively eliminate Cdr1as in mice, and the results showed that the expression of most miRNAs was changed in the four studied brain regions^[16]. Many circRNAs are highly brain tissue-specific (e.g., specific circRNAs have been found in excitatory neurons^[17], and 80% of the circRNAs that were highly expressed in mouse neuron were also detected in human brain. Many of the genes that encode circRNAs are specifically expressed in the brain. When host genes are expressed in other tissues, the proportion of transcripts that encode circRNAs in the brain increases remarkably, suggesting neuro-specific regulation of circRNA production. Genes that carry brain-specific circRNAs, including CircPHF21A, CircirulP4, Circelf2, Circrims2, and circ-MyST4, have key regulatory functions during neuronal and brain development^[18]. These results suggest that the circRNA host genes are involved in biological processes relevant to synaptic functions (such as nervous system development, neurogenesis, and neuronal differentiation), including synapses, presynaptic active regions, presynaptic membranes, and postsynaptic density^[19]. CircRNA expression profiles were found to change significantly in mouse neural tissues after transient focal ischemia^[20].

Environmental factors such as various teratogenic agents, maternal obesity and diabetes, increased oxidative stress may influence neural tube closure (NTC) through a direct effect on metabolic regulation in the context of embryonic development^[4]. Deficiency of folates at the cellular level may be responsible in NTDs as a result of disturbed bioavailability of folates. In recent years, a series of RNA-sequencing analyses revealed that circRNAs dynamically expressed in various developmental stages and physiological conditions. It implicated the role of circRNAs in stress response and explored potential mechanisms underlying their regulations. However, the effects of folate metabolism disorder in the role of circRNAs in neural tube defects (NTDs) remains unclear. In this study, we statistically compared the circRNA expression profiles between NTDs and normal brain tissues in low folate diet and MTX induced NTDs mouse model by high‐throughput sequencing. We identified a small number of circRNAs that are up- or down-regulated. These circRNAs may exhibit an important function in NTDs and neural tube development. The results were confirmed in human NTDs fetuses with their maternal blood serum under low folate category, indicating that abnormal expression of Circzfp644 and Gas7 may disrupt neural development and contribute remarkably to NTDs risk.

C57BL/6 mice were purchased from JinMuYang Ltd aged from 8–10 weeks old, SPF grade. After adaptive feeding for 2 days, mice were fed with low folate diet for 4 weeks. Female mice and male mice were caged overnight at a ratio of 2:1, and the vaginal plug was found in the morning of the second day, which was defined as 0.5 days of gestation. NTDs group (case group) were induced by intraperitoneal injection of methotrexate (MTX) (Sigma, USA) at 1.5 mg/kg body weight on E7.5. The control (con group) group were fed with a normal diet. On E13.5 Embryos, Con (n = 3) and NTD (n = 3) embryos were dissected and placed in ice- cold tubes, DEPC-treated PBS, respectively. Brain vesicles and spinal cord tissue were isolated from E13.5 d mice, and the samples were divided into three tubes. NTDs Mouse Model were performed as previously mentioned. The details were performed as previously described^[21]. All animal experiments were conducted in compliance with the guidelines of the Institute for Laboratory Animal Research, Capital Institute of Pediatrics.

Cultured cells and homogenated animal tissues were prepared and total RNA was isolated using TRIzol Reagent (Invitrogen, USA) according to the manufacturer’s protocol. RNA was checked directly after extraction by Bioanalyzer 2200 (Aligent) and stored at –80 ℃ until machine analysis. The RNA with integrity number (RIN) > 6.0 is right for rRNA depletion. The RNAs were qualified with RIN > 6.0. The miRNAs were purified by miRNeasy Mini Kit (Qiagen) and agarose gel electrophoresis was used in the determination.

Pooled RNA samples were set for the cDNA libraries with the NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina by following the manufacturer’s instructions. Samples were fragmented into 150–200 bp after rRNA depletion. The MPure XP Beads (beckman coulter) was used to harvest the target bands. The final cDNA libraries were completed through product purification and enrichment, and then the Illumina HiSeqXTen was used to sequence 150 bp paired-end in a single lane.

MiRNA library was prepared with total RNA from mice brain. Sample isolation was completed by following to the user’s protocol. RNAs were extracted and quantified after purification. The HiSeq 2500 (illumina, USA) was used for library sequenced.

Clean reads were picked from the raw reads while the adaptor sequences were deleted, the clean reads were then adjusted to mouse genome (version: GRCh38 NCBI) with hisat2. The transcripts were collected to analyze their differential expression in TopHat and Cufflinks^[22].

CircRNA interacts with miRNA like a miRNA sponge and regulate posttranscriptional genes through binding with miRNA recognition elements. We proposed a gene co-expression network to detect gene crosstalk and construct a gene co-expression network based on the normalized signal intensity of specific expressed genes. The pearson correlation was used in every single pair of genes and the most significant correlation pairs were chosen for the network. Degree centrality is significant within a network to measure centrality of a gene. It determines the relative importance. Furthermore, w-cores in graph theory were used instead of graph topology analysis to study certain properties of the networks. The network structure analysis is to determine the relationship between core regulatory factors and its most adjacent genes. They own the biggest differences. Gene ontology illuminated the cellular component, molecular function, and biological process of mRNAs. (GO; http://www.geneontology.org). Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used for pathway enrichment analysis, annotation and visualization of target genes function.

Mouse teratocarcinoma cells (F9) were purchased from National Medical Centre of Cell Line Resource. F9 were fed in DMEM with 10% fetal bovine serum (Gibco) and cultured in an incubator containing 5% CO₂ and 37 °C humidified airs.

Proteins were separated in 10% SDS-PAGE. The isolated proteins were transferred from the gel to NC membrane and blocked with 5% skim milk. Incubation with primary antibodies Anti-Gas7 antibody (1:1,000 dilution; Santa Cruz), Anti-GAPDH antibody (1:1,000 dilution; Santa Cruz) overnight at 4 °C, and then incubated with secondary antibody (goat anti-rabbit IgG, 1:2,000, CST, U.S.A.) for 40 min at room temperature.

The siRNA targeting against circzfp644 and negative control, as well as mimics of miR-20b-5p were designed and synthesized by GenePharma (Shanghai, China). F9 cells were transfected using Lipofectamine 2000 (Invitrogen) according to the manufacturer’s instruction. The designed sequences were shown in Supplementary Table S1 (available in www.besjournal.com).

F9 cells were cultivated in 24 well plates for transfection. Three UTRs of circrnazfp644 and Gas7 were amplified by PCR, and then subcloned into luciferase reporter vector GP-miRGLO basic. GenePharma (Shanghai, China) Renilla luciferase was used as the control reporter. miR-20-5p mimics or negative controls (100 nm) were transfected with Lipofectamine 2000 (Invitrogen, California, USA) flowing to user's instructions. After 36 hours, the cells were lysed and luciferase activity was detected. The values were normalized to Renilla luciferase transfection control using the dual luciferase reporter gene detection kit (Promega) according to the manufacturer's instructions. The sequences were listed in Supporting Information Supplementary Table S1.

Total RNA was extracted using the Trizol reagent (Invitrogen, USA). The RNA to cDNA reverse transcription was performed using Revert Aid First Strand cDNA Synthesis Kit (ABM, Canada). Real-Time PCR was conducted in Quant StudioTM 7 Flex PCR System (Thermo). miRNA was conducted in Hairpin-itTM microRNA and U6 snRNA Normalization RT-PCR Quantitation Kit. The threshold cycle of each sample was recorded, and data were analyzed by normalization to GAPDH values using the 2^−ΔΔCt method. U6 was used as an endogenous control gene for miR-20-5p. The primer sequences were listed in Supplementary Table S1.

Clinical samples used in this study included 8 pairs of human NTDs and matched normal tissues were collected from Lvliang area in Shanxi Province with informed consent accepted by patients and their families. The International Classification of Diseases (ICD-10) was used to classify NTDs. Participants were diagnosed by local clinicians. All participates gave written informed consent. The details of surgical procedures are described in other studies^[23]. The epidemiological studies are described in detail in our publication^[24]. The Ethics Board of Capital Institute of Pediatrics approved the study protocol^[6,25].

In this study, NanoString nCounter was used to detect the specific number of transcripts in human NTD fetuses brain tissues. This technology was performed as previously described^[26].

Data were analyzed using GraphPad Prism (Version 9.2.0) and SPSS version 19.0 (SPSS Comp., Chicago, IL, USA). All results have been repeated for more than three times. Data were presented as means ± SD. Student’s t test or ANOVA analysis was used for group comparisons. When the data did not conform to a normal distribution, nonparametric tests were used. P < 0.05 was considered statistically significant and is presented as *P < 0.05 or **P < 0.01.

We used mice that were kept under a low folate diet for 4 weeks to construct a mouse model of NTDs by intraperitoneal injection of methotrexate (MTX) (1.5 mg/kg) at embryonic day 7.5 as reported previously^[26]. CircRNA transcripts were obtained by high‐throughput sequencing analysis of three pairs of normal and NTDs mouse brain tissues. Hierarchical clustering was performed to visualize the circRNA expression patterns in normal and NTDs mouse brain samples at embryonic day 13.5 (Figure 1A), and the volcano plots showed the variations in circRNA expression levels between the normal and NTDs brain samples (Figure 1B). We found 57 circRNAs that were differentially expressed in the NTDs brain samples compared with their expression in the normal samples; 16 of them were upregulated and 41 were downregulated (Supplementary Table S2, available in www.besjournal.com). To identify important circRNAs associated with NTDs, we constructed gene co-expression networks and predicted the roles of the target genes on the differentially expressed circRNAs by GO-Biology process (BP) analysis (Figure 1C, Supplementary Table S3, available in www.besjournal.com).

Figure 1.  Differentially expressed circRNAs in nervous tissues. (A) Hierarchical cluster analysis of differentially expressed circRNAs in NTDs and normal brain samples; each group includes 3repeats (greater than 1.5-fold changes in expression; P < 0.05). Expressions represented in different colors, indicating expression above and below the median across all samples. (B) Volcano plots were constructed in fold-change and p-values. The vertical lines mean 1.5-fold up- or down-regulation in normal and NTDs samples, and the horizontal lines correspond to P-value. The red nodes in the plot represent the differentially expressed circRNAs showing statistical significance. (C) Biology Process (BP) analysis for the host gene dysregulated genes (NTDs vs. Con).

Hierarchical clustering was performed to visualize mRNA expression patterns in normal and NTD brain samples (Figure 2A, Supplementary Table S4, available in www.besjournal.com), and the volcano plots showed the variation in mRNA expression levels between the normal and NTDs brain samples (Figure 2B). We performed Gene ontology (GO) and KEGG pathway enrichment analyses to predict the functions of the target genes (mRNAs). The differentially expressed genes were enriched for multiple GO terms under the three main GO categories, biological process, molecular function, and cellular component. Among them, signaling pathways were highly enriched (Supporting Information Supplementary Figure S1, Supplementary Table S5, available in www.besjournal.com). The enriched KEGG pathways included PI3K-Akt signaling pathway, cGMP-PKG signaling pathway, and protein digestion and absorption. (Figure 2C, Supplementary Table S6, available in www.besjournal.com).

Figure 2.  Differential expression and enrichment analysis of mRNA in Neural Tube Defects . (A) cluster analysis of mRNA between NTDs and normal brain samples; each group contains three repeats (greater than 1.5-fold difference in expression; P < 0.05). Results are represented in different colors, indicating expression levels above and below the median across all samples. (B) Volcano plots displayed in fold-change and P-values. The vertical lines mean 1.5-fold up- or down-regulation between normal and NTDs samples. and the horizontal line correspond to a P-value. The red-colored point in the plot represents the differentially expressed mRNA with statistical significance. (C) KEGG enrichment of biological process molecular function and cell component of mRNA. The –log₁₀ P value of enrichment is shown on x axis; the numbers represent the number of associated proteins for each term.

Figure S1.  Gene Ontology (GO) enrichment of biological process molecular function and cell component of mRNA. The -log₁₀ P value of enrichment is shown on X axis; the numbers represent the number of associated proteins for each term.

Hierarchical clustering was performed to visualize miRNA expression patterns in normal and NTD brain samples (Figure 3A, Supplementary Table S7, available in www.besjournal.com), and the volcano plots showed the variation in miRNA expression levels between the normal and NTDs brain samples (Figure 3B). We found 16 miRNAs that were differentially expressed in the NTDs brain samples compared with their expression in the normal samples; 7 of them were upregulated and 9 were downregulated. Overall, miRNA-20-5p had a statistically significant increase in NTDs mouse sample. Consistent with RNA-seq data, RT-qPCR showed that mmu-miR-20b-5p was significantly upregulated in brain and spine of NTDs mouse tissue samples (Figure 3C). Next, we performed ceRNA analysis. On previously the basis of the analysis strategy, the differentially expressed miRNAs were screened and collected from miRanda (Score ≥ 140 and Energy < −15). The negative correlation between miRNA and circRNA, miRNA and mRNA, miRNA and circRNA was selected^[27]. Then 50 circRNA-miRNA- mRNA pathways were constructed, including 5 circRNAs, 12 mRNAs, and 4 miRNAs (Figure 3D). The circRNA-miRNA-mRNA relationships were shown in Supplementary Table S8, available in www.besjournal.

Figure 3.  CeRNA networks analysis in the NTD mouse brain tissue. (A) Cluster analysis of small RNA between NTDs and normal brain samples; each group contains 3 individuals (greater than 1.5-fold difference in expression; P < 0.05). (B) Volcano plots were constructed using fold-change values and P-values. The vertical lines correspond to 1.5-fold up- and down-regulation between normal and NTDs samples. (C) miR-20-5p in cranial neural tissue of mouse NTDs was measured by RT-qPCR. Data are mean ± SD (n = 3). (D) circRNAs-miRNAs-mRNAs network consists of 12 genes. In the network, circle represents circRNA and rhombus represents mRNA and inverted triangle represents miRNA, and their relationship was represented by line.

The RNA-seq showed that miRNA-20-5p had a statistically significant increase in NTDs sample and RT-qPCR analysis acted consistently with the RNA-seq data. Above our ceRNA analysis, we found that circzfp644 - miRNA-20-5p target on four different mRNAs, Cngb3, Nfic, Chmp4c, Scara5, Cygb and Gas7. Gas7 is expressed primarily in terminally differentiated brain cells and predominantly in mature cerebellar Purkinje neurons. GAS7 plays a putative role in neuronal development in human. Gas7 is involved not only in neurogenesis but also in the regulation of neuronal cell death in murine. Gas7 is preferentially expressed in the nervous system and participates in the neurogenesis of mammals^[28]. To demonstrate the feasibility and value of the constructed ceRNA networks, we chose the circzfp644 - miRNA-20-5p - Gas7 network for further validation. Firstly, we examined the expression level of Gas7 and circzfp644 in the NTDs tissue samples. Gas7 was in one of the top 20 differentially expressed mRNAs (Supporting Information Supplementary Table S4), and found to be significantly downregulated in the NTDs tissue samples by RNA-seq (Figure 2B). Consistent with RNA-seq data, RT-qPCR assays were performed on cranial neural tissue of E13.5 mouse embryos, indicated that the mRNA level of Gas7 decreased significantly, in NTDs embryos compared to that in controls (Figure 4A). Western blots showed that the Gas7 level was remarkably decreased in cranial neural tissues in E13.5 embryos (Figure 4B). Next, as shown in Figure 4C, Circzfp644 expression was significantly downregulated, which is consistent with the results from the high-throughput sequencing data. We also analyzed that the exon structure of Circzfp644, which is derived from zinc finger protein 644 (zfp644) gene, was 3,539-nt long. The same results found as previously studies^[29]. The back-spliced junction of Circzfp644 was amplified by RT-PCR using divergent primers (Supporting Information Supplementary Figure S2, available in www.besjournal.com). To further detect the role of circRNA-zfp644, small interfering RNAs (siRNAs) were constructed and transfected into in F9 cell lines. Both of RT-qPCR and western blots showed that Circzfp644 knockdown significantly reduced the Gas7 levels in F9 cells (Figure 4D and 4E). Together, these results suggest that Gas7 is the main candidate target of Circzfp644 in the ceRNA networks.

Figure 4.  Circzfp644 targets Gas7 in the NTD mouse brain tissue. (A) Gas7 mRNA in cranial neural tissue of mouse NTDs was quantified by RT-qPCR. Data are mean ± SD (n = 3). (B) Cranial neural tissue of normal and NTDs mouse was achieved at E13.5, and analyzed by western blotting. Aliquots of total lysates were immunoblotted to Gas7 antibody. (C) Expression of circRNAzfp644 was quantified by qPCR. ΔCt was used to measure gene expression normalized against GAPDH expression levels. Normal group = 3, NTDs group = 3. The presented values are the means ± SD. (D, E) mRNA and protein levels of Gas7 in the F9 cells with circzfp644 knockdown.

CircRNAs function as miRNA sponges to regulate the expression of the target genes^[30]. Next, we explored whether Circzfp644 acted as a sponge to regulate Gas7 expression through the ceRNA networks. The 3′ UTR sequence of Circzfp644 complemented the sequences of mmu-miR-20b-5p. We selected miRNAs that were common for Circzfp644 and Gas7 based on the bioinformatics analysis and preformed AGO2 immunoprecipitation assays to determine whether Circzfp644 was an intermediate for AGO2 and miR-20-5p. As shown in Figure 5A, Circzfp644 was specifically enriched in miR-20-5p mimic cells. To determine whether Circzfp644 and Gas7 could be regulated by miR-20-5p, we constructed luciferase reporters containing wild-type and mutated binding sites of Circzfp644 and Gas7. The luciferase reporter assays showed that wild-type Circzfp644 and Gas7 were significantly reduced in miR-20-5p mimics, whereas the mutated Circzfp644 and Gas7 were not affected compared with the controls (Figure 5B, C). The Q-PCR analyses showed that miR-20-5p mimics significantly decreased both Gas7 and circRNA-zfp644 levels in F9 cells (Figure 5D, E). The western blot analyses showed that miR-20-5p mimics significantly decreased Gas7 protein levels in F9 cells (Figure 5F). These results suggest that Circzfp644 acts as a sponge for miR-20-5p to regulate Gas7 expression via a ceRNA network.

Figure 5.  Circzfp644 sponges with miR-20-5p and inhibits its activity for Gas7 expression. (A) Enrichment of circzfp644 detected through RIP in F9 cells transfected with miR-20-5p mimics or mimics NC. (B) Luciferase reporter assay of circzfp644 and Gas7-3′UTR in F9 cells co-transfected with miR-20-5p mimics or mimics NC. (C) Luciferase reporter assay of Gas7–3′UTR with circzfp644 knockdown in F9 cells. (D) circzfp644 was measured by RT-qPCR. Data are mean ± SD circzfp644expression in F9 cells transfected with miR-20-5p mimics.(E) Gas7 mRNA was measured by RT-qPCR. Data are mean ± SD Gas7 expression in F9 cells transfected with miR-20-5p mimics. (F) Cell lysate was detected Gas7 antibody. F9 cells were transfected with the miR-20-5p mimics or mimics NC.

CircRNAs have been shown to be highly enriched in neural tissues of mammals^[31-33]. We evaluated brain tissue of human NTDs fetuses and the folate levels in the corresponding maternal serum. Eight fetus brain and their folate level paired, including fetuses with spina bifida, with closely matched age and sex of the fetuses were selected, as summarized in Table 1. Maternal serum folate levels were lower in the pairs with fetuses with spina bifida than in the controls (Figure 6A, Supporting Information Supplementary Figure S3A, available in www.besjournal.com). Furthermore, we quantified Cirnzfp644, Gas7 and miRNA-20-5p expression using NanoString nCounter in brain tissue of fetuses with 8 NTDs patients. We test of normality tests’ between the NTD and control groups. However, the results of normality test do not conform to normal distribution. nCounter analysis showed that expression of Cirnzfp644 and Gas7 expression levels were downregulated in human NTDs fetuses brain samples, whereas miRNA-20-5p expression was upregulated with nonparametric test. (Figure 6B, Supporting Information Supplementary Figure S3B). Western blot analysis of the eight NTDs samples and eight controls with age and sex matched showed that Gas7 levels in most anencephaly samples decreased by 10%–80% (Figure 6C). These results indicated that low folate levels may decrease both Gas7 and Circzfp644 levels, thereby potentially increasing miRNA-20-5p expression levels, leading to NTDs.

Figure 6.  Down-regulated of Circzfp644 and Gas7 expression in low-folate NTDs. (A) Detection of maternal serum folate between normal fetus and NTDs fetus. Data displayed in IQ^R. (n = 8), ^*P < 0.05, by Mann-Whitney U test. (B) The mRNA expression of Circzfp644, GAS7, MiR-20-5p in NTDs fetuses’ brain, detected by Nanostring. Data displayed in IQ^R. (n = 8), ^*P < 0.05, by Mann-Whitney U test. (C) Detection of Gas7 in brain tissues from normal fetus and NTDs in low folate maternal serum. GAPDH was used as a loading control. Data are mean ± SD. (n = 3), ^*P < 0.05, by Student’s t-test.

Figure S3.  (A) Statistical analysis of f maternal serum folate concentrations between normal fetus and NTDs fetus. (B) Statistical analysis of Cirnzfp644, Gas7 and miRNA-20-5p expression between human con and NTDs tissues.

NTDs are congenital malformations mainly led by abnormal development of brain and/or spinal cord due to incomplete closure or disorders of segmentation of the neural tube during early embryonic development. Many studies on NTDs have focused on miRNAs and lncRNAs, whereas the role of circRNAs in the occurrence of NTDs is still largely unknown. In this study, whole transcriptome data were systematically analyzed to obtain the expression profiles of circRNAs, miRNAs, and mRNAs in embryonic day 13.5 NTDs mouse brain tissues. On the basis of the co-expression profiling analyses of circRNAs and mRNAs, we identified Circzfp644 as a significantly downregulated circRNA in NTDs tissues. Circzfp644 exerted its function through a ceRNA network by competitively binding to miR-20-5p, which had a suppressive effect on the target gene Gas7, implying that Circzfp644 had potential role in the development of NTDs via a ceRNA network (Supplementary Figure S4, available in www.besjournal.com).

Figure S4.  The flowchart delineates the steps for identifying and validating circRNAs in brain tissue of mouse NTDs.

Exon transcripts can form circRNAs by a back-splicing process. CircRNAs are widely expressed in human cells, and their expression levels can be ≥ 10 times more than their linear isomers^[13]. Compared with other noncoding RNAs, such as miRNA and lncRNAs, these characteristic circRNAs are naturally ideal biomarkers and potential therapeutic agents. CircRNAs are highly expressed in human brain, mouse brain, and rat brain. Examples of brain-specific circRNAs include Circrims2, Circtulp4, Circelf2, Circphf21a, and Circ-myst4, most of which are encoded in genes with key regulatory functions during neuronal and brain development^[32]. In this study, we established a low folate diet with MTX-induced mouse NTDs model at E13.5 and found 57 circRNAs that were differentially expressed in the NTDs brain samples compared with the normal samples; 16 of them were upregulated and 41 were downregulated (Figure 1A, B). These dysregulated circRNA may provide new insights into the role of CircRNA in neural development and NTDs occurrence. To our knowledge, this is the first reported study of circRNAs expression profile in NTDs. Our study also found the that low levels of circzfp644 observed in both brain tissue of mouse and human NTD cases in the corresponding NTDs samples with low-folate levels (Figure 4C and Figure 6B). Previously studies found that the circRNAs involved in various stress such as redox homeostasis and metabolism^[34]. CircRNAs may function to maintain homeostasis in different stress response to adapt response pathways. Many circRNAs is regulated as a product of increased transcription. circRNAs may play important roles response to stress, such as function on transcription of stress related genes, protecting linear RNAs. Next, we also found that the differentially expressed mRNAs were enriched in KEGG pathways associated with nervous system development (Figure 2C). We explored the downregulation of Gas7 in brain tissues of fetuses with NTDs that were paired with the corresponding maternal serum folate levels (Figure 6C). These results suggest that folate deficiency may have contributed to the onset of NTDs by altering Circzfp644 levels, thereby affecting the expression of Gas7. This sequential regulation of Circzfp644 and Gas7 may be a potential risk factor for NTDs in response to folate deficiency.

The ceRNA hypothesis states that RNA transcripts such as mRNAs, lncRNAs, and circRNAs share miRNA response elements, compete to bind to miRNAs, and form complex post-transcriptional regulatory networks^[35]. We performed bioinformatic analyses to select miRNAs that shared common binding sites with Circzfp644 and Gas7. Gas7 is expressed preferentially in the brain and is involved in morphological differentiation and neurogenesis^[36]. In this study, we showed that Circzfp644 upregulated Gas7 expression (Figure 4D). Then, we constructed Circzfp644 luciferase reporters and found that miR-20-5p reduced the luciferase activity of the Circzfp644 luciferase reporter (Figure 5B). These results showed that the strongest binding was between Circzfp644 and miR-20-5p, which confirmed that miR-20-5p was the binding target of Circzfp644 (Figure 5B and 5D). We also showed that Gas7 mRNA was the direct target of miR-20-5p (Figure 5C, 5E and 5F). Circzfp644 has been validated as a sponge of miR-93-5p, and was shown to function in a ceRNA network to upregulate Limk1 expression by sequestering miR-93-5p in angiotensin II-treated cardiomyocytes^[37]. CircZFP644 was found to be negatively correlated with miR-21-3p in Severe acute pancreatitis (SAP)^[29]. Together, these results indicate that miR-20-5p can bind with both Circzfp644 and Gas7.

Many studies have shown that overexpression of a miRNA in different cell environments can activate nerve cell differentiation. For example, miR-124 may be the best example of a miRNA that controls the fate of nerve cells^[38], and miR-20 was found to play a significant role in post-transcriptional regulation of cell proliferation and induced cardiomyocyte hypertrophy^[39]. However, the roles of miR-20-5p in NTDs remain unclear. The abundance of miR-20-5p was higher both in the mouse and fetal NTD brain tissues (Figure 3C and Figure 6B). The RNA immunoprecipitation assays confirmed the binding of Circzfp644 to miR-20-5p (Figure 5A), resulting in less miR-20-5p to target Gas7. Previous findings have shown that miR-362-5p promotes cell proliferation and cell cycle progression by targeting GAS7 in acute myeloid leukemia^[40]. The identification of a Circzfp644–miR-20-5p–Gas7 ceRNA axis in this study expands the understanding of the underlying mechanism of the progression of NTDs. A number of studies have demonstrated that ncRNAs are important regulators of human brain function. In brain development, some specific miRNAs function on neurons, cerebral cortex, and dendritic branching. Some studies showed that circRNA is differentially expressed in different brain regions, which suggest its potential role in brain development^[19]. However, due to the more complexity of human brain function, our research will investigate the interaction between miRNA and circRNA, and the molecular mechanism of transcription factors regulation and chromatin-modifying enzymes in future.

In this study, we found that dysregulated circRNAs identified candidate circRNAs associated with NTDs. Furthermore, we have shown that Circzfp644 competitively binds to miR-20-5p to inhibit its suppressive effect on Gas7 expression. These findings link the abnormal expression of Circzfp644 and Gas7 to the development of NTDs under low-folate conditions, and contribute to the understanding of nervous system development and indicate a potential approach for the diagnosis and treatment of NTDs.

SW conceived and designed this study. SW, YZ, and PP participated in laboratory work. SW, YZ, FL, and XH performed the data analysis. SW and YW participated in the writing of the manuscript. SW and TZ participated in advising and revising the manuscript critically. All authors read and approved the final manuscript.

The study was conducted according to the guidelines of the management and utility of experimental animals, and approved by the Ethics Committee of Capital Institute of Pediatrics, Beijing and all animal manipulations were strictly performed following the relevant laws of China. (SHERLLM2021024, DWLL2021011).