Maternal Undernutrition Induces Cell Signalling and Metabolic Dysfunction in Undifferentiated Mouse Embryonic Stem Cells

Peri-conceptional environment can induce permanent changes in embryo phenotype which alter development and associate with later disease susceptibility. Thus, mouse maternal low protein diet (LPD) fed exclusively during preimplantation is sufficient to lead to cardiovascular, metabolic and neurological dysfunction in adult offspring. Embryonic stem cell (ESC) lines were generated from LPD and control NPD C57BL/6 blastocysts and characterised by transcriptomics, metabolomics, bioinformatics and molecular/cellular studies to assess early potential mechanisms in dietary environmental programming. Previously, we showed these lines retain cellular and epigenetic characteristics of LPD and NPD embryos after several passages. Here, three main changes were identified in LPD ESC lines. First, their derivation capacity was reduced but pluripotency marker expression was similar to controls. Second, LPD lines had impaired Mitogen-activated protein kinase (MAPK) pathway with altered gene expression of several regulators (e.g., Maff, Rassf1, JunD), reduced ERK1/2 signalling capacity and poorer cell survival characteristics which may contribute to reduced derivation. Third, LPD lines had impaired glucose metabolism comprising reduced upstream enzyme expression (e.g., Gpi, Mpi) and accumulation of metabolites (e.g., glucose-6-P, fructose-6-P) above the phosphofructokinase (PFK) gateway with PFK enzyme activity reduced. ESC lines may therefore permit investigation of peri-conceptional programming mechanisms with reduced need for animal experimentation. Graphical Abstract Supplementary Information The online version contains supplementary material available at 10.1007/s12015-022-10490-1.


Introduction
The origin of adult chronic non-communicable diseases across the cardiometabolic and neurological spectrum can arise from the prenatal environment according to the Developmental Origins of Health and Disease (DOHaD) hypothesis [1]. In particular, the period around conception is vulnerable to diverse stress conditions including sub-optimal parental nutrition and assisted reproductive treatments (ART) [2][3][4]. Adverse peri-conceptional environment can alter multiple molecular, cellular and physiological processes in the developmental programme that persist through gestation affecting growth, metabolism and later disease risk, across mammalian species, including humans [2,5,6]. Whilst challenging environmental conditions can evoke developmental plasticity to aid survival, altered programming can become maladaptive where conditions later change leading to pathophysiology [1,2,4].
We have used a mouse maternal protein restriction model whereby low protein diet (LPD), administered exclusively during preimplantation development with normal nutrition thereafter and postnatally, is sufficient to induce adult offspring cardiovascular, metabolic, structural and neurological dysfunction [7][8][9][10]. This transient restriction alters uterine metabolite and amino acid levels sensed directly Pooja Khurana, Andrew Cox and Barira Islam contributed equally by the preimplantation embryo through mTOR signalling [11] leading to an array of cellular responses affecting both embryonic [8][9][10]12] and extra-embryonic [11,[13][14][15][16] lineages. Given the limited size and inaccessibility of the preimplantation embryo, we have used embryonic stem cell (ESC) lines derived from blastocysts to probe effects of maternal LPD in early embryonic programming. This novel approach has shown common molecular and epigenetic changes induced by maternal diet in both early passage ESC lines and embryos in vivo [13,14]. Moreover, ESC lines have been employed for mechanistic analysis of other peri-conceptional conditions associated with DOHaD [17,18].
Here, we provide a detailed characterisation of undifferentiated LPD and control (normal protein diet, NPD) ESC lines for further evidence of early steps in adverse programming and phenotypic plasticity. Through RNA sequencing, metabolomic, integrated bioinformatics, qRT-PCR, cellular and physiological analyses, we identify that maternal undernutrition leads to dysfunction in MAPK cell signalling and glycolytic metabolic pathways. We propose that these cells provide a valuable resource for DOHaD mechanistic studies, and can be used to generate in vitro models for studying developmental programming events in early embryogenesis or in specific organs through directed differentiation or organoid cultures.

Animal Treatment, ESC Line Derivation
Animal procedures were conducted under UK Home Office project license and local ethics approval (University of Southampton). C57BL/6 females (7-9 week) were mated naturally overnight with C57BL/6 males and plug positive females were housed individually the following morning and assigned randomly to either control normal protein diet (18% casein, NPD) or isocaloric low protein diet (9% casein, LPD) supplied by Special Dietary Services Ltd, UK, and fed ad libitum until embryonic day 3.5 (E3.5); diet composition has been described elsewhere [7]. Embryos were collected at the blastocyst stage after cervical dislocation and uterine flushing with H6 medium with 4 mg/ml BSA (H6 + BSA) [19]. Blastocysts were placed onto a freshly-prepared feeder layer of mitotically inactivated mouse embryonic fibroblasts (MEFs) prepared using standard procedures [20] 10,378) and 1000 U/ ml of leukaemia inhibitory factor (LIF)), incubated at 37 °C in 5% CO 2 . After substantial outgrowth of the inner cell mass (ICM; 2.5 days culture), ICM-derived clumps from individual embryos were dislodged using a drawn capillary tube and mechanically dissociated into single-cell and smallcell aggregates in 0.05% trypsin-EDTA. Disaggregated cells were seeded onto fresh feeder layers on 96 well plates and cultured in ES medium. Viable ESC colonies were primarily passaged at day 3.5, expanded by trypsin passaging at 3-day intervals to reach confluence on a 100 mm dish. ESC lines were frozen after ~ 24 days in culture, passage number (PN) 5-7. Sex determination of ESC lines was carried out after feeder cell depletion, by amplification of Y chromosome Sry and Zfy and control X chromosome DXNds3 [21,22]. Male ESC lines were predominant (LPD 70.6%, NPD 80.0%). ESC lines were karyotyped after at least 10 passages by fixation and Giemsa staining with most showing euploid karyotype (40 chromosomes in > 50% metaphase spreads; (LPD 71.4%; NPD 80.0%). A total of 18 ESC lines were generated from LPD blastocysts from 11 mothers and 38 lines from NPD blastocysts from 11 mothers. Male normal karyotype cell lines, all from different mothers, were used in all analyses.

RNA Sequencing and Data Analysis
mRNA sequencing of five LPD and five NPD ESC lines (PN 9-12) was performed by an Illumina platform. RNA was extracted from ESCs cultured in fibroblast-free media using an RNeasy mini kit (Qiagen, UK). RNA was quality checked using nanodrop, agarose gel electrophoresis and Agilent 2100 Bioanalyzer followed by enrichment using oligo(dT) beads. The resulting mRNA was fragmented randomly and used as the template for synthesis of first strand of cDNA using random hexamers primer. The second strand synthesis was carried out by removing the RNA template and synthesizing a replacement strand to generate double-stranded cDNA. A series of terminal repair was carried out by ligation and adapter-ligation. The double-stranded cDNA library was completed through size selection and PCR enrichment. The quality control and concentration of the library was assessed using Qubit 2.0, Agilent 2100 bioanalyzer and qRT-PCR. The qualified libraries were fed into Illumina sequencers after pooling according to their effective concentration and expected data volume. The original data obtained from the high throughput sequencing platforms were transformed into sequenced reads by base calling. Reads were then filtered for adaptor contamination using cutadapt [23] and filtered such that at least 90% of bases of each read had a quality score > 20. Raw RNA sequencing data has been submitted to GEO repository (Accession number GSE216108).
The program Salmon was used for quantifying abundances of transcripts from RNA-sequencing data, or 1 3 more generally of target sequences using high-throughput sequencing reads [24]. The index for Mus_musculus. GRCm38.88_transcripts.fa was built in Salmon using the default parameters. The program Wasabi was used to prepare Salmon output for downstream processing by Sleuth. Bootstrapping of the quantification was performed iteratively for 100 rounds. Resulting counts were then normalized through the transcripts per kilobase million (TPM) method. Differential expression analysis was carried out using the software package Sleuth (v0.30.0) using TPM normalized counts for each cell type. Gene Set Enrichment Analysis (GSEA) was carried out on the Broad Institute website (http:// www. broad. mit. edu/ gsea) using Hallmark Gene Sets of Molecular Signatures Database (MSigDB) collection because it represents 50 well-defined, large-scale biological processes that are applicable to the widest range of potential cellular responses [25]. The mouse genes were mapped to human genes before carrying out the analysis. A ranking metric was assigned to all the annotated genes using their P-values and fold change [26]. GSEA calculates unique enrichment score for each gene set that reflects the degree to which a gene set is overrepresented at either end of this continuum.

Metabolomics and Data Analysis
Seven ESC lines were cultured (in duplicate) from each dietary group (PN [10][11][12], and expanded to 70-80% confluency in serum-free ES medium (KnockOut™ DMEM, LIF, 15% KO/SR), counted using trypan blue, flash frozen in liquid nitrogen as 50-100 μl packed cell pellets and stored at -80 °C. Samples were further processed using an automated MicroLab STAR® system (Hamilton) and analysed commercially (Metabolon Inc., NC USA) using ultra-performance liquid chromatography (UPLC, Waters ACQUITY) and a high-resolution tandem mass spectrometry (MS/MS, Thermo Scientific Q-Exactive) platform. Curated relative ion intensity data were normalised by both protein concentration (Bradford assay) and cell counts (trypan blue assay) analysis, and missing values were imputed to ensure improved data quality and to preserve biological variability of the samples. Metabolite values were normalised (i.e. rescaled) according to the median and log-transformed and back-transformed (as simple geometric means) and significance was evaluated by the Welch's 2-sample t-test. The dataset was also adjusted for false discovery rate (FDR, q-value).

Integrated Metabolomics and Transcriptomics Analysis
Paintomics version 4 was used to carry out integrated transcriptomics and metabolomics analysis using our two datasets [27]. The program integrates metabolite and transcript abundances and significances and maps them onto organism-specific Kyoto Encyclopedia of Genes and Genomes (KEGG) maps. All the annotated transcripts and metabolites with their fold-change values were imported and a list of the significant transcripts (P < 0.05) and metabolites (P < 0.05) were used as feature for respective analysis. The program computed P-values for individual transcriptomics and metabolomics analysis as well as combined P-values. The joint pathway enrichment P-value considering both the omics data was computed by applying the Fisher combined probability test. This is done by assigning a weight to transcriptomics and metabolomics that is the percentage of features mapped to pathways in that omic.

Quantitative Real-time PCR
RNA was extracted as for above RNA sequencing for quantitative real-time PCR (qRT-PCR). cDNA was prepared using the Verso cDNA synthesis Kit (Thermo Fisher Scientific, Waltham, USA). A mix comprising 2 μl 5 × cDNA synthesis buffer, 1 μl deoxyribonucleotides triphosphate (dNTPs), 0.5 μl random hexamers, 0.5 μl oligo (dT) primers and 0.5 μl Verso enzyme mix was added to 150 ng RNA. The mix was then incubated at 42 °C for 60 min, followed by enzyme inactivation at 95 °C for 2 min as for manual instructions. The resulting cDNA was diluted 1:5 in Tris-EDTA (TE) buffer (Thermo Fisher Scientific) and stored at -20 °C until use. Primers used for qRT-PCR were designed using the online tools found in the Primer3 website (http:// bioto ols. umass med. edu/ bioap ps/ prime r3_ www. cgi) and were designed to produce amplification products < 300 bp, and to have a melting temperature of 60 ± 1 °C (Supplementary Table 1). Primers were manufactured by Eurofins Genomics (Edersberg, Germany). For each sample, 1.5 ng of cDNA was used for qRT-PCR using a CFX96 RT instrument (Bio-Rad, Hercules, USA). Analysis of samples was performed in duplicate with each 10 μl reaction composed by 3 μl iTaq Universal SYBR Green Supermix (Bio-Rad), 300 nM of each forward and reverse primer and 5 μl diluted cDNA (0.3 ng/μl). The qRT-PCR reaction involved polymerase activation and cDNA denaturation at 95 °C for 5 min, followed by 40 thermal cycles consisting of denaturation at 95 °C for 5 s and primers annealing and extension at 60 °C for 30 s. Fluorescence was measured at the end of every cycle. Upon completion of thermal cycling, melting curve analysis was performed to confirm reaction specificity. Bio-Rad CFX Maestro software (Bio-Rad) was used to perform baseline subtraction and determine the threshold cycle (Ct). Normalization of expression of gene of interest to the selected reference genes (Sdha, Tbp and Tuba-1; GeNorm) was performed applying the formula for the calculation of the ΔCt. Data normalization and ΔΔCt calculations were performed using the CFX Maestro software (Bio-Rad).

3
The expression fold change difference was calculated as 2 ^(− ∆∆Ct) For pluripotency gene expression analysis, cDNA was synthesised by random priming from 250 ng RNA using the Improm-IITM Reverse Transcriptase System (Promega, UK) kit, amplification was conducted using a Chromo4 Real-time Detector (BioRad, UK) with Opticon Monitor v3.1 software, and amplified target genes were quantified following the SYBR-green precision mastermix (PrimerDesign, UK) protocol and using the same reference genes as above.

ESC Survival and Apoptosis
ESCs (4 × 10 4 cells) on MEFs in ES medium were harvested and trypsinized after 24, 48 and 72 h culture, stained with Annexin V-FITC and PI (100 µg/ml) in Annexin V binding buffer (10 mM HEPES pH 7.4, 140 mM NaCl, 2.5 mM CaCl 2 ) for 30 min at 4 °C in the dark. For flow cytometric analysis, 20,000 cell events were recorded using a FACS-Calibur flow cytometer (BD Biosciences). Data was assessed using Cell Quest Pro software (BD Biosciences). Forward scatter (indicating cell size) and side scatter (indicating cell granularity) gates were set so as to exclude MEF cells and debris. The percentage of apoptotic cells was determined by generating a dual-colour dot plot (Annexin V-FITC (FL-1) vs. PI (FL-2)) and then setting a quadrant marker based on unstained and single-labelled control samples.

Western Blotting
ESCs were grown feeder-free (3 × 10 5 cells) and cell lysates prepared at 24, 48 and 72 h by incubation for 20 min with Onyx lysis buffer (20 mM Tris-HCl, pH 7.4, 135 mM NaCl, 1.5 mM MgCl 2 , 1 mM ethyleneglycoltetraacetic acid (EGTA), 1% TX-100, 10% glycerol) supplemented with both a protease inhibitor cocktail (Sigma. 1:100) and phosphatase inhibitors (2 mM Na 3 VO 4 , and 50 mM NaF). Reduced samples were fractionated by SDS-PAGE on pre-cast 4-12% NuPAGE® Bis-Tris polyacrylamide gels (Invitrogen). After separation, proteins were transferred to a methanol-activated Polyvinylidene Fluoride membrane (Immobilon-FL; Millipore) by overnight electro-elution. Membranes were blocked with 5% (w/v) non-fat dry milk in Tris-Buffered Saline (TBS) + 0.05% Tween20 (TBS-T) for 1 h at room temperature with constant agitation. Membranes were washed with TBS-T and probed with primary antibodies (Supplementary Table 2) overnight at 4 °C, washing and incubation with secondary antibody conjugated to IRDye™ 800CW (Rockland Inc. Molecular Probes) near-infrared (IR) fluorochrome for 1 h at room temperature in the dark. Following further washing, blots were scanned using the LI-COR Odyssey IR imaging system (LI-COR Biosciences). Integrated density values (IDVs) were obtained using the Odyssey software based on pixel intensity, and normalised to α-tubulin loading controls. Where membrane stripping was required, blots were incubated in stripping buffer (200 mM glycine, 3.5 mM SDS, 1% (v/v) Tween20, pH 2.2) for 10 min at room temperature, and then subsequently washed in PBS for 20 min and TBS-T for 10 min.

Statistical Analyses
Statistical analysis of RNA sequencing, metabolomics and integrated Paintomics comprised tools available within the software programmes used and are referred to in the respective sections for generation of P values. The qRT-PCR, enzyme assays, and cellular studies were performed a minimum of three times and Minitab 16.1.0 or GraphPad Prism version 8.0 were used for statistical testing. Normality was tested using the Anderson-Darling or Shapiro-Wilk tests and considered normal when P > 0.05. Variance homogeneity was analysed using the F-test. Differences between dietary treatment groups were performed using the independent Student's t-test or one way ANOVA (for normally distributed data), or Mann-Whitney rank sum test (for non-normal data sets). For the metabolic enzyme assay, differences were calculated using Welch's t-test. All data are expressed as means ± SEM. P < 0.05 was considered statistically significant.
To investigate biological distinctions between established LPD and NPD cell lines at low passage number, we used a multi-omics approach comprising RNA sequencing (summarised Supplementary Fig. 2; Table 1), metabolomics (summarised Supplementary Fig. 3; Table 2) and integrated analysis using PaintOmics (summarised Supplementary Fig. 4; Table 3), together with supportive molecular and cellular studies. Two main distinctions were revealed as presented below.

LPD Impairs MAPK Network Signalling in ESCs Associated with Enhanced Cellular Stress and Apoptosis
RNA sequencing indicated an overall similar phenotype between LPD and NPD ESC lines (Supplemetary Fig. 2A) but with differences in transcript expression identified in volcano plot (Supplementary Fig. 2B) and with 24 transcripts passing false discovery rate (FDR) at < 0.1 (listed Table 1) and 8 passing FDR at < 0.05 further visualised by heat map analysis ( Supplementary Fig. 2C). Gene set enrichment analysis (GSEA) indicated enrichment of several pathways in LPD related to inflammation, energy metabolism and survival, collectively suggesting broad stress effects ( Supplementary Fig. 2D). Three out of eight of the FDR < 0.05 transcripts, Maff, Atf4, JunD, were increased in LPD (all P < 0.0001) and are associated with the Mitogenactivated protein kinase (MAPK) pathway, a network of serial kinases with downstream effectors that target genes to coordinate cell proliferation, differentiation, apoptosis and stress response [28,29]. The MAPK network was also the top significant class identified in PaintOmics, an integrated analysis of our transcriptomics and metabolomics datasets (Supplementary Fig. 4; Table 3).
Since MAPK signalling affects cell proliferation, this was measured for the ESC lines firstly during initial cell expansion at early passage number (PN 4-6) before cells were frozen, and secondly after cell expansion and freezing (PN 9-12). Viable cell counts comparing the change in cell numbers over 24, 48, 72 and 96 h after passage were determined and  (Fig. 2C). The rate of BrdU incorporation at 24 h to identify S-phase cells was also similar between LPD and NPD lines ( Supplementary Fig. 5B). Lastly, ESCs synchronised using nocodazole and subsequently released displayed similar cell cycle characteristics across diet treatments ( Supplementary Fig. 5C).
MAPK signalling is also critical in regulating cell survival and apoptosis [28,29] and so this was evaluated next. Trypan blue staining of cells at 24-96 h after plating showed greater proportions of dead (trypan +) cells in LPD ESC clones (Fig. 2D) suggesting reduced survival. This was most evident at 48 h (LPD 10.01 ± 0.331%, NPD 6.93 ± 0.552%; P = 0.003) and 96 h (LPD 10.71 ± 1.256%, NPD 7.64 ± 0.572%; P < 0.05). This was substantiated by flow cytometric assessment of cell lines following FITC-Annexin V staining which binds to phosphatidylserine, externalised during early apoptosis, and PI for loss of cell integrity (late apoptosis). A representative density plot of AnnexinV-PI staining by flow cytometry is shown in Supplementary   Fig. 5D. LPD ESC lines exhibited fewer live cells and increased early and late apoptosis and total dead cells over 72 h culture compared with NPD lines (P < 0.05 -0.01) (Fig. 2E). Quantitative immunoblotting of ESC lines over 72 h culture was used to assess cell signalling activity that may explain the dietary-induced change in cell apoptosis. ESCs maintain their pluripotent state through LIF signalling which activates the STAT3 cascade to retain undifferentiated self-renewal and regulate cell survival, through activation of PKB/AKT signalling [31,32]. We analysed activation of AKT (phosphorylation at ser493) and STAT3 (phosphorylation at tyr705) but neither indicated a change in activity or total levels between diet groups ( Supplementary Fig. 6A, B). MAPK signalling is known to be induced in ESCs through ERK1/2 and the stress-related p38 pathways [31,[33][34][35][36][37]. We tested the activation of ERK1/2 (phosphorylation at thr202/ tyr204) and p38 (phosphorylation at thr180/tyr182) and found that whilst both had similar total protein levels across diet groups, ERK1/2 activation, but not that of p38 (Supplementary Fig. 6C, D), was reduced by ~ 50% in LPD ESCs (P < 0.0016-< 0.00001 over 24-72 h culture) (Fig. 2F, G). Lastly, analyses of several components in the apoptotic pathway (PUMA, BimEL, Caspase 3, Cleaved caspase 3; [38]) showed no difference across diet treatment (Supplementary   Fig. 6C, D). These data show that altered MAPK pathway activity through reduced ERK1/2 signalling in LPD ESC lines coincides with their increase in apoptosis and cell death rate and reduced derivation efficiency.
The changes identified in glycolytic metabolites and associated gene expression following LPD, with accumulation of upstream metabolites coupled with reduced enzyme expression but with depleted metabolites downstream, is suggestive of a bottleneck in processing at the PFK step. Thus, enzyme activity assays were performed on both hexokinase (HK) glucose to G6-P conversion, and on phosphofructokinase (PFK) F6-P to F-1,6-bP conversion. Although no difference was observed in HK activity (P = 0.2), PFK activity in LPD lines was reduced (P = 0.03, Fig. 3F).

Discussion
Our characterisation of ESC lines derived from mothers fed LPD during the preimplantation period was conducted to identify early mechanisms of adverse developmental programming that can lead to increased risk of non-communicable diseases in later life [2,4]. Analyses of LPD ESC lines compared with control NPD lines using RNA sequencing, metabolomics, integrated bioinformatics, qRT-PCR and cellular studies revealed two major dysfunctional phenotypic states comprising altered gene expression and signalling activity of members of the MAPK signalling network (summarised Fig. 4A) and changes in the expression and activity of the glycolytic energy metabolism pathway (summarised Fig. 4B). We also found reduced efficiency in the derivation of LPD ESC lines alongside more subtle and transient changes in expression of pluripotency markers in undifferentiated cells which may be linked to the MAPK phenotype. Also, omics approaches identified additional changes in signalling pathways (see, for example, Supplementary Fig. 2D  and 4) that may further modulate phenotype which, while not corroborated here with additional studies, are discussed below to contextualise our outcomes.
MAPK signalling is a highly conserved network of serial kinases belonging to JNK, p38 and ERK classes responsive to environmental stimuli which phosphorylate downstream effectors to activate or repress expression of target genes that regulate fundamental cellular activities such as proliferation, differentiation, apoptosis and stress response [28,29]. Whilst our RNA sequencing and qRT-PCR analyses showed only minimal changes in expression of specific MAPK kinases and accessory proteins (eg, Map3k14; Map-kapk2), pronounced and consistent changes were identified in downstream effectors regulating transcriptional activity. Thus, LPD lines displayed a significant upregulation of several MAPK transcripts including JunD, Atf4, Fosl1, Fosl2 and Maff which collectively encode subcomponents of the activating protein-1 (AP-1) homo-or hetero-dimer transcription factor complex comprising the integrated stress response [39,40]. The AP-1 dimerised elements bind to cAMP response element (CRE) sites to effect downstream expression associated with cell survival or apoptosis in response to various stresses including amino acid starvation, ER stress and oxidative stress [41][42][43]. We note that amino acid deficiency in the uterine fluid has been identified as an initiator of LPD programming in the mouse preimplantation embryo leading to later life health risk [11,44] and so may be responsible causally for AP-1 transcript changes in LPD ESC lines. The combination and abundance Fig. 2 LPD causes impaired MAPK signalling with reduced survival and increased apoptosis in ESC lines. A RNAseq above and qRT-PCR analysis below of most differentially expressed MAPK components in undifferentiated LPD and NPD ESC lines (PN 8-12; n = 5) cultured in standard mESC medium (KnockOut DMEM + LIF + KO/ SR). qPCR gene expression is normalised with 3 house keeping genes Sdha, Tbp and Tuba-1. Here, * = P < 0.05, ** = P < 0.01, *** = P < 0.001 and **** = P < 0.0001. B Proliferation rate of LPD and NPD ESC lines (P9-12) over 24 h intervals using trypan blue exclusion and haemocytometer counting. Data from 8 LPD or 9 NPD lines with mean ± SEM, each analysed in triplicate. C Cell cycle distribution of exponentially growing ESC lines (P10-12) 24 h after seeding, assayed by DNA quantitation using flow cytometry following PI staining and FlowJo software. Data from 10 ESC lines per treatment with mean ± SEM, each analysed in triplicate. D Cell viability assay. Mean percentage (± SEM) of dead cells from Trypan blue exclusion assay over 96 h culture (PN 11-13). LPD n = 8 ESC lines from 8 mothers; NPD n = 9 ESC lines from 7 mothers; ***P < 0.001, *P < 0.05. E Cell death assay. Percentage (± SEM) ESC lines (PN 11-12) after culture, annexin V and PI staining and flow cytometry discriminated as live, apoptotic and necrotic. LPD n = 8, 6, 7 lines at 24, 48 and 72 h; NPD n = 10, 9, 6 lines at 24, 48 and 72 h; **P < 0.01, *P < 0.05, #P < 0.1. F, G ERK1/2 immunoblots (F) and quantitation (p/total) (G) of ESC lines cultured up to 72 h. 13 + ESC lines per diet from minimum of 3 blots (independent experiments); α-tubulin used as loading control; ****P < 0.0001, **P = 0.0016. For antibody sources, see Supplementary  of AP-1 dimerised partners can influence whether positive or negative regulation of downstream targets occurs leading to apoptotic or proliferative cell survival outcomes. For example, FOSL1 is expressed in tumour tissue where it promotes growth and differentiation rather than apoptosis [43]. In addition to increased expression of AP-1 mRNAs, LPD lines also expressed increased Dusp8 and Dusp14 mRNAs which encode members of the dual specificity phosphatase (DUSP) protein family that dephosphorylate MAPK kinases including ERKs as part of the regulation of stress response towards survival or apoptosis [45,46]. LPD lines were also characterised by increased Rassf1 expression which belongs to the RASSF protein family that links RAS to the MAPK pathway. As in above examples, RASSF has been shown to contribute to diverse cellular functions across the range from apoptotic to proliferative activities [47]. Further, LPD lines were downregulated in Ddx17 which encodes multitasking RNA helicases that similarly modulate cellular activity from apoptosis to proliferation [48]. Lastly, LPD lines exhibited upregulation of Nppb encoding B-type natriuretic peptide (BNP) which contributes to MAPK-mediated stress response in cardiomyocytes [30]. BNP plays a role in cardiac development and elevated plasma levels associate with heart failure and myocardial infarction [49]. Maternal undernutrition can also cause elevated BNP and increased cardiac dysfunction during development [50] but effects mediated through proliferation or apoptosis are unknown.
To evaluate the functional implications of MAPK pathway modulation in LPD lines, we first examined cell proliferation and cell cycle control. In detailed studies, we were not able to identify any changes in the rate of proliferation over 96 h using early or later passage cells, nor were differences found in cell cycling stages using flow cytometry, BrdU incorporation, and synchronised cell populations. In contrast, analysis of cell survival and apoptosis showed LPD ESCs had increased proportion of dead cells within cultures over time, and, following flow cytometry, had fewer live cells and increased early and late apoptotic and dead cells. Thus, MAPK signalling changes in LPD lines more likely relate to poorer survival than regulation of cell proliferation. Moreover, the poorer rate of establishment of LPD lines is consistent with this interpretation.
To assess directly MAPK signalling, we quantified by immunoblotting the posttranslational state of MAPK classes, p38 and ERK, known to promote apoptosis in mESCs [35]. Whilst no change was found in p38 signalling, ERK1/2 was significantly dephosphorylated in LPD lines, representing relative ERK1/2 inactivation. This altered state is consistent with the increase in apoptosis and cell death found in LPD lines by flow cytometry but is likely to be measured since the posttranslational state of selected apoptotic markers in LPD lines was not altered. This distinction between the flow cytometry and apoptosis immunoblotting data may reflect the relatively small changes in apoptotic cell percentages identified (Fig. 2E), the relative extent of apoptotic protein expression and posttranslational processing, or the timing of effector activation of proteins (eg, caspase-3) not synchronising with later morphological status of cells screened by flow cytometry (eg, Annexin V binding; PI internalisation). The altered MAPK and ERK1/2 phenotype is also consistent with our previous findings that LPD neural stem cells during fetal brain development are reduced in number partly through increased apoptosis [9].
Here, the upstream steps in the 10 enzyme pathway of cytoplasmic glycolytic metabolism from glucose to pyruvate were affected (summarised Fig. 4B). Whilst glucose levels were unchanged, the glycolytic intermediates G6-P and F6-P were both significantly increased in LPD ESC lines; similarly, the intermediates from mannose metabolism, M6-P which interconverts with F6-P, and its end-product G.amine 6-P, were also increased at trend level in LPD cells. RNA sequencing and qRT-PCR both showed that the expression of the key enzymes leading to F6-P production, namely Gpi1 from G6-P and Mpi from M6-P, were both significantly reduced in LPD lines. In contrast to the elevated levels of upstream sugars, glycolytic metabolites downstream from F6-P showed some evidence of reduction in LPD. Thus, F1,6 BP, DHAP, 3-PGlyc, PEP and pyruvate were all reduced in LPD but not individually reaching significance. This glycolytic pattern (Fig. 4B) suggests a 'bottleneck' may occur in the pathway at F6-P, limiting downstream processing. F6-P metabolism to F1,6 BP is a rate-limiting step for glycolysis mediated through phosphofructokinase (PFK) whose activity is primarily regulated through dimer and tetramer formation and stabilisation of the protein [55]. Moreover, PFK expression and activity has been shown to be developmentally regulated from the embryo onwards [56,57]. Genes encoding PFK and downstream enzymes in glucose metabolism were mainly unaltered in RNA sequencing or qRT-PCR although minimal Pfk isoforms were divergently changed and pyruvate kinase (Pkm; converting PEP to pyruvate) were downregulated in LPD, possibly contributing to impaired glycolysis. Thus, the putative F6-P bottleneck may reflect either transcriptional (Pfkp), translational, and/or posttranslational modification in PFK enzyme activity. PFK activity was indeed shown to be significantly reduced in LPD cell lines whilst HK activity which acts on glucose conversion to G6-P was unchanged.
These data therefore indicate LPD ESC lines show deficiency in glycolytic metabolism through reduced PFK activity possibly combined with reduced expression of Gpi1, Mpi and Pkm. There is also suggestion that the accumulation of upstream intermediates coupled with reduced PFK results in increased redirection of F6-P to G.amine 6-P and 6-P Gluc in the pentose pathway. How might these changes in glycolytic metabolism be mediated? Collectively, the impairment in PFK activity and glycolysis may reflect the altered maternal metabolism induced by LPD affecting the early embryo and ESC phenotype. LPD causes maternal hyperglycaemia through reduced circulating insulin and raised glucose at the time of blastocyst formation [11,58]. Reduction in insulin signalling in different diabetic models has been shown to disturb glycolysis through reduction in PFK protein but not transcription [59,60]. Moreover, metabolomic screening of adult cardiomyocytes under hypoinsulinaemic conditions revealed increased glycolytic intermediates upstream of PFK  Fig. 2 and text for abbreviations. B. Glycolysis. Changes in relevant RNA expression, metabolite levels, and downstream cellular / physiological phenotype are shown. Red indicates increased and green reduced levels in LPD based upon transcriptomic, metabolomic and cellular studies indicative of bottleneck at PFK and shown by reduced activity. Some metabolites at mean rather than significant levels, see text. See Fig. 3 and text for abbreviations 1 3 activity (F6-P, G6-P) and reduced downstream metabolites (fructose-1,6,bisphosphate [F 1,6 bP], 3-phosphoglycerate ) [59], closely mirroring the phenotype of LPD ESCs. Finally, maternal LPD has been reported to induce metabolic and mitochondrial dysfunction in early mouse embryos [61], a phenotype consistent with our data here. An increase in LPD fatty acids (listed Table 2) may be a further consequence of glycolytic impairment mediated by low insulin whereby maintaining energy supply requires a switch towards fatty acid metabolism as found in type 2 diabetic conditions [62].
Whilst we have focused on the phenotype changes in MAPK signalling and glycolytic metabolism, the RNA sequencing and Paintomics integrated analyses also found changes in other related signalling pathways, notably mammalian target of rapamycin complex 1 (mTORC1) signalling identified as an enriched gene set following GSEA (Supplementary Fig. 2D) and HIPPO signalling following string interaction Paintomics (Supplementary Fig. 4). Substantial crosstalk occurs across several signalling pathways in cells to fine-tune responses to specific environmental cues with, for example, integrated activity across MAPK, mTORC1 and HIPPO signalling following LIF stimulation [63]. Changes in mTORC1 activity commonly occur in DOHaD developmental models to coordinate growth with nutrient availability [64] including the current mouse LPD model [4]. An integrator role for HIPPO signalling has also been recognised coordinating input from metabolism and nutrient levels with cell contact patterns to regulate growth and differentiation [65]. Moreover, HIPPO signalling is directly sensitive to glycolytic activity and energy availability [65,66]. Thus, the phenotype we describe for LPD ESC lines may reflect a combinatorial response across signalling pathways but further studies will be required to confirm.

Conclusion
Our detailed characterisation of LPD ESC lines provides further support for their use in identifying early mechanisms of adverse developmental programming. We highlight two major changes in phenotype with (i) impairment of MAPK signalling associated with reduced cell survival and increased apoptosis and (ii) impairment in upstream glycolytic metabolism mediated in particular by reduced PFK activity. In the sequence of induction of adverse programming mediated through maternal protein restriction, these changes in ESC phenotype derive from the altered maternal metabolite environment in the vicinity of the embryo before implantation. However, the ubiquity of MAPK signalling and energy metabolism in regulating cellular function is such that future studies will be required to identify precisely how these conditions contribute to the step-wise programming of disease phenotype. For example, maternal LPD treatment is known to alter fetal and postnatal growth trajectory which in turn associates with adult disease risk [4] and influences the extent of cellular differentiation and apoptosis during development [9,11] which may reflect altered MAPK signalling. Energy metabolism is also recognised as central to embryo survival and viability and a contributor to DOHaD mechanisms [67]. The combination of these signalling and metabolic effects may also be causally related. For example, in cancer cells, a non-glycolytic role for Pfkp has been identified to alter MAPK signalling activity [68]. These datasets therefore provide the basis for more directed in vivo and in vitro studies to understand mechanisms and preventative approaches to the developmental origin of disease. thank the technical support provided by Metabolon staff (www. metab olon. com ) in the analysis of provided ESC lines. We would also like to thank members of the Centre for Biomarker Research (CeBioR) team at University of Huddersfield.
Author Contribution PK, AC and BI contributed equally and performed the experimental study, analysed the data, wrote draft sections of the manuscript and corrected the final draft. JJE, SWM, JMG provided specific technical and experimental support and edited the paper; PCM lead the RNA sequencing project, performed experiments, analysed qRT-PCR data, supervised researchers and edited the paper. NRS and TPF designed the study, supervised researchers, acquired main funding and wrote the final manuscript.
Funding This work was supported through the Biotechnology and Biological Sciences Research Council, UK (BBSRC; BB/F007450/1) to TPF, the European Union FP7-PEOPLE-2012-ITN EpiHealthNet programme (317146) and FP7-CP-FP Epihealth programme (278418) to TPF, a BBSRC DTA studentship to NRS, TPF and AC; the Centre for Biomarker Research, University of Huddersfield to PCM and BI; the Wessex Medical Research and Rosetrees Trust (M327-CD1) to SWM. For the purpose of open access, the author has applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.

Data Availability
The RNA sequencing dataset has a GEO accession number GSE216108; the metabolomics dataset and other data presented are available from the corresponding author on reasonable request.
Code Availability Not applicable.

Declarations
Ethics Approval Animal procedures were conducted under UK Home Office project licence in accordance with Animal (Scientific Procedures) Act 1986 and local ethics committee at the University of Southampton (UoS).

Consent for Publication
This manuscript has been approved by all authors and is solely the work of the authors named.

3
Conflicts of Interest/Competing Interests The authors have no conflicting or competing interests to declare that are relevant to the content of this article.

Consent to Participate Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.