Genetic regulatory subnetworks and key regulating genes in rat hippocampus perturbed by prenatal malnutrition: implications for major brain disorders

Demographic characteristics of RLP50 and neonatal rats

Prenatal malnutrition did not significantly alter birth numbers. However, the neonatal birth weights in both the LP and RLP50 groups were lower than those in the control group, with the RLP50 group having the lowest birth weights (p < 0.001). Meanwhile, the RLP50 group showed significant maternal weight percent gain compared to the other two groups (p < 0.001). As previously reported in our studies [12, 13], we observed a tendency for the LP group to build a smaller nest and a complete disappearance of nest building in the RLP50 group after gestational exposure to simulated famine (p < 0.001), indicating that gestational exposure to malnutrition and the stress of starvation resulted in maternal behavior disruption. Behavior tests of neonatal rats between groups did not find significant differences. The sole purpose of this low protein group was to provide a basis to calculate the appropriate diet for the "RLP50" group, (Figure 1) which was used to simulate the prenatal famine situation, thus only rats in RLP50 and control groups were selected to perform the gene expression profiling analysis.

Figure 1. Flow chart of the analytic procedure.

DEGs in the hippocampus and prefrontal cortex

After microarray data pre-processing, final data sets included 17 hippocampal tissues of rats (10 prenatal nutrition deficiency offspring groups and 7 normal offspring groups) and 17 prefrontal cortex tissues of rats (8 prenatal nutrition deficiency offspring groups and 9 normal offspring groups). The distribution of gender is not significantly different between the RLP50 group and normal nourished offspring groups (Table 1).

We analyzed rat hippocampus tissue including 10 prenatal nutrition deficiency offspring individuals and 7 normal offspring individuals, as well as prefrontal cortex including 7 prenatal nutrition deficiency offspring individuals and 9 normal offspring individuals, respectively. We found that 1844 genes were differentially expressed between prenatal nutrition deficiency and control offspring individuals, 209 of which were significantly different after FDR adjustment (FDR-corrected P < 0.05) (Supplementary Figure 2, Supplementary Table 3) in the hippocampus. In the prefrontal cortex, 717 genes were differentially expressed between prenatal nutrition deficiency and control offspring individuals, while, we did not find any differentially expressed genes with significance after FDR adjustment (FDR-corrected *P < 0.05) (Supplementary Table 3), which is consistent with our previous report [12] and suggested the hippocampus is more sensitive to the exposure to prenatal nutrition deficiency. In this current study, we identified a much smaller number of DEGs compared with the previous one, in which 2987 DEGs in the hippocampus and 415 DEGs in the prefrontal cortex were reported, because there was no correction for multiple test performed in that previous study.

Co-expression modules related to prenatal nutritional status

We conducted a weighted gene co-expression network analysis (WGCNA). Before using the WGCNA R package to construct the network, we used cluster analysis to check the quality of hippocampal and prefrontal cortex samples further and no samples were removed (Supplementary Figure 3). Here, the power of β equaling 15 was selected as the soft threshold for constructing a scale-free network in the hippocampus (Supplementary Figure 4A, 4B), while the power of β equaling 18 is selected for prefrontal cortex (Supplementary Figure 4E, 4F).

In all 17 subjects with hippocampal expression profiling data, 16 modules were identified, 3 of which (blue, pink and salmon modules) showed significant association “Factor (e.g. prenatal nutritional status)” in overall expression patterns (represented by MEs) (Figure 2A), while no module was significantly associated with “Gender” (Figure 2B). In all 16 subjects with prefrontal cortex expression profiling data, 12 modules with co-expressed genes were identified, but no module was significantly correlated with “Factor” (Figure 2C) or “Gender” (Figure 2D).

Figure 2. Association test of modules with prenatal nutritional status and gender (FDR-corrected *P < 0.05). (A) Module-level differential expression related to prenatal nutritional status in the hippocampus. (B) Module-level differential expression related to prenatal nutritional status in the prefrontal cortex. (C) Module-level differential expression related to gender status in the hippocampus. (D) Module-level differential expression related to gender status in the prefrontal cortex.

Module preservation test

We performed preservation analysis of the expression profiles of the hippocampus and prefrontal cortex and found that there were some weak preserved modules between these two different brain tissues (Supplementary Figure 5), and the findings indicated by medianRank and Zsummary statistics were consistent, which indicated that module size has little effect on preservation analysis. Six hippocampal modules were highly preserved with the prefrontal cortex, including the blue module in the hippocampus, seven modules, including a salmon module in the hippocampus, were moderately preserved in the prefrontal cortex, and two modules, including the pink module in the hippocampus, were weakly preserved with the prefrontal cortex. These findings demonstrated that the gene expression patterns between the hippocampus and prefrontal cortex are different to a large extent, and prenatal malnutrition may easily trigger the changes of gene expression patterns in the hippocampus.

Enrichment analysis of interesting modules

We conducted gene ontology enrichment analyses of the three modules in the hippocampus significantly related to prenatal nutrition deficiency and found that the module genes related to prenatal nutrition deficiency were involved in neuronal development and function. In the blue module, we found most genes are related to synaptic development, as their cellular component was associated with postsynaptic specialization, neuron to neuron synapse, postsynaptic density, asymmetric synapse, and synaptic membrane. And their genetic functions are related to learning and social behavior (Figure 3A, Supplementary Table 5). In the pink module, we found the cellular components of most genes were related to the synaptic membrane, postsynaptic membrane, and neuron to neuron synapse, etc. Meanwhile, their biological function was involved in learning, memory, and cognition (Figure 3B, Supplementary Table 5). In the salmon module, cellular components of most genes were regarding neuron projection terminus, synaptic vesicle, and neuron projection terminus, etc. Their biological function was major related to synaptic vesicles, such as synaptic vesicle cycle, synaptic vesicle transport, the establishment of synaptic vesicle localization, and synaptic vesicle exocytosis (Figure 3C, Supplementary Table 5). In summary, GO functional enrichment analysis demonstrated that prenatal nutrition deficiency affected the transcriptional expression of neuronal development-related genes.

Figure 3. Gene Ontology and KEGG pathway enrichment analysis of the modules significantly related to prenatal nutrition deficiency in the hippocampus. (A) Top 5 terms significantly enriched in the blue module; (B) Top 5 terms significantly enriched in the pink module; (C) Top 5 terms significantly enriched in the salmon module.

KEGG enrichment analysis demonstrated that the genes within the blue, pink, and salmon modules were involved in the ribosome, circadian entrainment and synaptic vesicle cycle, respectively (Figure 3, Supplementary Table 6). This further indicated that the salmon module genes are more closely related to synaptic development.

The rat brain hippocampus cell-types enrichment analysis showed that the salmon module was enriched with genes specifically expressed in pyramidal neurons, which are the most common excitable neurons (Figure 4, Supplementary Table 7) sending and receiving nerve impulses within the cerebral cortex, hippocampus, and the amygdala. In mammals, it is thought that pyramidal neurons play a key role in cognitive functions, such as perception, reasoning, remembering, thinking and understanding. We did not find significant enrichment of cell-types in blue and pink modules.

Figure 4. Cell-type enrichment analysis of the blue, pink and salmon modules significantly responding to prenatal nutrition deficiency in the hippocampus.

To further evaluate the effects of prenatal nutrition deficiency on the risk of neurodevelopmental diseases, we conducted disease-related gene set enrichment analysis based on the DisGeNET database, in which 24166 human diseases and their related-genes were included. We concentrated on 10 major brain diseases, including attention deficit hyperactivity disorder (ADHD), autism, anoxia, bipolar disorder, Parkinson’s disease, major depressive disorder, Alzheimer’s disease, schizophrenia, and glioma. The results demonstrated that the blue module was significantly overlapped with genes related to Alzheimer's disease, autism, and schizophrenia, while the pink and salmon modules were significantly overlapped with schizophrenia-related genes. These findings suggest that prenatal nutrition deficiency may increase the risk of neurodevelopmental diseases, including Alzheimer's disease, autism, and schizophrenia (Table 2).

Identification of hub genes

As shown in Figure 5A, 5D and 5G, the Factor-associated down-regulated blue, pink and salmon modules were significant positive correlations between GS and MM in the hippocampus. This showed that these three modules have high quality. We selected the top 50 genes in each of the three modules (Figure 5B, 5E, and 5H) with high GS and MM from blue, pink and salmon modules, (Supplementary Table 4) respectively. And we also showed the top 50 high connectivity genes based on the co-expression network, which may provide a basis for other studies (Figure 7A–7C, Supplementary Figure 6A–6C). We constructed PPI networks of these three module genes using STRING and highlighted the top 50 high connectivity genes (Figure 7D–7F, Supplementary Figure 6D–6F). Hub genes selection was based on features including significantly different between prenatal nutrition deficiency and control offspring, high GS and MM and connection degree ≥ 5 in the PPI network (Figure 5C, 5F, and 5I). We found twelve hub genes in the blue module, including Hmgn1, Ssbp1, LOC684988, Rpl23, Gga1, Rhobtb2, Dhcr24, Atg9a, Dlgap3, Grm5, Scn2b and Furin, two hub genes in the pink module, including Sh3kbp1 and Ubqln1, and one hub gene (Unc13a) in the salmon module (Supplementary Figure 7A).

Figure 5. Identification of final hub genes from the hippocampus tissue. Scatter plotting of the correlation between gene significance (GS) and module membership (MM), top 50 module genes with high GS and MM, and Venn diagram of the overlapping genes belonging to the PPI network with degree ≥ 5, top 50 module genes with high GS and MM, and the differentially expressed genes (DEGs) in the blue (A–C), pink(D–F), and salmon (G–I) modules.

qPCR validation of the hub genes

To verify the main conclusion drawn from the microarray results, the relative expression levels of the 15 key genes, including Hmgn1, Ssbp1, LOC684988, Rpl23, Gga1, Rhobtb2, Dhcr24, Atg9a, Dlgap3, Grm5, Scn2b, Furin, Sh3kbp1, Ubqln1, and Unc13a, were determined using qPCR. No significant difference in body weight and gender were detected between the two groups of rats. The qPCR analysis results indicated that the expression levels of Hmgn1, Ssbp1, LOC684988, Rpl23, Gga1, Rhobtb2, Dhcr24, Dlgap3, Scn2b, Furin, Ubqln1, and Unc13a were significantly up-regulated in samples from RLP50 group compared with the prenatal normal nourished group, whereas Grm5 and Sh3kbp1 were significantly down-regulated. Meanwhile, there were no significant differences in the expression level of Atg9a between the two groups.

Inference of miRNA–mRNA-transcription factor interactions

The final “real” 15 hub genes were simultaneously input into the miRDB database to identify the microRNAs target sites, and were combined with DiRE database to identify transcription factors regulating these hub genes (Figure 6, Supplementary Table 10). We found a transcription factor, Cebpa, was differentially expressed with statistical significance (Supplementary Figure 7B). Cebpa may transcribe the hub gene Rhobtt2 in the blue module, whereas miR-3569 and miR-18a-3p may suppress the translation of Rhobtt2. Meanwhile, miR-18a-3p may also suppress the translation of hub gene Furin in the blue module, which may cause the indirect correlated expression between Rhobtt2 and Furin in the blue module. All hub genes are co-regulated through miRNAs, or transcription factors, or other potential novel regulators of transcription and post-transcription dominating prenatal nutritional deficiency-induced gene expression changes. For example, miR-let-7g-3p can regulate both Hmgn1 and Dlgap3; meanwhile the transcription factor NXC may regulate both Rpl23 and Gga1 in the miRNA–mRNA-transcription factor interaction network. Detailed illustration of the inferred miRNA–mRNA-transcription factor interactions for the 15 hub genes is present in Figure 6.

Figure 6. Illustration of miRNA–mRNA-transcription factor interaction for hub genes identified in the blue, pink and salmon modules in the hippocampus. Each hub gene is denoted as solid circle; each transcription factor is denoted as solid square; and each miRNA is denoted as a solid delta.

In the study of hub genes, we additionally screened the Top50 connectivity gene (Figure 7A–7C, Supplementary Figure 6A–6C, Supplementary Tables 8 and 9) in the co-expression network and the Top50 connectivity gene (Figure 7D–7F and Supplementary Figure 6D–6F) in the PPI network. These contents are not highly related to the main ideas of this article, but could provide some basis for the follow-up study.”

Figure 7. Top 50 genes with high network connectivity. (A–C) Top 50 genes determined through co-expression networks in the blue, pink, and salmon modules, respectively; (D–F) Top 50 genes determined through PPI networks in the blue, pink, and salmon modules, respectively.

Animal model

Sprague-Dawley rats (Shanghai Slack Laboratory Animal Co. Ltd.), were used in all experiments under the Guide for the Care and Use of Laboratory Animals (Washington, DC: National Academic Press). All experimental procedures and protocols were approved by the Institutional Animal Care and Use Committee at the Shanghai Jiao Tong University.

Seventeen Sprague–Dawley female rats were randomly separated into three groups on the first day after mating, and were fed their respective diets until E18 (embryo, 18 days), two days before parturition (control: 5 dams; malnourished group [LP]: 5 dams; famine group [RLP50]: 7 dams). The control group was given a standard rodent diet (20% protein, Research Diets, Inc. D12450B) and water ad libitum and the LP group was given a low-protein (6% protein, Research Diets, Inc. D06022301) diet and water ad libitum. The RLP50 group was given 50% of the LP group’s low-protein diet, reflecting both the protein malnutrition and food-deficiency likely to prevail during the famine. All the female rats were fed with the standard rodent diet after E18. After parturition, all litters were culled to 8 pups, which were fostered by their own mothers. On day 21, all pups were weaned and placed on a standard rodent diet. Rats were raised to be 10 weeks old and were anesthetized before being killed. The PFC and hippocampus from the adult offspring of both control (n = 60) and RLP50 (n = 60) rats were then isolated and stored in RNAlater® Solution Life technologies). The treatment of the pregnant rats in RLP50 was according to our previous model [12, 13]

Microarray data pre-processing and identification of DEGs in the hippocampus and prefrontal cortex

Total RNA was extracted from the hippocampus and prefrontal cortex tissues derived from 10 control and 10 RLP50 rats using mirVana™ miRNA Isolation Kit. The RNA integrity number (RIN) was detected by Agilent Bioanalyzer and the samples with RIN > 7.5 were selected for microarray gene expression analysis using NimbleGen Rat Gene Expression 12 x 135K Arrays. (Supplementary Table 1) Quantile normalization of microarray data was performed to obtain standardized using preprocessCore package in R [38] (Supplementary Figure 1A, 1B, 1H, 1G and Supplementary Table 2). A total of 26419 microarray probes were converted into 9532 Ensembl gene IDs by using the biomaRt package [39]. Samples defined as outliers based on standardized sample network connectivity Z-scores < -2 were removed [40] (Supplementary Figure 1C and 1I). The batch effect of expression profiles was treated by using the ComBat function of the sva package in R (Supplementary Figure 1D, 1E, 1J, 1K, and Supplementary Table 2). We balanced maternal prenatal nutrition deficiency/control status across gender and there was no significant difference in gender between the two groups (p-value > 0.05) (Supplementary Figure 1F and 1I). Differential expression analysis using a linear empirical Bayes model based on the limma package [41]. Differentially expressed genes with statistical significance were identified after p-value correction with a false discovery rate (FDR).

Weighted gene co-expression network construction

To assess the inter-correlation of the intensities of the 9532 Ensembl gene IDs (that remained after preprocessing) and the relationship between them and feature data, signed co-expression analysis was performed using the WGCNA R package. The pair-wise correlation matrix was calculated by biweight midcorrelation, and an adjacency matrix is calculated by raising the correlation matrix to power [42]. Based on a fit to scale-free topology (R2 > 0.8) [40] (Supplementary Figure 4D and 4H), the powers of 15 (Supplementary Figure 4A and 4B) and 18 (Supplementary Figure 4E and 4F) were used for the construction of hippocampus and prefrontal cortex networks, respectively. Next, the dynamic tree cut methods were used to construct the basic network and the topological overlap dissimilarity matrix (1-TOM) was then used as input average linkage hierarchical clustering to create the network dendrogram [43] (Supplementary Figure 4C and 4G). For the basic network, we set the minimum number of modules to be 100 and the deep split parameter in 2 for a medium sensitivity of cluster splitting. Finally, modules were defined as branches of the dendrogram, which was cut based on the hybrid dynamic tree-cutting method [40, 42, 43].

Modules were summarized by their first principal component (ME, module eigengene) and the mergeCloseModules function in the WGCNA R package was used to merged modules whose eigengenes were more than 0.85 correlated [40]. We define module membership of each gene in each module (kME) by correlating its gene expression values with the module eigengene of a given module [44]. Pearson correlation values between MEs and “Factor” (famine offspring vs control offspring) or “Gender” were calculated, and p-values were adjusted with FDR. For each module, the correlation of two parameters [the gene significance (GS) and the module membership (MM)] was evaluated. GS stood for the magnitude of the correlation between the intensities of individual expression genes in the module and “Factor”, and MM meant the magnitude of the correlation between the intensities of individual expression genes in the module and MEs of the module.

Module preservation test between the hippocampus and prefrontal cortex

To examine whether the co-expression structure (density and connectivity) is similar between the hippocampus and prefrontal cortex samples, we conducted module preservation analysis based on the modulePreservation function in the WGCNA R package. We used a composite preservation statistics method to analyze the conservativeness of the module [45]. For each module, the Z-summary statistic was used to measure module density and intramodular connectivity metrics. In correlation networks, the 4 density preservation statistics are summarized by Zdensity (function 1), the 3 connectivity based statistics were summarized by Zconnectivity (function 2), and for each module, the Zsummary (function 3) was used to measure that combines module density and intramodular connectivity metrics defined as follows:

Z_density = median (Z_meanCor, Z_meanAdj, Z_propVarExpl, Z_meanKME) (function 1)

Z_connectivity = median (Z_cor.kIM, Z_cor.kME, Z_cor.cor) (function 2)

Z_summary = (Z_density + Z_connectivity)/2 (function 3)

We usually defined: if Z_summary > 10 indicated high preservation among modules; if 2 < Z_summary < 10 indicated weak to moderate preservation; if Z_summary < 2 indicated no evidence that the module preserved.

The module size has a great influence on Z statistics. Therefore, when comparing the preservation statistics of different sized modules, we need to conduct medianRank for preservation analysis. It shows that modules with lower median rank tend to show stronger preservation statistics than those with higher median rank.

Functional enrichment analysis of modules associated with prenatal nutritional status

Gene ontology (GO) terms, including biological processes, molecular functions and cellular components, and KEGG pathway enrichment analyses in the modules, were conducted with the R package, clusterProfiler (version 3.5) [45]. We used the annotation function in clusterProfiler to group significant GO terms and pathways based on co-associated genes to remove redundant terms. GO terms and pathway groups were considered significant when the FDR-corrected P-value of the clustered GO terms and pathways is below 0.05.

Cell-type enrichment analysis for interesting modules was conducted using the R package anRichment [44]. For enrichment analysis of cell-type-specific genes, we used cell type-specific genes identified in gene expression datasets from neurons, astrocytes, myelinating oligodendrocytes, microglia, and endothelial cells [46]. Genes were considered cell-type specific, if their expression levels had more than ten-fold higher compared to the mean levels of expression in the other cell types. The statistical significance of enrichment was determined by Fisher’s exact test.

The R package GeneOverlap(version 1.22.0) [47] was employed for overlap analysis of the module genes with disease-related genes annotated in a database named DisGeNET (http://www.disgenet.org) to evaluate if the modules are enriched with neuropsychiatric genes. This package uses Fisher's exact test to calculate the significance of each pair of gene lists in comparison to the genomic background. The function returns the number of intersecting genes between the two lists, the P-value, and the estimated odds ratio. The null hypothesis is that the odds ratio is no larger than 1, and values larger than this indicate a positive association between two lists. Therefore, if a disease gene set in DisGeNET is significantly overlapped with the gene set of a module responding to prenatal nutritional deficiency, it indicates that this disease is closely related to hippocampus dysregulation responding to prenatal malnutrition. All analyses with an FDR adjusted P-value of < 0.05 were selected as statistical significance.

Protein-protein interaction network construction and identification of hub genes in interesting modules

The Protein-Protein Interaction (PPI) network of the genes in the selected co-expression module was constructed according to the STRING v11.0 database [48] and was visualized using Cytoscape.

Hub genes have high connectivity within a gene module and are functionally significant [44]. In our study, after an interesting module was chosen, hub genes were defined by top50 module genes with high gene significance (GS) and module membership (MM). In the PPI network, genes with a combined score of ≥0.6 and a connectivity degree of ≥5 were also defined as hub genes. The shared hub genes in both the co-expression network and PPI network were regarded as “real” hub genes for further analyses. A Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) was constructed to identify these “real” hub genes [49].

Verification of hub genes with qPCR based on independent samples

Hub genes were selected for validation using quantitative real-time polymerase chain reaction (qPCR) in the different samples including 40 prenatal malnourished rats(RLP50 group) and 40 prenatal normal nourished rats(control group), as used for the microarray experiments. For the quantification of mRNAs, 1 μg of the total RNA was converted into cDNA using a reverse transcription kit (Promega) following the manufacturer's instruction. The TaqMan® Gene Expression Assay was then used to detect the gene expression of the hub genes. The qPCR was conducted using SYBR Premix Ex Taq (Takara, Japan) in 20 μl reaction solution containing 10 μl SYBR Premix Ex Taq (2X), 0.4 μl forward primer and 0.4 μl reverse primer, 0.4 μl ROX reference dye, 1μl cDNA and 8.2 μl ddH2O. The PCR amplification procedure was carried out at 95°C for 10 s; 40 cycles of 95°C for 5 s and 60°C for 34 s; followed by disassociation curve analysis in an ABI 7500 fast real-time PCR system (Applied Biosystems, USA). The amplification reaction without the template was used as a no template control. All reactions were performed in triplicate. The β-actin gene was used as an internal reference. Primers used for gene amplification are available on request. The relative mRNA abundance was calculated using the 2-ΔΔCt method. Statistical comparison of the levels was analyzed using two-tail unpaired Student's t-test, and differences were considered significant if P < 0.05.

Modeling miRNA–mRNA-transcription factor interactions for hub genes

To explore transcription factors (TFs) and miRNA–mRNA interaction in particular modules and to predict the possible regulatory relationships between them that dominate prenatal nutritional deficiency-induced gene expression changes, we used two common online databases, respectively. The TFs that bind onto hub genes were screened using DiRE [50] (https://dire.dcode.org), a web server can predict common REs for any set of input genes. microRNAs of hub genes were screened using miRDB [51] (http://mirdb.org/mining.html). This database collected miRNA-genes regulatory relationships.