Zijian Li
Juan Feng
Zhengwei Yuan2*†- 1Department of Neurology, Shengjing Hospital of China Medical University, Shenyang, China
- 2Key Laboratory of Health Ministry for Congenital Malformation, Shengjing Hospital of China Medical University, Shenyang, China
Spina bifida is a common neural tube defect (NTD) accounting for 5–10% of perinatal mortalities. As a polygenic disease, spina bifida is caused by a combination of genetic and environmental factors, for which the precise molecular pathogenesis is still not systemically understood. In the present study, we aimed to identify the related gene module that might play a vital role in the occurrence and development of spina bifida by using weighted gene co-expression network analysis (WGCNA). Transcription profiling according to an array of human amniocytes from patients with spina bifida and healthy controls was downloaded from the Gene Expression Omnibus database. First, outliers were identified and removed by principal component analysis (PCA) and sample clustering. Then, genes in the top 25% of variance in the GSE4182 dataset were then determined in order to explore candidate genes in potential hub modules using WGCNA. After data preprocessing, 5407 genes were obtained for further WGCNA. Highly correlated genes were divided into nineteen modules. Combined with a co-expression network and significant differentially expressed genes, 967 candidate genes were identified that may be involved in the pathological processes of spina bifida. Combined with our previous microRNA (miRNA) microarray results, we constructed an miRNA–mRNA network including four miRNAs and 39 mRNA among which three key genes were, respectively, linked to two miRNA-associated gene networks. Following the verification of qRT-PCR and KCND3 was upregulated in the spina bifida. KCND3 and its related miR-765 and miR-142-3p are worthy of further study. These findings may be conducive for early detection and intervention in spina bifida, as well as be of great significance to pregnant women and clinical staff.
Introduction
Spina bifida is one of the most common central nervous system (CNS) malformations found in fetuses: the neural tube shows incomplete closure, and is usually complicated by other neurological abnormalities and comorbidities such as hydrocephalus. Spina bifida can be broadly classified into two groups: spina bifida aperta and spina bifida occulta. The former refers to meninges or nerve tissue bulging out of the spinal canal through spina bifida and forms a cystic mass, so it also called spina bifida cystica, including myelomeningocele, meningocele, myelocele, etc. In general, there is no need to treat spina bifida occulta. Spina bifida mentioned in this study refers to spina bifida aperta. Spina bifida is often either treated after birth or an intrauterine defect repair is undertaken (Alabi et al., 2018). However, complications such as hydrocephalus, fecal dysfunction, voiding dysfunction, and limb mobility disorders cannot be completely avoided. To date, the etiology of spina bifida is not fully understood and may be related to genetic and environmental factors. Elucidation of the pathogenesis of spina bifida is urgently required in order to reduce harm to the fetus, as well as that to society and the family.
Low folic acid levels in pregnant women, or antiepileptic drugs such as valproic acid taken during pregnancy, and a history of obesity and diabetes in pregnant women are environmental risk factors for the disease. Folic acid supplementation in pregnant women, a prenatal diagnosis, and fetal in utero therapy can effectively reduce the number of children born with birth defects (Molloy et al., 2017). Evidence has shown that in spite of sufficient folic acid intake by pregnant women, impaired uptake and utilization of folic acid can also contribute to neural tube malformation (Gonseth et al., 2015). Gene mutations in methylene–tetrahydrofolate reductase and methionine synthase reductase lead to a folate metabolic disturbance. Reduced folic acid uptake causes the down-regulated expression of transcription factor AP-2 alpha, which may increase the level of DNA methylation in the fetus. All of these are likely genetic risk factors that affect the incidence of neural tube malformation.
The establishment and improvement of disease animal models are important for revealing the pathogenesis of human spina bifida, and play key roles in testing novel interventions (Smith et al., 1978; Duru et al., 2001; Danzer et al., 2005). Recent years, researchers have explored the effects of some risk factors by using these animal models (Zhao et al., 2008; Qin et al., 2016), and take further measures such as gene and cell therapy strategies to intervene to reduce the occurrence of perinatal birth defects (Melo-Filho et al., 2009; Turner et al., 2013; Wei et al., 2020a). It is indicated that genetic therapy may become a novel strategy for treating birth defects. Therefore, it is necessary to analyze the gene expression profiles of fetuses with neural tube defects (NTDs), which allows us to further understand the genetic determinants and epigenetic factors involved in this pathological mechanism. Nagy and colleagues (Nagy et al., 2006) contributed a dataset of whole genome mRNA expression profiles from amniocytes in amniotic fluid samples of nine people (GSE4182: five spina bifida and four healthy control samples). They revealed three novel candidate genes: Src like adaptor (SLAP), leukocyte specific transcript 1 (LST1), and mal, T cell differentiation protein like (BENE), which can play an important role in the pathogenesis of NTDs.
However, more than a decade has passed since their results were reported, and many new and effective bioinformatics methods have appeared. Many questions are worthy of further exploration such as a study on gene co-expression modules in human amniocytes from patients with spina bifida. In this study, we used genes in a GSE4182 dataset that were in the top 25% of variance from spina bifida and healthy control amniocytes to construct a co-expression network by WGCNA. We attempted to identify promising candidate biomarkers or potential therapeutic targets of spina bifida from modules in which highly correlated genes clustered. Furthermore, a study of specific molecular and biological functions of these hub genes may be better for understanding the underlying mechanisms of the disease. Based on our previous findings of diagnostic miRNAs in the serum of pregnant women with fetuses that have NTDs, we reveal the candidate genes in a hub module and a miRNA–mRNA interaction network for spina bifida.
Materials and Methods
Microarray Dataset Collection
The Gene Expression Omnibus (GEO) database1 was used to obtain gene expression profiles for spina bifida. Exclusion criteria: (1) Profiles based on cell lines and animal models were excluded. Inclusion criteria: (1) Only homo sapiens species samples with spina bifida and healthy controls were included in this study. The GSE4182 dataset was detected on a platform of Affymetrix GeneChip Human Genome U133 Plus 2.0 [HG-U133_Plus_2]. GSE4182 included nine amniotic fluid samples from pregnant women with spina bifida (N = 4) or healthy (N = 5) fetuses. Amniocytes in the amniotic fluid were collected by amniocentesis, from which fetal mRNA was isolated and analyzed. In this study, the expression profile was acquired directly from a public database.
Data Preprocessing
The RAW data of the expression dataset, GSE4182, in .CEL format were obtained from a GEO database. The “affy” package in R was used to conduct the normalization and background correction of data. Probe level data were then converted into gene expression values. For multiple probes corresponding to a gene, the average expression value was taken as the gene expression value in this study. The distribution patterns of disease and control samples (before and after clustering analysis and outliers removement) were observed by principal component analysis (PCA).
Identification of Differentially Expressed Genes in Spina Bifida
The “limma” package in R was used to obtained Differentially Expressed Genes (DEGs) between spina bifida and healthy control samples in the expression data. We then carried out a significance analysis of microarrays and set the selection criteria as a false discovery rate (FDR) value <0.05 and log2| fold change| >1 (fold change >2) for further network construction.
Construction of Co-expression Network
The “WGCNA” package in R was used to construct the co-expression network based on the expression data profile of these top 25% variant genes. The microarray quality was checked by the “impute” package in R, which could detect whether the genes had missing values and ensure they were good samples. We performed sample clustering to plot the sample tree, and to detect and delete outliers. We then performed Pearson’s correlation matrices for pair-wise genes and found a soft thresholding power β value by using the pickSoftThreshold function of WGCNA.
Identification of Hub Genes
First, we defined the module with the largest absolute value of Pearson’s correlation of module membership (MM) and a p-value < 0.05 as the hub module. Furthermore, we defined hub genes in co-expression networks as genes that satisfy two conditions: the absolute value of the Pearson’s correlation of MM >0.8 and the absolute value of the Pearson’s correlation of gene trait (GS) relationship >0.2, which represented high module connectivity and high clinical significance, respectively. We subsequently obtained real hub genes by taking the intersection of hub genes in a co-expression network and significantly DEGs. On the basis of this, such real hub genes were used to construct an miRNA–mRNA regulatory network.
Functional Enrichment Annotation
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of real hub genes were performed by using the online program: Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Dennis et al., 2003). The GO terms were divided into biological process (BP), cellular component (CC), and molecular function (MF). Statistical significance was considered as a p-value < 0.05.
Prediction of miRNA Target Genes and Construction of miRNA–mRNA Regulatory Network
TargetScan2 is a miRNA target prediction database that is used for mammals by searching for the presence of conserved 8, 7, and 6 mer sites that match the seed region of each miRNA (Agarwal et al., 2015). As a database, miRTarBase3 contains more than 360,000 miRNA–target interactions that have been validated experimentally by reporter assay, western blot, microarray, and high-throughput sequencing (Chou et al., 2018). These two online tools were used to predict target mRNAs in our previous study of diagnostic miRNAs (miR-142-3p, miR-144, miR-720, miR-575, and miR-765) in the serum of pregnant women with fetuses that had spina bifida (Gu et al., 2012). We took the intersection of predicted target mRNAs and real hub genes, and used the results to construct an miRNA–mRNA regulatory network. Furthermore, the network was visualized by Cytoscape 3.7.1 (Shannon et al., 2003).
Samples Collection
The samples used in this study are from the sample bank of our department, which were collected as previously described (Zhang H. et al., 2019). The study was approved by the Ethics Committee of Shengjing Hospital, China Medical University. The ethics number is 2015PS264K. All consent was obtained to use the samples for testing. Six cases of spinal cord tissue were obtained from fetuses induced labor by spina bifida, and eight cases of spinal cord tissue were obtained from fetuses induced labor by non CNS congenital malformations which were approximately similar age used as normal controls, as shown in the Supplementary Table 1. Tissue samples from spina bifida and healthy control fetuses were collected at pregnancy termination by experienced pathologists and stored at −80°C until analysis.
Quantitative Real Time Polymerase Chain Reaction
Subsequently, quantitative real time polymerase chain reaction (qRT-PCR) was used to verify the expression of the key genes in spinal cord tissue of clinical samples. Total RNA was isolated from each sample using Invitrogen (15596026) Trizol reagent by our research team member according to the manufacturer’s instructions. Reverse transcription from total RNA to cDNA and qRT-PCR were performed using the Takara PrimeScript RT Master Mix (RR036A) and SYBR Green Premix (RR420A), respectively. The results were analyzed using the 2-ΔΔCt method and represented as fold changes, normalized to GAPDH. The PCR primers used in this study were shown in Supplementary Table 2. Statistically significant was considered as the p-value < 0.05.
Results
Data Preprocessing and DEG Identification
The raw data of a GSE4182 dataset that included four spina bifida and five healthy control samples was downloaded. Such data were then preprocessed by conducting format transformation, the filling in of missing data, background correction, and data standardization. The expression matrixes of a total 21,626 genes of these nine samples were obtained. After performing sample clustering to plot the sample tree, it was found that disease sample GSM94601 was clustered with other healthy samples. Then sample GSM94601 was removed as an outlier from a subsequent analysis (Figure 1A). We then matched the disease state of samples with their expression matrixes. The remaining eight samples were re-clustered and a sample dendrogram and trait heatmap were plotted (Figure 1B). At the same time, we performed PCA which showed that samples of spina bifida (except sample GSM94601) were distributed on the left side, while healthy control samples were distributed on the right side. The PCA results were identical with those by clustering analysis. The PCA results before and after outlier removement are shown in Figures 1C,D. The hierarchical cluster analysis heatmap showed significantly different distributions of gene expression patterns between spina bifida fetuses and healthy control samples (Figure 2A). Under the threshold of FDR < 0.05 and fold change >2, a total of 2,634 DEGs (1,277 up-regulated and 1,357 down-regulated) were chosen for further analysis, as shown in Figure 2B.
Figure 1. Data Preprocessing. (A) Samples clustering of total nine samples in the GSE4182 to detect outliers. (B) Re-clustering of the remaining eight samples: sample dendrogram and trait heatmap. The clustering was based on the expression data of DEGs between healthy controls and spina bifida samples. In the disease state, white color means healthy controls and red color means spina bifida. (C) PCA for spina bifida and healthy control samples before outlier identification and removal. (D) PCA for spina bifida and healthy control samples after outlier identification and removal.
Figure 2. Heatmap and volcano plot of DEGs. (A) The heatmap of DEGs. (B) The volcano plot of DEGs. Red represents upregulated genes and green indicates down-regulated genes. Black means nondifferentially expressed genes.
Co-expression Network Construction
Based on variance analysis, the top 25% of genes (5,407 genes) was obtained from GSE4182 with eight samples. These 5,407 genes were further analyzed and screened using WGCNA. We explored the value of the weight parameter β of an adjacency matrix by setting up a set of soft-thresholding powers (from 0 to 20). In this study, the soft-thresholding power was chosen as β = 8 where the curve first reached R^2 = 0.88, to construct a weighted network based on a scale-free network distribution (Supplementary Figures 1A–D).
The minimum number of genes was set as 30 for each module, and a clustering of module eigengenes obtained. We then calculated the dissimilarity coefficients between the genes and plotted the system cluster tree. Modules with a similarity above 70% were merged (Supplementary Figure 2A), and a dynamic tree dendrogram was redrawn (Figure 3A). Finally, a total of 19 modules was obtained using a dynamic tree-cutting method. The number of genes in each module is listed in Table 1.
Figure 3. (A) Dendrogram of the top 25% of genes clustered based on a dissimilarity measure (1-TOM). (B) Heatmap of the correlation between module eigengenes and spina bifida. (C) Genetic network heatmap based on topological overlap. (D) Scatter diagrams for module membership vs. gene significance of disease state in turquoise modules. (E) Heatmap of gene expression in the module and feature vector histogram.
Table 1. The number of genes in the 19 modules.
Identification of Clinically Significant Module
After relating modules to traits, high correlations were observed in traits of the disease state (healthy control or spina bifida) in the turquoise module (Supplementary Figures 2B,C). A clinically significant module was defined as a module with a maximum correlation coefficient (meanwhile p-value < 0.05). In this study, a turquoise module (correlation coefficient = 0.99, p < 1–200) was considered a clinically significant module and eligible for further analysis, as shown in Figure 3B. A genetic network heat map was obtained by calculating topological overlap between the top 25% of genes (Figure 3C). Each row and column of the heat map corresponds to a gene, with the darker the color, the higher the topological overlap, and the higher the density of genes. In Figure 3D, scatter diagrams for MM versus gene significance of disease state are shown in the turquoise modules. We then drew a heat map of gene expression in the turquoise module and its feature vector histogram (Figure 3E).
Identification of Hub Genes
First, the 1000 real hub genes were obtained by screening genes which | GS| >0.2 and | MM| >0.8 in the turquoise module (with total 1996 genes). Then, we overlapped these real hub genes with DEGs with fold change >2 (with total 2634 genes). Finally, the 967 candidate genes (as shown in Figure 4) were obtained for subsequent functional enrichment analysis and miRNA-mRNA network construction.
Figure 4. Venn plot of real hub genes in turquoise module and DEGs. The intersection in deep red represents the 967 genes that are common between the real hub genes in the turquoise module and DEGs.
Functional Enrichment Annotation
The results of GO and KEGG analysis of these candidate genes were shown in the Figure 5. It is worth noting that biological processes related to central nervous system development (GO:0007417), cell migration (GO:0016477), keratinization (GO:0031424), innate immune response (GO0045087), inflammatory response (GO:0006954), establishment of skin barrier (GO:0061436), keratinocyte differentiation (GO:0030216), and epidermis development (GO:0008544) were enriched in the top 20 of the results. It may be associated with the skin defect, spinal dysraphism (with or without myelomeningocele) of the spina bifida fetus. The chip samples were collected from amniotic fluid cells, and amniotic fluid cells are derived from fetal exfoliated skin. In a fetus with spina bifida, the edge of the skin of spina bifida may be exfoliated into amniotic fluid and detected. Therefore, the expression of genes in amniotic fluid cells is thought to reflect fetal conditions. The KEGG pathways analysis of these candidate genes indicated that pathways for neuroactive ligand-receptor interaction (hsa04080), NF-kappa B signaling pathway (has04064), osteoclast differentiation (hsa04380), cytokine-cytokine receptor interaction (has04060), and cell adhesion molecules (CAMs) (hsa04514) may play vital roles in the pathological processes of spina bifida.
Figure 5. Gene Ontology annotation and pathway enrichment analysis of significantly different genes in turquoise module. (A) Top 20 biological process (BP) terms of GO analysis. (B) Top 20 cellular component (CC) terms of GO analysis. (C) Top 20 molecular function (MF) terms of GO analysis. (D) KEGG pathway enrichment analysis. Red color means –log10 (p-value), blue color means the gene count.
Construction of miRNA–mRNA Regulatory Network
Our previous study (Gu et al., 2012) explored circulating miRNAs in pregnant women’s sera as potential biomarkers for spina bifida fetuses. We revealed that the expression of five types of miRNAs (miR-142-3p, miR-144, miR-720, miR-575, and miR-765) in the sera of spina bifida fetuses and pregnant women was up-regulated more than two times. Such miRNAs decreased in maternal serum 24 h after delivery, suggesting that maternal serum levels of these five types of dysregulated miRNAs are associated with pregnancy, and which could also reflect fetal conditions to some extent. MicroRNA target genes were predicted by TargetScan and miRTarBase, of which two or more miRNAs might target the same mRNA. The mRNAs appeared simultaneously in the 967 candidate genes and simultaneously predicted ranges were considered as credible mRNAs and shown in the miRNA–mRNA regulatory network (Figure 6). According to the above principles, miRNA-720 did not have credible predicted target genes and was absent in this network. The details of DEGs in credible miRNA–mRNA regulatory networks in spina bifida are shown in Table 2. Three key genes, YOD1, TSPAN6, and KCND3 were, respectively, linked to two miRNA-associated gene networks, of which their potential biological function and underlying mechanism are worthy of further discussion and research.
Figure 6. The miRNA-mRNA regulatory network in the spina bifida. The network consisted of four miRNAs (yellow) and 39 mRNAs (light pink). YOD1, TSPAN6, and KCND3 were considered as the three key genes connecting two miRNA gene networks (turquoise).
Table 2. DEGs in the confirmed miRNA-mRNA regulatory networks.
Validation in the Clinical Samples
To further validate the results of the microarray analysis, we examined the expression of the three key genes which have dysregulated expression in spina bifida by qRT-PCR using spinal cord samples from spina bifida and control fetus with induced labor with non CNS congenital malformations. Subsequently, we found that KCND3 and YOD1 were upregulated and TSPAN6 was downregulated in spina bifida samples. The unpaired t-test was used for statistical analysis of the data shown in Figure 7. However, only the difference of KCND3 between the two groups had statistical significance (p < 0.05), which was consistent with our previous analysis. TSPAN6 had the same decreased expression trend with our bioinformatics analysis results, but did not show a statistically significantly difference between the two groups. However, the changing trend of YOD1 was in the opposite direction.
Figure 7. Validation in the clinical samples by qRT-PCR. KCND3 (A) upregulated in spina bifida samples. TSPAN6 (B) and YOD1 (C) did not show significantly statistical difference between two groups. *p < 0.05 compared to the control group.
Discussion
Spina bifida is one of the most common serious birth defects found worldwide for which prenatal treatment options remain limited. Researchers have thus recognized the urgency for improvements in early preclinical diagnosis and treatment levels. In recent years, our research group has been devoted to the study of the genetic and epigenetic etiology and pathogenesis of spina bifida (Fan et al., 2011; Shan et al., 2012; Wei et al., 2013; Zhang H. et al., 2019; Zhang H.N. et al., 2019; An et al., 2020; Liu et al., 2020), as well as the diagnosis and prenatal diagnosis of diseases (Gu et al., 2012; An et al., 2015), gene therapy (Ma et al., 2020) and stem cell therapy (Li et al., 2012; Wei et al., 2020a,b) and other etiological treatments, which are conducive to improving the level of disease prevention and treatment. However, we still believe more needs to be done. In this study, WGCNA was used to explore pathological processes and marker genes in amniotic fluid cells of the spina bifida fetus. In this study, we identified a miRNA–mRNA regulatory network that consisted of four miRNAs and 39 mRNAs. TSPAN6, YOD1, and KCND3 are considered key genes in this network.
As a member of the tetraspanin family, tetraspanin-6 (TSPAN6) was predicted to be co-regulated by miR-142-3p and miR-144 in spina bifida. Guix and colleagues (Guix et al., 2017) reported that TSPAN6 was increased in Alzheimer’s disease brains and promoted Aβ-peptide accumulation by affecting the autophagosome–lysosomal pathway and slowing down the degradation of Amyloid precursor protein (APP)–C-terminal fragments. At the same time, TSPAN6 increased exosome-mediated secretion of APP–C-terminal fragments. TSPAN6 was found to be important for cognition and to affect properties of the postsynaptic terminal (Salas et al., 2017). Deletions at TSPAN6 cause epilepsy female-restricted with intellectual disability (Vincent et al., 2012). In addition, TSPAN6 affected mitochondrial antiviral signaling (MAVS) formation in a ubiquitination-dependent manner, which further negatively regulated the retinoic acid–inducible gene I-like receptor (RLR) pathway and host antiviral immune response (Wang et al., 2012). In general, all-trans retinoic acid was employed as a spina bifida aperta–inducing agent by our research group (Xue et al., 2018; Zhang H.N. et al., 2019) and other international peers (Oria et al., 2018). TSPAN6 may be tractable as a therapeutic target for spina bifida based on the results obtained from a literature search and data analysis.
As a gene encoding deubiquitinating enzyme, YOD1 was predicted to be co-regulated by miR-142-3p and miR-144 in spina bifida. Previous studies suggested that YOD1 correlated with another congenital malformed disease: nonsyndromic cleft lip, with or without cleft palate (NSCL/P). RNA interference and overexpression experiments indicated that YOD1 could enhance cell migration during lip and palate formation through the transforming growth factor (TGF)-β3 signaling pathway (Zhou et al., 2018). A mutation of YOD1 may lead to impeded cell migration resulting in NSCL/P (Ju et al., 2018). We therefore have reason to speculate that YOD1 plays a similar role in spina bifida. A few studies have highlighted the significant role of YOD1 in neurodegenerative disease. Tanji et al. demonstrated that the neurogenic proteins that cause Huntington and Parkinson’s diseases induced upregulation of the YOD1 level (Tanji et al., 2018). Ubiquitin-directed AAA-ATPase p97 cooperating with YOD1 played an essential role in the clearance of ruptured lysosomes by autophagy, in which a mutation caused inclusion body myopathy and neurodegeneration (Papadopoulos et al., 2017). Overall, such data suggest that the deubiquitinase YOD1 contributes to the pathogenesis of neurodegenerative disease by affecting the ubiquitin–proteasome system and autophagy–lysosome pathway. Data has also suggested that YOD1 plays an important role in the development of an antiviral immune response. Liu and colleagues reported YOD1-mediated K63-linked deubiquitination could activate an innate antiviral immune response against viral infection, and the aggregation of MAVS (Liu et al., 2019). In addition, a catalytically inactive mutant of the deubiquitinase YOD1 enhances antigen cross-presentation (Sehrawat et al., 2013). YOD1 antagonizes ubiquitin ligase TRAF6/p62-dependent interleukin-1 signaling to NF-κB (Schimmack et al., 2017). Past studies have suggested that YOD1 and TSPAN6 both affect MAVS formation and an active antiviral immune response by regulating the ubiquitination pathway. More speculatively, spina bifida may be related to a type of viral infection.
As a gene encoding potassium voltage-gated channel subfamily D member 3, KCND3 was predicted to be co-regulated by miR-142-3p and miR-765 in spina bifida. It was shown that a mutation in KCND3 was associated with the autosomal dominant inherited neurodegenerative disorder, spinocerebellar ataxia types 19 and 22 (SCA19/22) (Duarri et al., 2012; Lee et al., 2012; Paucar et al., 2018; Hsiao et al., 2019). Moreover, KCND3 pathogenic variants may be responsible for a wider phenotypic spectrum than previously thought, including autosomal recessive early-onset sporadic cerebellar ataxia (Kurihara et al., 2018), childhood epileptic encephalopathy (Wang et al., 2019), Parkinsonism, and cognitive impairment (Huin et al., 2017). However, further experimental studies are needed to validate a dysfunction of the potassium voltage-gated channel exists in spina bifida.
In this study, novel biomarkers of spina bifida have been provisionally identified from the viewpoint of coexpression network analysis. Our findings also provide novel insights into an miRNA–mRNA regulatory network in spina bifida. The genes we identified as key genes were derived from the hub module (with highest correlation with disease, shown in Figure 3B) by WGCNA. Furthermore, the miRNAs we used for constructing the miRNA-mRNA network were obtained from our previous verified dysregulated microarray results in spina bifida. Based on our preliminary experimental validation results, KCND3 was upregulated in the spina bifida group, which supports that KCND3 may play a role in the development of spina bifida. This result corresponded to that of our bioinformatics analysis. According to Figure 6, miR-765 and miR-142-3p may play a role by regulating the expression of KCND3 involved in the development of spina bifida. But how they participate in the disease onset remains unclear. These details need to be confirmed by further in depth experimental studies, and this will be an important research direction to uncover the molecular mechanisms of abnormal expression of these gene and miRNAs in spina bifida. Furthermore, future research and validation studies need to be based on larger sample sizes.
Conclusions and Future Prospects
The incidence of NTDs is still high in some underdeveloped countries and regions. According to reports, folic acid–preventable spina bifida–related stillbirths and child mortalities are still widespread in Ethiopia (Dixon et al., 2019). Recently, our team and international peers have reported that in utero amniotic fluid stem cell therapy techniques may promised new hope to sufferers of NTDs (Abe et al., 2019; Huang et al., 2020). Our research group has made some important improvement in the field of transamniotic bone marrow mesenchymal stem cell transplantation therapy as well (Wei et al., 2020a,b). In this study, we identified a miRNA–mRNA regulatory network including four miRNAs and 39 mRNAs in spina bifida and implicate TSPAN6, YOD1, and KCND3, ubiquitylation pathways, and an antiviral immune response as modulators of neural tube malformations. Following the verification of qRT-PCR and KCND3 was upregulated in spina bifida. KCND3 and its related miR-765 and miR-142-3p appear worthy of further study. However, it remains unclear whether amniocytes reflect the genetic, epigenetic and/or environmental influences related to spina bifida. Such findings and proposed mechanisms obtained from our bioinformatics analyses require validation in spinal cord and adjacent tissues from spina bifida samples by further experimental research using larger sample sizes.
Data Availability Statement
Raw data is available at GEO database under accession number GSE4182 here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE4182. The data used to support the findings of this study are available from the corresponding author upon request.
Ethics Statement
The studies involving human participants were reviewed and approved by Ethics Committee Shengjing Hospital of China Medical University (No. 2015PS264K). Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.
Author Contributions
ZL analyzed the data, wrote the manuscript, and created the figures and tables. As the postdoctoral co-supervisors of ZL, JF, and ZY provided critical revisions and conceptual support. All authors approved of the final submission.
Funding
This work is financially supported by grants from the National Natural Science Foundation of China (Grant Number: 81671469); the National Key Research and Development Program (2016YFC1000505); the National Basic Research Program of China (973 Program, No. 2013CB945402) to ZY.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
Thanks to the teachers Hui Gu and Tianchu Huang and students Henan Zhang and Wanqi Huang in Key Laboratory of Health Ministry for Congenital Malformation for their help.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.583316/full#supplementary-material
Supplementary Figure 1 | Determination of soft-thresholding power in the WGCNA. (A) Analysis of the scale-free index for a set of soft-thresholding powers (β). (B) Analysis of the mean connectivity for a set of soft-thresholding powers. (C) Histogram of connectivity distribution when β = 8. (D) Checking the scale free topology when β = 8.
Supplementary Figure 2 | (A) Screening of modules needed to be merged. (B) Cluster dendrogram of relation between modules and clinical traits. (C) Clustering heatmap of relation between modules and clinical traits.
Supplementary Table 1 | Summary of clinical features of samples in this study.
Supplementary Table 2 | Sequences of primers used for qRT-PCR.
Footnotes
- ^ http://www.ncbi.nlm.nih.gov/geo
- ^ http://www.targetscan.org/vert_72/
- ^ http://mirtarbase.mbc.nctu.edu.tw/
References
Abe, Y., Ochiai, D., Masuda, H., Sato, Y., Otani, T., Fukutake, M.,, et al. (2019). In utero amniotic fluid stem cell therapy protects against myelomeningocele via spinal cord coverage and hepatocyte growth factor secretion. Stem Cells Transl. Med. 8, 1170–1179.
Agarwal, V., Bell, G. W., Nam, J. W., and Bartel, D. P. (2015). Predicting effective microRNA target sites in mammalian mRNAs. eLife 4:e05005.
Alabi, N. B., Thibadeau, J., Wiener, J. S., Conklin, M. J., Dias, M. S., Sawin, K. J.,, et al. (2018). Surgeries and health outcomes among patients with spina bifida. Pediatrics 142:e20173730. doi: 10.1542/peds.2017-3730
An, D., Wei, X., Li, H., Gu, H., Huang, T., Zhao, G.,, et al. (2015). Identification of PCSK9 as a novel serum biomarker for the prenatal diagnosis of neural tube defects using iTRAQ quantitative proteomics. Sci. Rep. 5: 17559.
An, D., Wei, X. W., Zhang, H. N., Liu, D., Ma, W., and Yuan, Z. W. (2020). Spatiotemporal expression of leukemia inhibitory factor receptor protein during neural tube development in embryos with neural tube defects. Neural Regenerat. Res. 15, 705–711. doi: 10.4103/1673-5374.266921
Chou, C. H., Shrestha, S., Yang, C. D., Chang, N. W., Lin, Y. L., Liao, K. W.,, et al. (2018). miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucl. Acids Res. 46, D296-D302.
Danzer, E., Schwarz, U., Wehrli, S., Radu, A., Adzick, N. S., and Flake, A. W. (2005). Retinoic acid induced myelomeningocele in fetal rats: characterization by histopathological analysis and magnetic resonance imaging. Exp. Neurol. 194, 467–475. doi: 10.1016/j.expneurol.2005.03.011
Dennis, G.Jr., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., et al. (2003). DAVID: database for annotation., visualization., and integrated discovery. Genome Biol. 4:P3.
Dixon, M., Kancherla, V., Magana, T., Mulugeta, A., and Oakley, G. P. Jr (2019). High potential for reducing folic acid-preventable spina bifida and anencephaly., and related stillbirth and child mortality., in Ethiopia. Birth Defects Res. 111, 1513–1519.
Duarri, A., Jezierska, J., Fokkens, M., Meijer, M., Schelhaas, H. J., den Dunnen, W. F.,, et al. (2012). Mutations in potassium channel kcnd3 cause spinocerebellar ataxia type 19. Ann. Neurol. 72, 870–880.
Duru, S., Ceylan, S., and Ceylan, S. (2001). Comparative effects of valproic acid sodium for Chiari-like malformation at 9 and 10 days of gestation in the rat. Childs Nerv. Syst. 7, 399–404. doi: 10.1007/s003810000417
Fan, Y., Wang, L., Zhou, F., Zhang, Y., Li, H., Shan, L.,, et al. (2011). Comparative proteomics of spinal cords of rat fetuses with spina bifida aperta. J. Proteom. 75, 668–676. doi: 10.1016/j.jprot.2011.09.006
Gonseth, S., Roy, R., Houseman, E. A., de Smith, A. J., Zhou, M., Lee, S. T.,, et al. (2015). Periconceptional folate consumption is associated with neonatal DNA methylation modifications in neural crest regulatory and cancer development genes. Epigenetics. 10, 1166–1176. doi: 10.1080/15592294.2015.1117889
Gu, H., Li, H., Zhang, L., Luan, H., Huang, T., Wang, L.,, et al. (2012). Diagnostic role of microRNA expression profile in the serum of pregnant women with fetuses with neural tube defects. J. Neurochem. 122, 641–649. doi: 10.1111/j.1471-4159.2012.07812.x
Guix, F. X., Sannerud, R., Berditchevski, F., Arranz, A. M., Horre, K., Snellinx, A.,, et al. (2017). Tetraspanin 6: a pivotal protein of the multiple vesicular body determining exosome release and lysosomal degradation of amyloid precursor protein fragments. Mol. Neurodegen. 12:25.
Hsiao, C. T., Fu, S. J., Liu, Y. T., Lu, Y. H., Zhong, C. Y., Tang, C. Y.,, et al. (2019). Novel SCA19/22-associated KCND3 mutations disrupt human KV 4.3 protein biosynthesis and channel gating. Hum. Mut. 40, 2088–2107.
Huang, J., Wei, X., Ma, W., and Yuan, Z. (2020). The miR-532-3p/Chrdl1 axis regulates the proliferation and migration of amniotic fluid-derived mesenchymal stromal cells. Biochem. Biophys. Res. Commun. 527, 187–193. doi: 10.1016/j.bbrc.2020.04.099
Huin, V., Strubi-Vuillaume, I., Dujardin, K., Brion, M., Delliaux, M., Dellacherie, D.,, et al. (2017). Expanding the phenotype of SCA19/22: parkinsonism., cognitive impairment and epilepsy. Parkinsonism Relat. Disord. 45, 85–89. doi: 10.1016/j.parkreldis.2017.09.014
Ju, Q., Li, M. X., Chen, G., Wang, H. X., Shi, Q. M., Ge, X.,, et al. (2018). Overexpression of YOD1 promotes the migration of human oral keratinocytes by enhancing TGF-beta3 signaling. Biomed. Environ. Sci. 31, 499–506.
Kurihara, M., Ishiura, H., Sasaki, T., Otsuka, J., Hayashi, T., Terao, Y.,, et al. (2018). Novel De Novo KCND3 Mutation in a Japanese Patient with Intellectual Disability., Cerebellar Ataxia., Myoclonus., and Dystonia. Cerebellum. 17, 237–242. doi: 10.1007/s12311-017-0883-4
Lee, Y. C., Durr, A., Majczenko, K., Huang, Y. H., Liu, Y. C., Lien, C. C.,, et al. (2012). Mutations in KCND3 cause spinocerebellar ataxia type 22. Ann. Neurol. 72, 859–869.
Li, H., Gao, F., Ma, L., Jiang, J., Miao, J., Jiang, M.,, et al. (2012). Therapeutic potential of in utero mesenchymal stem cell (MSCs) transplantation in rat foetuses with spina bifida aperta. J. Cell Mol. Med. 16, 1606–1617. doi: 10.1111/j.1582-4934.2011.01470.x
Liu, C., Huang, S., Wang, X., Wen, M., Zheng, J., Wang, W.,, et al. (2019). The Otubain YOD1 Suppresses Aggregation and Activation of the Signaling Adaptor MAVS through Lys63-Linked Deubiquitination. J. Immunol. 202, 2957–2970. doi: 10.4049/jimmunol.1800656
Liu, Y. S., Gu, H., Huang, T. C., Wei, X. W., Ma, W., Liu, D.,, et al. (2020). miR-322 treatment rescues cell apoptosis and neural tube defect formation through silencing NADPH oxidase 4. CNS Neurosci. Ther. 26, 902–912.
Ma, W., Wei, X., Gu, H., Liu, D., Luo, W., An, D.,, et al. (2020). Therapeutic potential of adenovirus-encoding brain-derived neurotrophic factor for spina bifida aperta by intra-amniotic delivery in a rat model. Gene Ther. 110, 2448–2449.
Melo-Filho, A. A., Weber Guimaraes Barreto, M., Capelli Nassr, A. C., Rogerio, F., Langone, F., Pereira, LA, et al. (2009). Corticosteroids reduce glial fibrillary acidic protein expression in response to spinal cord injury in a fetal rat model of dysraphism. Pediatr. Neurosurg. 45, 198–204. doi: 10.1159/000222670
Molloy, A. M., Pangilinan, F., and Brody, L. C. (2017). Genetic risk factors for folate-responsive neural tube defects. Annu. Rev. Nutr. 37, 269–291. doi: 10.1146/annurev-nutr-071714-034235
Nagy, G. R., Gyorffy, B., Galamb, O., Molnar, B., Nagy, B., and Papp, Z. (2006). Use of routinely collected amniotic fluid for whole-genome expression analysis of polygenic disorders. Clin. Chem. 52, 2013–2020. doi: 10.1373/clinchem.2006.074971
Oria, M., Figueira, R. L., Scorletti, F., Sbragia, L., Owens, K., Li, Z.,, et al. (2018). CD200-CD200R imbalance correlates with microglia and pro-inflammatory activation in rat spinal cords exposed to amniotic fluid in retinoic acid-induced spina bifida. Sci. Rep. 8:10638.
Papadopoulos, C., Kirchner, P., Bug, M., Grum, D., Koerver, L., Schulze, N., et al. (2017). VCP/p97 cooperates with YOD1., UBXD1 and PLAA to drive clearance of ruptured lysosomes by autophagy. EMBO J. 36, 135–150. doi: 10.15252/embj.201695148
Paucar, M., Bergendal, A., Gustavsson, P., Nordenskjold, M., Laffita-Mesa, J., Savitcheva, I., et al. (2018). Novel features and abnormal pattern of cerebral glucose metabolism in spinocerebellar ataxia 19. Cerebellum 17, 465–476. doi: 10.1007/s12311-018-0927-4
Qin, P., Li, L., Zhang, D., Liu, Q. L., Chen, X. R., Yang, H. Y.,, et al. (2016). Altered microRNA expression profiles in a rat model of spina bifida. Neural Regen. Res. 11, 502–507. doi: 10.4103/1673-5374.179070
Salas, I. H., Callaerts-Vegh, Z., Arranz, A. M., Guix, F. X., D’Hooge, R., Esteban, J. A.,, et al. (2017). Tetraspanin 6: A novel regulator of hippocampal synaptic transmission and long term plasticity. PLoS One 12:e0171968. doi: 10.1371/journal.pone.0171968
Schimmack, G., Schorpp, K., Kutzner, K., Gehring, T., Brenke, J. K., Hadian, K.,, et al. (2017). YOD1/TRAF6 association balances p62-dependent IL-1 signaling to NF-kappaB. eLife 6:e22416.
Sehrawat, S., Koenig, P. A., Kirak, O., Schlieker, C., Fankhauser, M., and Ploegh, H. L. (2013). A catalytically inactive mutant of the deubiquitylase YOD-1 enhances antigen cross-presentation. Blood 121, 1145–1156. doi: 10.1182/blood-2012-08-447409
Shan, L., Fan, Y., Li, H., Liu, W., Gu, H., Zhou, F., et al. (2012). Proteomic analysis of amniotic fluid of pregnant rats with spina bifida aperta. J. Prot. 75, 1181–1189. doi: 10.1016/j.jprot.2011.10.033
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D.,, et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303
Smith, M. T., Wissinger, J. P., Smith, C. G., and Huntington, H. W. (1978). Experimental dysraphism in the rat. J. Neurosurg. 49, 725–729. doi: 10.3171/jns.1978.49.5.0725
Tanji, K., Mori, F., Miki, Y., Utsumi, J., Sasaki, H., Kakita, A.,, et al. (2018). YOD1 attenuates neurogenic proteotoxicity through its deubiquitinating activity. Neurobiol. Dis. 112, 14–23. doi: 10.1016/j.nbd.2018.01.006
Turner, C. G., Pennington, E. C., Gray, F. L., Ahmed, A., Teng, Y. D., and Fauza, D. O. (2013). Intra-amniotic delivery of amniotic-derived neural stem cells in a syngeneic model of spina bifida. Fetal Diagn. Ther. 34, 38–43. doi: 10.1159/000350267
Vincent, A. K., Noor, A., Janson, A., Minassian, B. A., Ayub, M., Vincent, J. B., et al. (2012). Identification of genomic deletions spanning the PCDH19 gene in two unrelated girls with intellectual disability and seizures. Clin. Genet. 82, 540–545. doi: 10.1111/j.1399-0004.2011.01812.x
Wang, J., Wen, Y., Zhang, Q., Yu, S., Chen, Y., Wu, X.,, et al. (2019). Gene mutational analysis in a cohort of Chinese children with unexplained epilepsy: identification of a new KCND3 phenotype and novel genes causing Dravet syndrome. Seizure 66, 26–30. doi: 10.1016/j.seizure.2019.01.025
Wang, Y., Tong, X., Omoregie, E. S., Liu, W., Meng, S., and Ye, X. (2012). Tetraspanin 6 (TSPAN6) negatively regulates retinoic acid-inducible gene I-like receptor-mediated immune signaling in a ubiquitination-dependent manner. J. Biol. Chem. 287, 34626–34634. doi: 10.1074/jbc.m112.390401
Wei, X., Cao, S., Ma, W., Zhang, C., Gu, H., Liu, D.,, et al. (2020a). Intra-Amniotic Delivery of CRMP4 siRNA improves mesenchymal stem cell therapy in a rat spina bifida model. Mol. Ther. Nucl. Acids 20, 502–517. doi: 10.1016/j.omtn.2020.03.007
Wei, X., Li, H., Miao, J., Liu, B., Zhan, Y., Wu, D., et al. (2013). miR-9∗- and miR-124a-Mediated switching of chromatin remodeling complexes is altered in rat spina bifida aperta. Neurochem. Res. 38, 1605–1615. doi: 10.1007/s11064-013-1062-8
Wei, X., Ma, W., Gu, H., Liu, D., Luo, W., Bai, Y.,, et al. (2020b). Transamniotic mesenchymal stem cell therapy for neural tube defects preserves neural function through lesion-specific engraftment and regeneration. Cell Death Dis. 11:523.
Xue, J., Gu, H., Liu, D., Ma, W., Wei, X., Zhao, L.,, et al. (2018). Mitochondrial dysfunction is implicated in retinoic acid-induced spina bifida aperta in rat fetuses. Int. J. Dev. Neurosci. 68, 39–44. doi: 10.1016/j.ijdevneu.2018.04.003
Zhang, H., Guo, Y., Gu, H., Wei, X., Ma, W., Liu, D., et al. (2019). TRIM4 is associated with neural tube defects based on genome-wide DNA methylation analysis. Clin. Epigenet. 11:17.
Zhang, H. N., Guo, Y., Ma, W., Xue, J., Wang, W. L., and Yuan, Z. W. (2019). MGMT is down-regulated independently of promoter DNA methylation in rats with all-trans retinoic acid-induced spina bifida aperta. Neural Regen. Res. 14, 361–368. doi: 10.4103/1673-5374.244799
Zhao, J. J., Sun, D. G., Wang, J., Liu, S. R., Zhang, C. Y., Zhu, M. X., et al. (2008). Retinoic acid downregulates microRNAs to induce abnormal development of spinal cord in spina bifida rat model. Childs Nerv. Syst. 24, 485–492. doi: 10.1007/s00381-007-0520-5
Keywords: spina bifida, weighted gene co-expression network analysis, bioinformatics analysis, pathological process, hub genes
Citation: Li Z, Feng J and Yuan Z (2020) Key Modules and Hub Genes Identified by Coexpression Network Analysis for Revealing Novel Biomarkers for Spina Bifida. Front. Genet. 11:583316. doi: 10.3389/fgene.2020.583316
Received: 14 July 2020; Accepted: 09 November 2020;
Published: 02 December 2020.
Edited by:
Qihua Fu, Shanghai Children’s Medical Center, ChinaReviewed by:
David D. Eisenstat, University of Alberta, CanadaBo Wang, Shanghai Children’s Medical Center, China
Copyright © 2020 Li, Feng and Yuan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Juan Feng, anVhbmZlbmdAY211LmVkdS5jbg==; Zhengwei Yuan, eXVhbnp3QGhvdG1haWwuY29t
†These authors have contributed equally to this work