Introduction

An increasing number of experimental studies have shown striking similarities between early embryonic cells and tumor cells in gene and protein expression, metabolic pathways, and important biological processes of ontogeny1,2,3. Embryonic development is the process in which cells continue to divide, proliferate, and differentiate from fertilized eggs, gradually forming tissues and organs, and finally forming a complete and mature individual. All cell groups and organs has gradually formed through various stages of early embryonic development. Cancer is a disease which cells grow uncontrollably, characterized by infinite proliferation, transmutation, and potential metastasis. Moreover, these cells can destroy normal cells. In 1892, French biologists Lobstein and Recamier hypothesized the embryonic origin of tumors. In the 1970s, Pierce proposed that tumorigenesis is largely related to developmental biology4. Aggressive tumor cells and pluripotent embryonic cells are similar in their ability on self-renew, proliferation, and plasticity. In 2019, Coorens et.al found that childhood kidney cancer has embryonal precursors5. In 2020, Coorens et.al proved that the normal genomic structure of the placenta is different from that of other human organs, and more similar to tumors, and the genetic mutation site is the same as that of many childhood cancers6. In 2020, a single-cell sequencing analysis of approximately 212,000 representative cells from human fetuses, hepatocellular carcinoma (HCC), and mouse livers by the Singapore Genetics Institute showed that there is a shared fetal ecosystem between fetal liver and liver cancer, suggesting that the tumour microenvironment may undergo regression remodelling7. In addition, at the miRNA level, tumor occurrence and embryonic development also have the same mechanism8,9.

The neuroectoderm is the central ectoderm on the dorsal surface of early vertebrate embryos. Under induction of the notochord and paraxial mesoderm, the ectoderm thickens to form neural plate. The neural plate widens and elongates with the lengthening of the notochord. Concomitantly, the central depression forms the neural groove, and the two sides of the uplift form a neural fold. On the 22nd day of human embryonic development, the folds begin to close and form a neural tube. In mouse, this period is from embryonic day (E) 8.5 to 10.5. Failure to close the neural tube leads to severe malformations such as spina bifida, and anencephaly, collectively known as NTDs10. The causes of NTDs are numerous and highly complex, including genetic and environmental factors, including exposure to high sugar, high fever or multiple toxic substances. However, the pathogenesis of NTDs caused by complex genetic and environmental factors has not been fully elucidated11.

In addition, numerous studies have shown that various types of tumors originate from normal tissue. Tumors that develop from the ectoderm of early embryos are often referred to as primitive neuroectodermal tumors, including glioblastoma (GBM). GBM is classified as a grade IV (the most severe) astrocytoma and develops from a lineage of astroglial cells (astrocytes) that support neural cells. The molecular mechanism of GBM remains unclear, and more reliable and effective treatment are needed12.

Our previous research demonstrated that the expressions of some genes may be involved in NTDs formation and GBM tumorigenesis through miRNA regulation13. This study try to evaluate the key genes and underlying mechanism involved in NTDs formation and GBM tumorigenesis to identify small-molecule drugs with possible clinical value.

Materials and methods

Mouse mRNA sequencing data

Retinoic acid (RA) plays major roles in the development of the nervous system. mRNA-sequencing data of the RA-induced NTDs mouse model were obtained as reported previously14. We established the NTDs mouse model by intragastric administration of 28 mg/kg RA on pregnant mouse at E7.5, the control pregnant mice were treated with the equivalent volume of sesame oil. The brain vesicle of mouse embryos were collected at E8.5, E9.5 and E10.5 respectively, and sequenced by using Illumina HiSeqTM 2000 sequencing platform (SRP070626).

Human NTDs sample chip data

mRNA-sequencing data of human NTDs samples were obtained in our previous study. Three cases of spina bifida detected by ultrasonography at 17–24 gestational weeks and three control embryos from spontaneous abortions at the corresponding gestational weeks were collected. RNA was extracted from embryonic brain tissue, and transcriptome microarrays were performed using AffymetrixHG-U133A2.0 gene chip technology15.

Human GBM tumor database data

The GBM tumor transcriptome sequencing data in this study were obtained from the Gene Expression Profiling Interactive Analysis (GEPIA) database16,17. The GEPIA database enables a dynamic analysis of gene expression profile data. Additionally, it is an open public database developed for cancer and normal gene expression profile analysis. GEPIA include the RNA sequencing expression data of 9736 tumors and 8587 normal samples from the Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases. These expression data were re-calculated using the unified standard process, which can directly conduct a comprehensive expression analysis.

Animals

Institute of Cancer Research (ICR) mice (9–11 weeks, 19–25 g) were provided by the Animal Center of Shanxi Medical University and were housed in pathogen-free cages at the animal center with a 12 h light/dark cycle, maintaining the temperature and humidity within 20–26°℃ and 40–70%, respectively. NTDs mice were induced by intragastric administration of 28 mg/kg RA to pregnant mice at E7.5. Control mice were treated with the equivalent volume of sesame oil. The brain tissues of mouse embryos were collected at E9.5 and E10.5. Total RNA was extracted using TRIzol reagent (Ambion, USA), and cDNA was synthesized using the Revert Aid First Strand cDNA Synthesis Kit (Takara, Japan). Pregnant mice were killed by cervical dislocation, the experiments were conducted in accordance with the guidelines for animal care of the European Union (ARRIVE). In addition to this, all animal experiments are conducted in accordance with relevant laws and regulations on experimental animals and the rules and regulations related to experimental animals of Shanxi Medical University, and have been approved by the Ethics Committee of Shanxi Medical University (SYDL2021022).

GBM samples

All GBM samples were collected from the neurosurgery department of the Shanxi Cancer Hospital. Normal adjacent tissues were obtained from the extended excision ___domain of the GBM neurosurgical resection. All the human tissue samples and information used in this study were from the previous cases, which will not pose any risk to the subjects, we have received the exemption from informed consent. Ethical approval was obtained from the corresponding review board of the Shanxi Cancer Hospital (0002102).

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)

The Database for Annotation, Visualization, and Integrated Discovery (Version: v2023q1) was used for gene functional enrichment, including GO18 and KEGG (Version: Release 104.1) analysis19.

Protein–protein interaction (PPI) network analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database was used for PPI analysis20. Based on the mutual relationship results calculated from the STRING database, the Cytoscape3.42 software was used to construct the visual network topology, and the maximal clique centrality (MCC) algorithm of the cytoHubba plug-in was used to screen the hub genes.

Quantitative reverse transcription polymerase chain reaction (RT-qPCR)

The mRNA expression levels of candidate genes in the NTDs mouse and human GBM samples were analyzed using RT-qPCR. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), and cDNA was synthesized using a Revert Aid First Strand cDNA Synthesis Kit (Thermo, USA). qPCR was performed using Maxima SYBR Green/ROX qPCR Master Mix (Takara, Japan). The data were analyzed using the 2–∆∆Ct method. All experiments were repeated three times. The mRNA levels were normalized with GAPDH. The primer sequences are listed in Supplementary Table 1.

Western blot analysis

Proteins from the corresponding samples were extracted and subjected to SDS–polyacrylamide gel electrophoresis. The proteins were then transferred to polyvinylidene fluoride (PVDF) membranes and treated with the corresponding primary and secondary antibodies. Finally, the membrane was rinsed with TBST and treated with enhanced chemiluminescence solution. The primary antibodies for NTDs mouse samples were antibodies against stanniocalcin 1 (Stc1, 1:1000, Abclonal, China), dihydropyrimidinase like 4 (Dpysl4, 1:1000, Thermofisher, USA), DnaJ heat shock protein family (Hsp40) member C6 (Danjc6, 1:5000, Proteintech,USA), calpain-3 (Capn3, 1:1000,Proteintech,USA), paired like homeobox 2A (Phox2a, 1:1500, BIOSS, China), Polymerase iota (Poli, 1:1000, Proteintech,USA), Fibroblast growth factor 1 (Fgf1, 1:1000,Abclonal,China), and secondary antibodies were Horseradish Peroxidase (HRP)-conjugated AffiniPure Goat Anti-Rabbit IgG (1:4000, ZSGB‐BIO, China) and Goat Anti-Mouse IgG(1:4000, ZSGB‐BIO, China). The primary antibodies for human GBM samples were antibodies against Dnajc6 (1:1000, Abcam, UK), Poli (1:2000, Thermofisher, USA), Fgf1 (1:1000, Abcam, UK), Stc1 (1:1000, Abcam, UK), and secondary antibodies were HRP-conjugated goat anti-rabbit IgG (1:5000, Abcam, UK); HRP-conjugated rabbit anti-mouse IgG (1:5000, Abcam, UK). Images were obtained using a chemiluminescent blot detection system (ChemiDoc™ Imaging Systems, BIO-RAD, USA) and the results were further analysed using Image-J (Version 1.44p) software (open source from National Institutes of Health,USA)21. By using the Gels tool in ImageJ to measure the grayscale of the bands, peak maps can be obtained, and the Area value of each peak is the quantized value of each band. Normalizing the grayscale values of target gene/reference genes for each sample, performing t-test analysis on the treat and control group, and visualization was performed using Graphpad software.

Immunohistochemistry

Immunohistochemistry (IHC) was performed for the NTDs and control. We started the embedding machine to melt the paraffin used for immersion fixation of dehydrated tissues. After the paraffin was soaked and embedded, the slices were sectioned with a slicer (HistoCore BIOCUT, Leica, Germany). Subsequently, the sections were dewaxed, rehydrated, and placed in a pan with repair solution (1 mmol/L Tris–EDTA) and heated to 96 ℃. After cooling, the sections were treated with endogenous enzyme inactivation and blocked with 10% goat serum blocking solution. Primary antibodies were then added and incubated overnight at 4 ℃. After washing with PBS, the secondary antibodies were incubated for 1 h at room temperature. All sections were stained with 3, 3′-diaminobenzidine (DAB) solution and hematoxylin. The primary antibodies were antibodies against Stc1(1:100, Abclonal, China), Danjc6 (1:100, Proteintech,USA), Capn3 (1:100, Proteintech,USA), Phox2a (1:100, BIOSS, China), Poli (1:100, Proteintech,USA), Fgf1 (1:100,Abclonal,China), and secondary antibody was goat anti-mouse/rabbit IgG (ZSGB-BIO, China).After dehydration and gel sealing, the sections were examined and photographed under a microscope (ECLIPSE Ci, Nikon, Japan). After image acquisition, ImageJ analysis was performed. Images were analyzed semi-quantitatively based on AverageOptical Density (AOD) value. AOD=integrated optical density/area, and higher AOD indicate higher immunohistochemical staining intensity. AOD in control sample was set to 1 to measure the differential expression of target proteins. t-test analysis and visualization was performed using Graphpad software.

Cell culture and treatment

The mouse hippocampal cell line (HT-22) was purchased from BeiNaBio (No. BNCC358041), and the cells were cultured in DMEM high glucose medium (GIBCO, No. C11995500BT) containing 1% penicillin–streptomycin (Solarbio, No. P1400), and 10% foetal bovine serum (Cellmax, No. SA211.02) at 37 ℃ in a 5% CO2 incubator. Cells were passaged every two days. Human GBM cell lines U87 and SWO-68 were obtained from the American Type Culture Collection ATCC (Manassas, VA, USA), and were cultured with the DMEM containing 10% fetal bovine serum (FBS, Invitrogen; Thermo Fisher Scientific, Inc. Waltham, MA, USA) at 37 ℃ with 5% CO2. The cells were treated with 100 μmol/L retinoic acid (Sigma, Merck KGaA, USA), and 120 μmol/L pazopanib (TargetMol, USA) for 48 h.

CCK-8 analysis

HT-22 cells were treated with increasing doses of pazopanib for 24 h. At the end point of this assay, 10 μL of CCK-8 solution (meilunbio, No.MA0225) were added into the 100 μL medium and incubated for additional 2–3 h at 37 ℃. The absorbance of each well were detected at 450 nm using the microplate reader (Bio-Tek Instruments, Inc., Winooski, VT, USA). GBM cells were treated with increasing doses of pazopanib for 48 h. At the end point of this assay, 10 μL of CCK-8 solution (C0038, Beyotime, shanghai, China) was added into the medium and incubated for additional 3 h at 37 ℃. The absorbance of each well was detected at 450 nm using the microplate reader (Bio-Tek Instruments, Inc., Winooski, VT, USA). Three replicates were set at each point. Cell viability = (absorbance of experimental wells—absorbance of blank wells)/ (absorbance of control wells—absorbance of blank wells) × 100%. P values were determined by the Two-way ANOVA analysis.

Flow cytometry analysis

HT-22 cells were seeded into six-well dishes with 1.8 × 105 cells per well. Cells were treated with 100 μmol/L retinoic acid, 100 μmol/L retinoic acid and 80 μmol/L pazopanib, 100 μmol/L retinoic acid and 120 μmol/L pazopanib for 24 h, respectively. After 24 h of intervention, cells were digested and cells were stained with EDTA-free trypsin (meilunbio, No. PEL061). Then digested with trypsin and collected after 24 h. Cells were stained with Annexin V-FITC/PI assay kit (meilunbio, item no. MA0220-1) and cell apoptosis were detected by flow cytometry (BD FACSAria III). Based on the fluorescence intensity values automatically collected by the flow cytometry analysis software BD FACSDiva Software v9.4, the software automatically calculates the proportion of dead, normal, apoptotic, and early apoptotic cells and visualizes the results.

GBM cells were seeded into the 35 mm dishes and treated with 10 μmol/L of pazopanib for 24 h. Cells were trypsinized and harvested for staining with the Annexin V-FITC/PI Detection kit (Sigma, Merck KGaA, USA) according to the manufacturer’s instructions. Cells were resuspended after centrifugation at 2,000 rpm for 5 min and then incubated with propidium (PI) and FITC for 15 min at room temperature (RT), respectively .Cell apoptosis at different conditions were detected by the BD Accuri C6 Plus flow cytometry (BD Biosciences, San Jose, CA, USA). The FITC + /PI- and FITC + /PI + fractions were considered as early and late apoptosis, respectively. Data were analyzed using the FlowJo software (Version X.10.0.7-1). P values were determined by the student’s t test, and all experiments were repeated three times.

Cell invasion assay

Transwell assay was performed to evaluate the cell migration. Glioma cells were seeded in the upper transwell chamber (Corning, USA) that were pre-coated with matrigel and contained serum-free medium. RPMI-1640 with 10% FBS were added into the lower chamber.Then cells were treated with 10 μmol/L of pazopanib. The chambers were placed in the CO2 incubator for 48 h. Cells in the upper chamber were fixed with 4% paraformaldehyde for 15 min at room temperature, and then stained with 0.1% crystal violet (Beyotime, shanghai, China) for 10 min. The invaded cells were counted under an inverted microscope (Olympus, Tokyo, Japan). The number of invaded cells was determined by counting stained cells from multiple randomly selected microscopic visual fields.

Prediction process

The protein crystal structure was obtained from the Protein Data Bank (PDB)22, and the protein crystal structure of small molecules were obtained from the ZINC database23, which is one of the largest organic small-molecule compound databases at present. The protein molecules were pretreated using the Autodock software24. Hydrogenation, charge calculation, atomic type addition, charge detection, and routable chemical bond detection were performed on the small molecules. The accurate search docking method was used to dock the protein crystal structure with small molecules. OpenBabel (Version 3.1.1)25 was used for format conversion and PyMol (Version 2.3.2) for visualization26.

Surface plasmon resonance assay by biacore

Based on surface plasmon resonance (SPR) technology, the Biacore platform was used to verify the docking between proteins and small-molecule drugs in vitro. The binding affinity of FGF1 protein and small molecule pazopanib were performed by Biacore X100 SPR system (Biacore, Cytiva). HBS-EP (0.01 M HEPES with 3 mmol/L EDTA and 150 mol/L NaCl, 5% DMSO, 0.05% surfactant P20, pH 7.4) were used as running buffer. FGF1 was diluted with 10 mmol/L NaAc buffer (pH 4.5) and then immobilized on CM5 chip using an amine coupling kit (Biacore, Cytiva). The final FGF1 immobilization level were 1941.7 RU. Subsequently, pazopanib was injected as analyte at various concentrations (0.34, 1.03, 3.09, 9.27, 27.8, 83.3 and 250 μmol/L) at a flow rate of 30 μL/min with a contact time of 60 s and a dissociation time of 60 s. 30 s injection of glycine–HCl (pH 2.5, 10 mmol/L) were performed as regeneration. Data were analyzed by Biacore evaluation software (X100 version 2.0) using affnity model.

Statistical analysis

Statistical analysis were performed using the GraphPad software, version 8.0. Data are expressed as the mean ± standard deviation (SD). The t-test and chi-square test were used for comparisons between two groups, while for multiple comparison tests, analysis of variance (ANOVA) followed by Tukey’s post hoc test were performed. Statistical significance was set at P < 0.05.

Results

Obtaining genes related to NTDs and performing functional enrichment analysis

In order to obtain the genes related to NTDs, the data from RA-induced NTDs mouse and human NTDs samples were performed to select the candidate genes. The results showed that 162 genes in both RA-induced NTDs mouse and human NTDs samples were differentially expressed according to the criterion of | log2FC |> 1, P < 0.05 (Fig. 1A,B). Among the 162 genes, 78 genes were upregulated and three genes were downregulated in all four sets (human samples, E8.5, E9.5, E10.5 mouse samples). To better explore the role of NTDs related DEGs, the analysis of NTDs related DEGs for enrichment function analysis GO, KEGG analysis were performed.

Fig. 1
figure 1

DEGs from human NTDs sample and RA-induced NTDs mouse embryo. (A) Venn diagram; blue represents genes differentially expressed in human NTD samples (n = 3), pink represents genes differentially expressed in NTD mouse E10.5 embryonic brain tissues (n = 9), orange represents genes differentially expressed in NTD mouse E8.5 embryonic brain tissues (n = 9), and yellow represents genes differentially expressed in NTD mouse E9.5 embryonic brain tissues (n = 9). (B) The heatmap cycle is of differentially expressed genes. Different colors represent log2FC values. The circles from outside to inside are: gene name, log2FC value in NTDs mouse E10.5, E9.5, E8.5 and human NTD samples.

Full size image

GO analysis categogries included biological process, cellular component and molecular function. As shown in Fig. 2, at the biological process level, the NTDs related DEGs were enriched in kinds of development processes such as animal organ development, system development, multicellular organismal processes, tissue development, anatomical structure morphogenesis and cellular developmental process. At the cellular component level, NTDs related DEGs were enriched in the neurons, synapse, plasma membrane, post synapse, and cell periphery. At the molecular function level, NTDs related DEGs are mainly enriched in DNA-binding activator activity, RNA polymerase II-specific, molecular transducer activity, signaling receptor activity, and calcium ion binding. The results of KEGG enrichment analysis showed that NTDs related DEGs were mainly enriched in metabolism, environmental information processing, organismal systems, and human disease (Fig. 3A). At metabolism level, the NTDs related DEGs mainly enriched on tyrosine metabolism and other amino acid metabolism process, at the level of environmental information processing, the NTDs related DEGs mainly enriched on MAPK, Neuroactive ligand-receptor interaction, ErbB , TNF and other signaling pathways, at organismal systems level, the NTDs related DEGs mainly related to endocrine, nervous and immune systems, at human disease level, genes mainly related to pathways in cancer (Fig. 3A, Table 1). The 162 gens enrichment function analysis results remind us that there may be a potential relationship between NTDs and cancer that we haven’t yet discovered. What about the hub genes function and the relationship with cancer is unknown, we performed PPI and hub gene finding analysis next.

Fig. 2
figure 2

GO enrichment of 162 genes relating to NTDs. (A) GO functional classification. Red represents biological process, green represents cellular component, and blue represents molecular function. (B–D) Top 20 GO Enrichment terms in biological process, cellular component and molecular function respectively.

Full size image
Fig. 3
figure 3

KEGG pathway enrichment analysis of NTDs relating genes. (A) Enrichment cycle diagram of Grade I categories. (B) Statistical charts of Grade II categories of each pathway.

Full size image
Table 1 KEGG pathway enrichment of NTDs relating genes.
Full size table

A protein interaction network of 162 genes were constructed at STRING database and visualized by Cytoscape, and we found a very close and complex network interaction between were found among the NTDs related DEGs (Fig. 4A). The key modules in the complex network composed of these genes are shown in Fig. 4B. We obtained the top 10 hub genes, such as epidermal growth factor receptor (Egfr), SRC proto-oncogene (Src), erb-b2 receptor tyrosine kinase 2 (Erbb2), catenin beta 1 (Ctnnb1), kinase insert ___domain receptor (Kdr), erb-b2 receptor tyrosine kinase 3 (Erbb3), Jun proto-oncogene (Jun), erb-b2 receptor tyrosine kinase 4 (Erbb4), cadherin 1 (Cdh1), fms related receptor tyrosine kinase 1 (Flt1). Two of them are proto-oncogenes, Src and Jun. Four of them are the member of ERBB family of receptor tyrosine kinases, consisting of Egfr, Erbb2, Erbb3 and Erbb4, were first implicated in cancer in the beginning of the 1980s27. As vascular endothelial growth factor (VEGF) receptors, Flt1 and Kdr also play a crucial role in cancer progression individually28. These suggest that there may be some unrecognized relationship between neural tube closure disorders during embryonic development and tumorigenesis. Since neural tube closure primarily involves the cell movements of neuroectoderm, we chose primitive neuroectodermal tumors, GBM, for our subsequent study.

Fig. 4
figure 4

Protein–protein interaction (PPI) network of 162 NTDs relating genes and top module. (A) PPI network of 162 NTDs relating genes, colors other than green represent the top 10 hub genes. (B) Interaction status of the top 10 hub genes and the ranking of interaction strength.

Full size image

Screening for genes aberrantly expressed in both NTDs and GBM

Previous studies showed that tumor occurrence and embryonic development also have the same mechanism1,2,3,4,5,6,7,8,9. GBM, as primitive neuroectodermal tumors, develop from the ectoderm of early embyros. To explore whether have the same mechanism between NTDs and GBM, we determined that 162 NTDs-related DEGs were involved in GBM. The results showed that 68 genes showed differential expression in GBM tumors with |log2FC|> 1 and P < 0.05. Of these, 38 genes were up-regulated and 30 genes were down-regulated in tumor tissues (Fig. 5). Moreover. 23 genes in all five sets (NTDs mouse embryos at E8.5, E9.5, E10.5, human NTDs, and human GBM tumor samples) showed consistent differential expression profiles (Fig. 6). Finally, these 23 genes were used as candidate genes for verification in subsequent experiments.

Fig. 5
figure 5

Gene expression in GBM samples. (A) Total of 68 genes out of 162 genes showed differential expression in GBM samples in the GEPIA database. 38 up-regulated genes in GBM samples (n (tumour) = 163; n (normal) = 207). (B) 30 down-regulated genes in GBM samples (n (tumour) = 163; n (normal) = 207). *P < 0.01 vs normal.

Full size image
Fig. 6
figure 6

Heatmap of 23 candidate genes based on the log2FC value. Different colors indicate the log2FC value of this gene in the 5 transcriptome data sets (human NTD n = 3, human GBM n = 163, NTD mouse E10.5 n = 9, NTD mouse E9.5 n = 9, NTD mouse E8.5 n = 9). The figure shows 23 genes with the same expression trend, that is, the expression is up-regulated or down-regulated in all 5 transcriptome data sets.

Full size image

Experimental verification proves that Poli and Fgf1 are aberrantly expressed in both NTDs mouse and human GBM samples

Of the 23 candidate genes, seven were consistent with sequencing data both at E9.5 and E10.5, of which four were upregulated in NTDs samples (Fig. 7A,B), namely Stc1, Fgf1, DNA Poli, and Phox2a, and three were downregulated in NTDs samples (Figs. 2, 3), namely Dpysl4, DNAJ heat shock protein family (Hsp40) Danjc6, and Capn3. In addition to this, six genes showed a consistent trend with the sequencing results in the verification of the E10.5 sample (Fig. 7C). Overall, the validation consistency rate of sample E10.5 was higher than that of the E9.5 sample.

Fig. 7
figure 7

The differential expression of candidate genes in mouse embryo. (A) The expressions of up-regulated genes in E9.5 and E10.5 mouse embryos were determined by RT-qPCR, n = 10, *P < 0.05, **P < 0.01, ***P < 0.001 vs control. (B) The expressions of down-regulated genes in E9.5 and E10.5 mouse embryos were determined by RT-qPCR, n = 10, *P < 0.05, **P < 0.01, ***P < 0.001. (C) The expressions of candidate genes in E10.5 mouse embryos were determined by RT-qPCR, n = 10, *P < 0.05, **P < 0.01, ***P < 0.001.

Full size image

To better determine the candidate genes, we used the western blot to measure the expressions of the seven genes with high consistency in NTDs samples, and six genes at the protein level were consistent with that at the mRNA level (Fig. 8A,B). Meanwhile, the expression of the proteins encoded by Fgf1, Phox2a, Poli, Stc1, Capn3, and Dnajc6 in E9.5 and E10.5 mouse embryos were performed by immunohistochemistry. The results showed that the expression of the proteins FGF1, PHOX2A, POLI, STC1 in the NTDs mouse all showed up-regulated expression trends, and CAPN3 and DNAJC6 showed down-regulated expression trends in NTDs mouse compared to control (Fig. 8C). Among them, there was no significant difference in the up-regulation expression of POLI in the E9.5 NTDs mouse, and there were no significant difference in the down-regulation expression of DNAJC6 in the NTDs mouse (Fig. 8C).

Fig. 8
figure 8

Protein expressions were detected by western blot and immunohistochemistry in E9.5 and E10.5 mouse embryos. (A, B) The levels of proteins were detected by western blot in E9.5 and E10.5 mouse embryos, n = 10, *P < 0.05, **P < 0.01, ***P < 0.001 vs control. (C) The levels of proteins were detected by immunohistochemistry in E10.5 mouse embryos. *P < 0.05, **P < 0.01, ***P < 0.001, ns P > 0.05 vs control. All experiments were repeated three times. Indicates a significant difference compared to other groups in a two-way ANOVA followed by Tukey’s test (P < 0.05).

Full size image

Eleven of the 23 candidate genes showed consistent trends in GBM tumor tissues, namely, signal transducer and activator of transcription 5A (Stat5a), Dpysl4, periostin (Postn), homeobox D9 (Hoxd9), homeobox A7 (Hoxa7), Fgf1, Poli, mitogen-activated protein kinase 6 (Map3k6), potassium voltage-gated channel interacting protein 1(Kcnip1), Capn3, and inositol monophosphatase 2 (Impa2). The results showed that Stat5a, Postn, Hoxd9, Hoxa7, Poli and Fgf1 significantly highly expressed in GBM tumor tissues (Fig. 9A, p< 0.001).

Fig. 9
figure 9

The differential expression of candidate genes in GBM samples at the mRNA and protein levels. (A) The mRNA levels of all genes in GBM and normal tissues were measured using RT-qPCR (n (tumour) = 10; n (normal) = 10). The data are shown as the mean (n = 3). *P < 0.05. **P < 0.01, ***P < 0.001, ns P > 0.05. (B) The levels of DNAJC6, STC1, FGF1, POLI were detected by western blot in GBM tumor and normal tissues, n = 5,*P < 0.05, **P < 0.01, ns P > 0.05. All experiments were repeated three times. Indicates a significant difference compared to other groups in a two-way ANOVA followed by Tukey’s test (P < 0.05).

Full size image

In GBM samples, four genes were selected for protein level verification: Poli, Fgf1 (at the mRNA level, the expression was consistent with the sequencing results), Stc1, and Dnajc6 (at the mRNA level, the expression was opposite to the sequencing results). The verification results are shown in Fig. 9B. At the protein level, these four genes were consistent with the mRNA level. POLI and FGF1 were highly expressed in tumor tissues at both the mRNA and protein levels (Fig. 9A,B).

To sum up, among the 23 candidate genes, two were well confirmed in both NTDs and GBM tumor samples at both mRNA level and protein level, namely, Poli and Fgf1. Poli is a member of the DNA polymerase Y family, which can repair DNA damage29, but with unstable expression30. Fgf1 regulates the proliferation and differentiation of neuronal cells31. As the function of them, we guessed that Fgf1 might directly regulated cell activity, so we selected Fgf1 for further study via Molecule Docking and cytological experiment.

Pazopanib is a small molecule that binds to Fgf1 and has the effect of inhibiting cell proliferation and promoting cell apoptosis

As Fgf1 gene was well confirmed and regulates the proliferation and differentiation of neuronal cells. Molecular docking was performed for FGF1. In the PDB, the file of the FGF1 protein crystal structure (2hw9) were downloaded. In the DrugBank32 database, there are 12 small molecule drugs that can target and bind to FGF1. Four of these drugs are already approved for clinical treatment, with pazopanib being the antineoplastic agent used in the treatment of kinds of cancers33. Given the fact that pazopanib exhibits anti-tumor activity, we speculated that it may have effects on cell activities and has potential clinical application value. So in our study, pazopanib was selected for further study. The minimum binding energy of 2hw9 to pazopanib among all binding positions was-2.49, and one hydrogen bond was formed at the LYS101 position. The results of in vitro docking experiments showed that FGF1and pazopanib had a target-binding relationship, KD value = 9.962E–7 (Fig. 10).

Fig. 10
figure 10

Prediction and validation of Fgf1 and Pazopanib targeting binding relationship. (A) Prediction result of targeting binding relationship between FGF1 and Pazopanib. (B) Validation result of targeting binding relationship between FGF1 and Pazopanib.

Full size image

To investigate the effects of Pazopanib on cell proliferation, apoptosis, and migration. CCK-8, flow cytometry, and transwell experiments were conducted. Through the CCK-8 assay, we quantitatively analyzed the impact of Pazopanib on the proliferation of U87 and SWO68 cells. The results demonstrated that as the concentration of Pazopanib increased, the proliferative activity of these two glioma cell lines was significantly inhibited, with a clear dose-dependent effect (Fig. 11A). Furthermore, we assessed the apoptosis levels of U87 and SWO68 cells using flow cytometry. It was found that after treatment with Pazopanib, the proportion of cells in early apoptosis (Q2 quadrant) and late apoptosis (Q3 quadrant) significantly increased compared to the control group (DMSO treatment), indicating that Pazopanib can induce apoptosis in these cell lines (Fig. 11B). In terms of cell migration ability, we evaluated the impact of Pazopanib on the migration of U87 and SWO68 cells through experiments. The results showed that compared to the control group, the migration ability of cells treated with Pazopanib did not significantly differ, suggesting that Pazopanib has little effect on the migration ability of these two cell lines (Fig. 11C). Additionally, we conducted similar experiments on the HT-22 cell line. The results indicated that Pazopanib also has an inhibitory effect on the proliferation of HT-22 cells, and this inhibitory effect is enhanced with increasing drug concentration (Fig. 11D). Furthermore, we explored the impact of the combined intervention of Pazopanib and retinoic acid on the apoptosis levels of HT-22 cells, and we found that the combined use significantly increased the apoptosis levels of HT-22 cells (Fig. 11E).

Fig. 11
figure 11

Effect of Pazopanib on proliferation, apoptosis and migration in GBM tumor cell lines and HT-22 cells. (A) Cell proliferation were detected by CCK-8 in lines-U87 and SW68, ***P < 0.001. (B) Cell apoptosis were detected by flow cytometry in U87 and SW68 cells, **P < 0.01, ***P < 0.001. (C) Cell migration were detected in U87 and SW68 cells, ns P > 0.05. (D) Cell proliferation were detected by CCK-8 in HT-22 cells. **P < 0.01. (E) Cell apoptosis were detected by flow cytometry in HT-22 cells under the co-intervention of RA and PAZ, ***P < 0.001. All experiments were repeated three times. Indicates a significant difference compared to other groups in a two-way ANOVA followed by Tukey’s test (P < 0.05).

Full size image

Discussion

The NTDs susceptibility genes screened in animal models cannot fully explain the mechanism in humans. Current studies have shown that more than 200 genes are associated with NTDs in mouse models, but many of them have not been validated in human NTDs samples34. There are two important reasons for this. First, many related genes found in animal models are evaluated using the knockout method, and this type of experiment often leads to a syndrome-like phenotype rather than a single defect phenotype. However, human NTDs are mostly sporadic and non-syndromic with high genetic complexity. Secondly, the severity of malformations caused by gene knockouts in animal models is relatively high, and the proportions of different malformation categories are quite different from those in humans. In this study, we screened NTDs susceptibility gene based on both the NTDs mouse transcriptome sequencing data and the microarray expression profiling data of human NTDs samples. The limitation of our results is that the NTDs mouse transcriptome data used was only based on a mouse model induced by RA.

RA-induced NTDs mouse model is only one of many NTDs mouse models. While it cannot represent the NTDs mouse model, it is particularly irreplaceable. Since 1986, RA has been reported to be the most powerful teratogenic form of vitamin A and has been shown to induce various types of developmental defects in mice35. A study in 1988 has already found that RA caused some damage to embryos, including changes in the position of the spinal cord and the neural tube shape36. Furthermore, the impact on the whole CNS is far-reaching37. In early embryonic development, either excess or lack of RA can lead to neural tube closure disorders38. Such as, excess RA can elicit NTDs in rats and mice, ranging from caudal spina bifida aperta to encephalocele and a combination of these NTDs39,40. By contrast, lack of RA synthesis elicits craniorachischisis in mice41, emphasizing the importance of optimum levels of RA in convergence and extension of the neural plate during early neural tube closure.

Some classic studies have concluded that RA plays two major roles in the development of the nervous system. Firstly, RA is essential for patterning neural progenitor cells in the hindbrain and spinal cord. It contributes to both the anteroposterior and dorsoventral patterning of the neural plate and neural tube42. Normal neural tube patterning requires a balance of RA, Wnt and FGF signaling. FGF and RA pathways act antagonistically, and they have opposing effects on neuronal differentiation and ventral patterning onset in the newly formed neural tube43,44,45, while the Wnt signaling pathway can play a regulatory role in this process46. The second role RA has in the nervous system is neuronal differentiation. This has been studied extensively in in vitro models, such as embryonic carcinoma cells and neuroblastoma cells. RA induces the differentiation of various types of neurons and glia42 by activating the transcription of many genes47. In a study of glioma tumors, the expression of most RA signaling molecules increased with malignancy and was associated with augmented intratumoral retinoid levels in high-grade gliomas. The aberrant up-regulation of RA signaling molecules may be a central component in glioma pathogenesis48.

Based on the importance of RA in the neural tube and CNS development, we constructed an RA-induced NTDs mouse model to study the incorrect spatio-temporal expression of target genes, which will help in a deeper understanding of neural tube development. In previous studies, our research group has adopted high-throughput RNA sequencing technology to conduct mRNA-seq transcriptome profiles for the RA-induced NTDs mouse model. Since this study was based on the above data, the differentially expressed genes identified in this study could only be verified again in the RA-induced NTDs mouse model. Totally 162 intersecting genes were screened based on the NTDs mouse transcriptome sequencing data and the microarray expression profiling data of human NTDs samples, which could be considered as susceptibility genes closely related to NTDs.

Metabolic pathway enrichment of these NTDs related genes showed that more genes were enriched in cancer, MAPK signaling, and neuroactive ligand-receptor interaction pathways. These showed that some genes related to NTDs overlap with cancer-related genes, which highlighted the close relationship between early embryonic development and tumorigenesis.

After constructing the PPI network of NTDs related genes, 10 key genes were screened: Egfr, Erbb2, Erbb3, Erbb4, Src, Jun, Kdr, Flt1, Ctnnb1, Cdh1. Four of these are members of the Erbb family of receptor tyrosine kinases: Egfr, Erbb2, Erbb3, and Erbb4. All these genes are required for normal development and are involved in normal cell function49,50. EGFR is a growth factor receptor that induces cell differentiation and proliferation through ligand binding upon activation. EGFR is closely related to cancer occurrence51.

Two of them are proto-oncogenes, Src and Jun. The SRC family of protein kinases are critical for cell morphology, motility, proliferation, and survival52. SRC genes are expressed during embryonic development and tumorigenesis but are rarely expressed in normal cells53. The JUN also plays an important role in functions such as cell proliferation and survival54,55. Two are members of the VEGF receptor family: Kdr and Flt1. KDR is the main mediator of endothelial cell proliferation, survival, migration, tubule morphogenesis and budding induced by VEGF. The signal transduction and transport of this receptor are regulated by various factors, and the mutation of this gene is associated with infantil capillary hemangioma56. The protein encoded by the FLT1 gene binds to VEGFR-A, VEGFR-B and placental growth factor of the vascular endothelial growth factor receptor family and plays an important role in angiogenesis and angiogenesis57.

The gene Ctnnb1 encodes a protein that makes up the adherens junctions (AJs) protein complex. AJs are required to maintain the epithelial cell layer and regulate cell growth and intercellular adhesion. The protein it encodes also anchors the actin cytoskeleton and may be responsible for transmitting contact inhibition signals that cause cells to stop dividing once the epithelial layer is complete. AJs also are closely associated with cancer58,59.

Another gene, CDH1, encodes canonical cadherin of the cadherin superfamily. It is required for physiological signaling pathways such as cell proliferation, maintenance of cell adhesion, cell polarity, and epithelial-mesenchymal transition. Abnormal expression leads to tumor proliferation, invasion, migration, and metastasis60. All the above findings indirectly provide a basis for shared regulatory genes between embryonic development and tumorigenesis.

Both KEGG metabolic pathway analysis and top10 NTDs susceptibility gene analysis showed that a potential relationship between NTDs and tumorigenesis. It was suggested that there might be common genes involved in pre-embryonic development, middle and late embryonic development, as well as in individual maturation. Neural tube and neuroectodermal tumors have the same tissue origin, and it is reasonable to speculate that there might be common genes involved in both neural tube development and neuroectodermal tumor formation.

In the present study, susceptibility genes closely related to NTDs also showed differential expression in GBM tumors, and there were 23 genes with the same differential expression trend as NTDs samples, indicating that these genes may have potential effects on the development and differentiation of neuroectoderm. Additionally, dysregulated expression leads to abnormal phenotypes, which can provide new target information for the treatment of these diseases.

In this study, there were two genes whose mRNA and protein levels had been consistently verified in NTDs samples and GBM tumor samples, namely, Poli and Fgf1. Poli is a member of the DNA polymerase Y family and has translesion synthesis activity, which can repair various types of DNA damage61. DNA repair genes are essential for maintaining DNA integrity and normal function during early fetal development62. In contrast, tumor cells may escape cell death by enhancing their capacity to repair DNA damage. Fgf1 regulates the proliferation and differentiation of neuronal precursor cells within the neural tube63. The SV40 T antigen driven by the 1B promoter (−540 to + 31) of the human Fgf1 gene causes tumorigenesis in the brains of transgenic mice. Fgf11/FGFR signaling is associated with the expression of Aurora A and activation of its kinase ___domain (Thr288 phosphorylation), thereby maintaining the stem cell properties of GBM cells64.

Starting from different purposes and using different analysis methods, this study also found that the high expression of Fgf1 in abnormal mouse embryonic brain tissue and adult brain GBM tumor tissue coincides with previous research findings. The current understanding of the role of FGF signaling in a variety of biological and developmental processes has improved considerably over the past few decades. The FGF/FGFR system affects the pathophysiology of many human diseases, including hereditary diseases, metabolic diseases, and cancer65. In this study, it was found that Fgf1 are aberrantly expressed in both NTDs mouse and human GBM. To further investigate the value of Fgf1 in clinical applications, we found small molecule drugs that can bind to Fgf1 and also demonstrated the impact of this small molecule drug on cell proliferation, migration and apoptosis at the cellular level.

The database showed that Fgf1 may target and bind to pazopanib, and its docking site was predicted by molecular docking in this study. Studies have confirmed that pazopanib can target the vascular endothelial growth factor receptor (VEGFR) and inhibit the formation of new blood vessels. FLT1, one of the 10 NTD-related key genes identified in the this study, encodes a protein that can bind to VEGFR-A and VEGFR-B of the vascular endothelial growth factor receptor family. Therefore, pazopanib may plays a complex role in the regulation of key NTDs genes. Considering this, interaction networks may play an important regulatory role. Pazopanib has been approved for the treatment of advanced renal cell carcinoma, soft tissue sarcoma, epithelial ovarian cancer, and non-small cell lung cancer33. In a study on GBM tumors, pazopanib was therapeutically effective against a class of GBM tumors (classification based on molecular mutational heterogeneity) by targeting CSF1R55.

In summary, Poli and Fgf1 aberrantly expressed in both NTDs mouse embryo and GBM human samples, they may influence cell activities both in embryonic development and adult development. They likely influence the development of the neuroectoderm. And might be involved in the occurrence of NTDs and GBM tumor cells. Pazopanib may affect the network interaction of NTDs-related genes and may be a new drug for the treatment of GBM tumors.

However, in our study NTDs samples and GBM tumor samples were collected separately, and it is unclear whether the NTDs mouse model in this study will develop GBM tumor. There are some reasons for collecting NTDs and GBM samples separately, firstly because mouse embryos develop in utero, in our laboratory setup, it is difficult to judge whether the NTDs model is successfully constructed in vitro and therefore cannot verify whether NTDs mouse that have been killed will suffer from GBM in the subsequent growth. Secondly, patients with both diseases are rare in clinic. This is because many children with NTDs have a short life span. One out of ten babies with NTDs will die before their first year66. For most of the fetuses diagnosed with NTDs, the parents choose to terminate the pregnancy before 28 weeks of pregnancy67. GBM is also a disease with a high mortality rate, and the 5-years survival rate is less than 10%68. The life cycle of patients with both diseases is very short, so it is difficult to obtain samples with both diseases in the same patient. Although we did not obtain such samples, it does not prevent us from studying the shared genes and potential mechanisms of the two diseases.

In addition, our team is currently collaborating with the Capital Institute of Pediatrics and Shanxi Cancer Hospital to closely track the health of children with NTDs, and conduct regular follow-up visits. At the same time, epidemiological investigations are also carried out for patients with GBM to obtain important information, such as their medical history. If samples with both diseases are obtained, we will start the study immediately.

Conclusion

In this study, there were two genes whose mRNA and protein levels have been well verified in both NTDs mouse samples and GBM human samples, namely, Poli and Fgf1. Fgf1 may target and bind pazopanib. In the present study, the docking site was predicted and verified in vitro. Cytological experiments also showed that pazopanib affects the proliferation and apoptosis of nerve and tumor cells. Overall, Poli and Fgf1 may play an undiscovered role in neuroectodermal development. Therefore, deregulation of these genes might be associated with the occurrence of NTDs during embryonic development and GBM tumors during adult development. Fgf1 may be a potential biomarker for NTDs. Pazopanib may be a new drug for the treatment of GBM tumors.