All animal protocols were approved by the Zhejiang University Animal Care and Use Committee. All experiments were performed with Institute of Cancer Research (ICR) mice. Eight-week-old virgin ICR females (n = 40) were mated with normal ICR males. Pregnancy was dated by the presence of a vaginal plug (designated day 0 [D0] of pregnancy). After a 12-h fast, the D0 females were randomly divided into a control group (Ctrl) and a GDM group. Mice in the GDM group were injected with a single intraperitoneal injection of streptozotocin (STZ; Sigma, St. Louis, MO) in 0.1 mmol/L citrate buffer (pH 4.5) at a dose of 150 mg/kg body weight. Control pregnant mice received an equal volume of citrate buffer. Blood glucose levels were measured at day 3, day 7 and day 18.5 of pregnancy via the tail vein, and diabetes was defined as a blood glucose level between 14 and 19 mmol/L as previously described (Ding et al., 2012). At embryonic day 18.5 (E18.5), the fetal mice were removed from 9 dams per group by abdominal cesarean section after immediate euthanization. Then, the meninges were carefully removed, and the hippocampi were dissected for electron microscopy (n = 4 per group) and analysis of the methylome and transcriptome (n = 3 per group) (Supplementary Figures S1A,C,D).

To further mimic the hyperglycemic environment in the middle and late trimesters of most clinical GDM, on day 6 and day 12 of pregnancy, females fasting for 8 h received an intraperitoneal injection of STZ as the late GDM group (n = 15). Pregnant females in the Ctrl group (n = 15) received an equal volume of citrate buffer. Blood glucose levels were measured via the tail vein, and diabetes was defined as a blood glucose level between 14 and 19 mmol/L (Zhu et al., 2019). At E18.5, the fetal mice brains were isolated from the dams by abdominal cesarean section after immediate euthanization. Then, untargeted metabolomics was performed on the whole brains of E18.5 fetal mice from Ctrl group (n = 7) and GDM group (n = 10) (Supplementary Figures S1B,D).

Samples (n = 4 per group) were fixed in 2.5% glutaraldehyde and postfixed with 3% osmium tetroxide (OsO4) for 2 h. The samples were dehydrated in a graded series of ethanol, embedded in Epon resin, sliced into 100 nm pieces and finally imaged with a transmission electron microscope at 80 kV (Hitachi, Japan, H7650 TEM).

Reduced representation bisulfite sequencing (RRBS) was performed in hippocampi obtained at embryonic day 18.5 from the Ctrl and GDM group (n = 3 per group) (Genergy Biotechnology Co., Ltd., Shanghai, China). Briefly, 5 μg genomic DNA was digested using the methylation-insensitive restriction enzyme MspI (New England Biolabs, Beverly, MA, United States). A Qiagen Mini Purification kit (Qiagen, Hilden, Germany) was used to purify the digested products. Then, the ends of each restriction fragment were filled in, and adenosine was added at the 3′-end. Methylated paired-end Illumina adapters were ligated to the ends of the DNA fragments using T4 DNA ligase, and fragments sized 100–200 bp were purified by agarose gel extraction. The purified fragments were treated with sodium bisulfite and then amplified by PCR. The final PCR products were sequenced on HiSeq 3,000 (Illumina Inc., San Diego, CA, United States).

Hippocampi from the Ctrl and GDM groups (n = 3 per group) were isolated from fetal mice and snap frozen in liquid nitrogen (for RNA-seq). Total RNA was isolated from hippocampi using Trizol. The purity of RNA was determined using a NanoDrop (Thermo, United States) spectrophotometer (OD260/280 is 1.8–2.2), and the integrity of RNA was assessed by electrophoresis (Tianneng, China) using agarose gel (28S:18S > 1.5). mRNA was enriched by magnetic beads with oligo (dT) and fragmented into 200 bp. RNA-seq libraries underwent 15 cycles of PCR amplification to obtain cDNA library. After RNA library construction, Qubit (Invitrogen, United States) was used to detect accurate RNA concentrations (≥500 ng/μl). Paired-end sequencing with a read length of 151 bp was conducted using the Illumina HiSeq3000 (Illumina Inc., San Diego, CA, United States) platform.

100 mg (±2%) of each sample was weighed and transferred into a 5 ml EP tube, which was prepared with 2 ml tissue extract (75% 9:1 methanol: chloroform, 25% ddH2O) (−20°C) and 3 steel beads, and then grinded for 60 s at 50 Hz in a high-throughput tissue grinder twice. Ultrasound-assisted extraction was processed at room temperature for 30 min and then the samples were put on ice for 30 min. After that, samples were centrifuged at 14,000 rpm for 20 min at 4°C and 1250 μl of the supernatant from each sample were transferred into another 2 ml centrifuge tube. Then, samples were concentrated to dry in vacuum and dissolved with 200 μl 2-chlorobenzalanine (4 ppm) 50% acetonitrile solution and the supernatant was filtered through 0.22 µm membrane to obtain the prepared samples for LC-MS. Every 20 ul of each sample was mixed as a quality control (QC) sample. The rest of the samples were used for LC-MS detection.

Chromatographic separation was accomplished in a Thermo Ultimate 3,000 system equipped with an ACQUITY UPLC^® HSS T3 (150 × 2.1 mm, 1.8 μm, Waters) column maintained at 40°C. The temperature of the autosampler was 8°C. Gradient elution of analytes was carried out with 0.1% formic acid in water (C) and 0.1% formic acid in acetonitrile (D) or 5 mM ammonium formate in water (A) and acetonitrile (B) at a flow rate of 0.25 ml/min. Injection of 2 μl of each sample was performed after equilibration. An increasing linear gradient of solvent B (v/v) was used as follows: 0–1 min, 2% B/D; 1–9 min, 2–50% B/D; 9–12 min, 50–98% B/D; 12–13.5 min, 98% B/D; 13.5–14 min, 98–2% B/D; 14–20 min, 2% D-positive model (14–17 min, 2% B-negative model).

The ESI-MSn experiments were executed on a Thermo Q Exactive mass spectrometer with spray voltages of 3.8 kV and −2.5 kV in positive and negative modes, respectively. Sheath gas and auxiliary gas were set at 30 and 10 arbitrary units, respectively. The capillary temperature was 325°C. The analyzer scanned over a mass range of m/z 81–1,000 for full scan at a mass resolution of 70,000. Data-dependent acquisition (DDA) MS/MS experiments were performed with HCD scans. The normalized collision energy was 30 eV. Dynamic exclusion was implemented to remove some unnecessary information in MS/MS spectra.

RRBS data were analyzed using the R package minfi for raw data preprocessing (quality control and normalization). RRBS reads were mapped to the reference mouse genome (mm10) by Bismark (version 0.16.3). DSS V2.30.1 in R was used to detect DMLs. Differentially methylated loci (DMLs) were analyzed based on a Bayesian approach (Feng et al., 2014) and were summarized as follows: two groups were modeled according to the Bayesian stratification model, and the Wald test was applied to each locus to obtain a p-value for each CpG site. For each CpG site, a difference in methylation value between two groups ≥5% and a posteriori probability of Wald test ≥0.95, or p-value < .01 was considered to be a differentially methylated loci (DML). A methylation region was defined as a differently methylated regions (DMR) when it met the following three criteria: 1) the length of this region was at least 50 bp; 2) the region contained no less than three CpG sites; and 3) the proportion of DMLs in this region was no less than 50%. When a DMR showed no less than 50% overlap with one element of the gene, it was defined as a differentially methylated gene (DMG).

Quality control of the raw and trimmed reads of RNA-seq was performed using FastQC (v0.11.5). The STAR (2.5.3a) software package performed the task that mapped high-throughput sequencing reads to a reference genome in RNA-seq data analysis (Dobin et al., 2013). The DEseq2 package v1.16.1 in R was used to identify the differentially expressed genes (DEGs) in GDM hippocampi compared with the Ctrl group (Love et al., 2014). The t-test and Benjamini–Hochberg method were used to calculate the p-value and FDR, respectively. By using a stringent threshold and significant criteria of FDR <0.05 and >1.2-fold change or <0.833-fold change, a total of 953 DEGs were screened out.

Metabolites were identified according to their exact molecular weight and the MS/MS fragmentation pattern by a comparison with those in the online Metlin database (http://metlin.scripps.edu), MoNA (https://mona.fiehnlab.ucdavis.edu//). The mass error was 15 ppm. Score plots of orthogonal partial least-squares discriminant analysis (OPLS-DA) clearly separated the GDM group from the control group. The variable importance in projection (VIP) generated in OPLS-DA processing represented the contribution to the discrimination of each metabolite ion between groups. The t-test was used to calculate the p-value. The significantly differential metabolites were identified using the threshold of p-value < .05 and VIP ≥1. MetaboAnalyst 4.0 (www.metaboanalyst.ca) was used for metabolic pathway analysis based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Enrichment for each identified pathway was calculated using Fisher’s exact test and p-value < .05 was considered statistically significant.

ClusterProfiler V3.8.0 in R was used to identify and visualize the gene ontology (GO) terms and KEGG pathways enriched by the DMGs and shared genes of DEGs and DMGs. To identify the function of the aforementioned DEGs, KEGG analyses were performed by Cytoscape ClueGO (v2.5.7) and CluePedia (v1.5.7). Gene enrichment for each identified pathway was calculated using Fisher’s exact test and p-value was set to .05.

To screen for functional networks among the differentially methylated and expressed genes, we employed the Functional Networks of Tissues in Mouse (FNTM) prediction tool for mouse tissue-specific protein interactions, which integrates genomic data and prior knowledge of gene function (Goya et al., 2015). We used the hippocampus tissue database and only kept edges with a relationship confidence greater than 0.75.

An interactome network model integrating transcriptomics and metabolomics was constructed with Cytoscape MetScape (v3.1.3).

Statistical analyses of results of metabolites content were performed with GraphPad Prism (version 8.0). Statistical significance was determined by two-tailed unpaired Student’s t-test. p-value < .05 was considered statistically significant.

We used transmission electron microscopy (TEM) to examine nuclei in the CA1 region of hippocampal neural cells. In a part of neurons, the chromatin aggregated into several clumps (Figure 1D, arrows); the shape of the nucleus was irregular; and the invaginations of nucleolemma became evident indicating the shrinkage of the nuclei (Figure 1D, asterisk), suggesting that the nuclei of neural cells were impaired in the GDM group compared to the Ctrl group (Figure 1).

We performed RRBS to search for DMRs in GDM group vs. Ctrl group. Initial principal component analysis (PCA) of DMLs revealed a clear distinction between GDM and control samples (Supplementary Figure S2; Supplementary Table S1). The distribution of DMRs in the different functional components of the genome was illustrated in Figure 2A in the upstream 2 k (5.64%), 5′-untranslated region (5′-UTR, 1.06%), coding sequence (CDS, 11.06%), introns (37.31%), 3′-UTR (3.11%), downstream 2 k (4.4%), and other elements (37.42%) of genes. In GDM group, a total of 429 DMR-related genes had lower methylation levels, and 321 DMR-related genes had higher methylation levels (Supplementary Table S2). The top 10 hypomethylated genes (Rbm19, Tox2, Erg, Irs1, Sox5, Sorcs3, Sgol2b, Specc1, Gm5441, Ptprn2) and hypermethylated genes (Unc45bos, Unc45b, Ash1L, Olah, Serpinf2, Tmco1, Lyst, Wdr66, Gm20388, Fendrr) were shown in Figure 2B. Among them, Sox5 controls neural differentiation and genesis (Lai et al., 2008). Sorcs3 is a stronger regulator of glutamate receptor functions compared to GABAergic mechanisms in the hippocampus (Christiansen et al., 2017). The top 10 terms of the GO biological process analysis were shown in Figure 2C, among which “dendrite development,” “regulation of dendrite development,” and “synapse organization” were thought to be important for cognition, learning and memory. In addition, the top 20 pathways of KEGG enrichment were shown in Figure 2D, and some of these enriched terms were involved in “axon guidance,” “glutamatergic synapse,” and “long-term depression,” which were closely related to cognitive function.

The gene expression profiles of fetal mice hippocampi at embryonic day 18.5 from the Ctrl and GDM groups were analyzed using RNA-seq. Initial principal component analysis (PCA) revealed a clear distinction between GDM and control samples (Figure 3A). By using a stringent threshold and significant criteria of FDR <0.05 and fold change >1.2, a total of 953 DEGs, including 378 up-regulated and 575 down-regulated genes, were identified in the GDM group compared to the Ctrl group (Figure 3B, Supplementary Table S3). Among the top 20 DEGs, 11 genes (Ndufaf4, Nrgn, Hap1, Sv2b, Slc12a5, Sncb, Caln1, Kcnma1, Abi1, Trim9, Shroom2) were down-regulated, and 9 genes (Hmgb2, Kctd20, H1f0, Rnf13, Atg7, Knstrn, Csnk1g1, Crk, Gm26786) were up-regulated (Figure 3C). Nrgn is expressed exclusively in the brain, particularly in dendritic spines and has been implicated in hippocampal plasticity. Hap1 depletion impaired postnatal neurogenesis in the dentate gyrus (DG) of the hippocampus (Xiang et al., 2015). Sv2b has been shown to be involved in the regulation of synaptic vesicle exocytosis and plays a role in neuronal excitability (Stout et al., 2019). Caln1 participates in the regulation of trans-synaptic signaling (Reshetnikov et al., 2020).

Functional enrichment analysis identified KEGG pathways including “GABAergic synapse,” “neuroactive ligand-receptor interaction,” “glutamatergic synapse,” “cAMP signaling pathway,” “synaptic vesicle cycle” and “long-term potentiation” related to brain cognition. The top 10 KEGG pathways were illustrated in Figure 3D. In addition, GO analysis revealed that the terms “synapse organization,” “neurotransmitter transport,” “signal release from synapse,” “neuron to neuron synapse” and “postsynaptic neurotransmitter receptor activity” were enriched (Figure 3E).

We further identified DNA methylation changes linked to gene expression differences in fetal mice hippocampi in the GDM group compared with the Ctrl group. A Venn diagram was constructed to illustrate 43 shared genes of the methylome and transcriptome (Figure 4A). A scatter diagram of the 4-quadrant plots showed that 16 genes (Cdipt, Mfhas1, Slc7a14, Gm38393, Pcdha6, Paqr9, Rexo2, Tbc1d30, Akap6, Kcnf1, Sv2c, Arhgap26, Ank1, Arhgap23, Shank2, Sez6l) were identified as hypermethylated and down-regulated, and 9 genes (Thtpa, Sema5b, Lima1, Gse1, Tmem98, Mpped2, Peg12, Mcm10, Mid1) were hypomethylated and up-regulated in the GDM group. The other 18 genes did not show a canonical association of gene expression and methylation status (Figure 4B).

GO analysis showed that the above overlapping genes were mainly enriched in “long-term synaptic depression,” “negative regulation of synaptic transmission” and “learning and memory” (Figure 4C). Functional enrichment analysis identified 17 KEGG pathways that were significantly different between the GDM and Ctrl groups. The top 10 KEGG pathways included “GABAergic synapse,” “glutamatergic synapse,” “serotonergic synapse” and “dopaminergic synapse,” mainly relating to information transfer and cognition (Supplementary Figure S3).

Among all 43 genes that were differentially methylated and expressed, we obtained the seven most functionally connected predicted interactor genes visualized using the FNTM tool (Pcdha6, Ank1, Sez6l, Kcnb1, Kcnj6, Prkcb, Sorcs3) (Figure 4D), all of which were down-regulated in transcriptome profiling. These genes are involved in synapse signal transmission, such as Kcnj6 and Prkcb, and function in cognition, such as Ank1, Kcnb1 and Sorcs3.

The UPLC-MS-based metabolic profiling method was performed to explore the potential biochemicals involved in the epigenetic regulation of the GDM group. The GDM group could be distinguished from the Ctrl group in the OPLS-DA score plot in both positive and negative ionization modes, which suggested that intrauterine hyperglycemia could induce distinct changes in fetal brain metabolites (Figures 5A,B).

In total, 123 significantly differential metabolites were identified in fetal brain tissues using the threshold of p < .05 and VIP ≥ 1, of which 8 up-regulated and 115 down-regulated in GDM group (Supplementary Table S4). Metabolites determined under the condition of |log₂FoldChange| > 1.5 are shown in Figure 5C. Differential metabolites were in various classes, including “amino acid” (39%), “nucleotide” (11%), “carbohydrate” (15%), “lipid” (14%), “cofactors and vitamins” (11%), “xenobiotics” (7%), “peptide” (1%) and “energy” (1%) (Figure 5D). The top 30 metabolites included nicotinamide riboside, AICAR, gamma-aminobutyric acid, 2-hydroxybutyric acid, methyl jasmonate, choline, 3-hydroxyanthranilic acid and adenosine, most of which were related to cognition (Figure 5E). Further target metabolic pathways were screened out according to enrichment analysis (p value < .05) and pathway topological analysis (impact value >0.10). A total of 51 metabolic pathways were altered in response to intrauterine hyperglycemia, including “cysteine and methionine metabolism,” “citrate cycle (TCA cycle),” “pantothenate and CoA biosynthesis,” and “one carbon pool by folate,” (Figure 5F, Supplementary Table S5). In the GDM fetal brains, metabolites related to epigenetic modification were significantly down-regulated, including betaine, fumaric acid, L-methionine, succinic acid, 5-methyltetrahydrofolic acid, and S-adenosylmethionine (SAM) (Figure 6).

An interactome network model integrating the transcriptome and metabolome was generated, which connected pathways via gene-metabolite interactions. Analysis of this interactome network demonstrated that several critical metabolic processes were altered in the brains of GDM fetal mice, including “TCA cycle,” “tryptophan metabolism,” “tyrosine metabolism” (Figure 7). Other pathways included “purine metabolism,” and “urea cycle” and “metabolism of arginine, proline, glutamate, aspartate and asparagine” (Supplementary Figures S4, S5). The DEGs and differential metabolites associated with these pathways are presented in Table 1. Predicted genes were mostly belong to gene families such as cytochromeP450 family, which was likely to play roles in “tryptophan metabolism” and “tyrosine metabolism.” Predicted compounds such as citrate, succinyl-CoA in TCA cycle may interact with input compound and genes and co-mediate cognition disorder in GDM offspring. The results suggested that amino acid metabolism and TCA cycle may be the key processes associated with cognition of GDM offspring.