Effects of low doses of methylmercury (MeHg) exposure on definitive endoderm cell differentiation in human embryonic stem cells

Fetal development is one of the most sensitive windows to methylmercury (MeHg) toxicity. Laboratory and epidemiological studies have shown a dose–response relationship between fetal MeHg exposure and neuro performance in different life stages from infants to adults. In addition, MeHg exposure has been reported to be associated with disorders in endoderm-derived organs, such as morphological changes in liver cells and pancreatic cell dysfunctions. However, the mechanisms of the effects of MeHg on non-neuronal organs or systems, especially during the early development of endoderm-derived organs, remain unclear. Here we determined the effects of low concentrations of MeHg exposure during the differentiation of definitive endoderm (DE) cells from human embryonic stem cells (hESCs). hESCs were exposed to MeHg (0, 10, 100, and 200 nM) that covers the range of Hg concentrations typically found in human maternal blood during DE cell induction. Transcriptomic analysis showed that sub-lethal doses of MeHg exposure could alter global gene expression patterns during hESC to DE cell differentiation, leading to increased expression of endodermal genes/proteins and the over-promotion of endodermal fate, mainly through disrupting calcium homeostasis and generating ROS. Bioinformatic analysis results suggested that MeHg exerts its developmental toxicity mainly by disrupting ribosome biogenesis during early cell lineage differentiation. This disruption could lead to aberrant growth or dysfunctions of the developing endoderm-derived organs, and it may be the underlying mechanism for the observed congenital diseases later in life. Based on the results, we proposed an adverse outcome pathway for the effects of MeHg exposure during human embryonic stem cells to definitive endoderm differentiation. Supplementary Information The online version contains supplementary material available at 10.1007/s00204-023-03580-7.


Introduction
Mercury (Hg) is a well-recognized environmental contaminant that is of global concern (Driscoll et al. 2013). In the environment, Hg (Hg 0 ) occurs in its elemental form, as inorganic salts (Hg + or Hg 2+ ), and as organic compounds such as methylmercury (MeHg) (Clarkson 2002). Anthropogenic effects have increased Hg emission dramatically into the environment (Streets et al. 2009), which is further converted to MeHg through anaerobic microorganisms, typically sulfur-reducing bacteria, in aquatic environments (Clarkson 2002;Parks et al. 2013). Methylmercury can then bioaccumulate along the aquatic food chains into top marine/ freshwater predators, making communities or populations that rely heavily on fish consumption in their diets susceptible to MeHg exposure (Driscoll et al. 2013;Van Oostdam et al. 2005). Also, rice which grows in moist, anaerobic soil that is ideal for microbial methylation of mercury, can also accumulate MeHg, and affect the rice-consuming populations (Feng et al. 2021).
The main target site for MeHg is the central nervous system, with clinical findings showing MeHg-induced neurobehavioral and cognitive function impairments in the victims of the Japanese and Iraqi MeHg poisoning events (Amin-Zaki et al. 1981;Bakir et al. 1973;Harada 1978Harada , 1995. Animal studies under well-controlled laboratory conditions have also confirmed the neurotoxic effects of MeHg, with the developing central nervous systems being especially susceptible, causing more severe effects in prenatal exposure compared to exposure in later life stages (Castoldi et al. 2008). This is in part attributed to MeHg-induced negative effects on neuronal migration and differentiation (Bose et al. 2012;Fujimura and Usuki 2015;Tamm et al. 2006), and also MeHg-induced cell death due to distortion of intracellular calcium and glutamate homeostasis, oxidative stress generation, as well as mitochondrial dysfunction (dos Santos et al. 2016;Nagashima 1997).
In addition to its well-characterized neurotoxic effects, MeHg exposure can also lead to adverse effects in tissues other than the brain. In victims of the MeHg poisoning event in Japan, erosive inflammation of the digestive tract, fatty degeneration of the liver and kidney, and disturbance of pancreatic islet cells were all observed in addition to neuronal disturbances (Eto 1997;Eto et al. 2010;Harada 1995). Furthermore, the hepatic disorder was also identified as a cause of death for MeHg poisoning victims, with comparable levels of MeHg accumulation in the liver compared to the brain (Harada 1995). Multiple cohort studies have shown a positive correlation between mercury exposure levels and elevated insulin levels accompanied by decreased pancreatic β-cell function (Chang et al. 2011;He et al. 2013). Animal studies also have corroborated epidemiology findings, showing that MeHg accumulated within the liver to comparable levels of the brain (Carneiro et al. 2014;Gonzalez et al. 2005), resulting in imbalanced redox, hepatic glycogen accumulation, disruption of energy metabolism pathways, as well as alteration of global gene expression profiles within the liver (da Rosa-Silva et al. 2020;Mela et al. 2007;Ung et al. 2010;Yadetie et al. 2013). Apoptosis and dysfunction in pancreatic β-cells through oxidative stress were also observed (Chen et al. 2006;Schumacher and Abbott 2017;Yang et al. 2022). In vitro studies on non-brain derived cell lines exposed to MeHg, though limited, have shown that MeHg-induced oxidative stress resulted in apoptosis in liverderived HepG2 cells (Cuello et al. 2010) and altered cell morphology in human-induced pluripotent stem-cell-derived hepatocytes (Jamalpoor et al. 2022). Despite these observations, the mechanism of the effects of MeHg exposure on non-brain-related organs is still poorly understood compared to the effects of MeHg on the brain, and studies examining the developmental effects of MeHg exposure on non-brain related organs is still lacking.
In recent years, there is a growing trend of reducing the use of animals for testing of hazardous compounds, and in vitro toxicity testing that relies on cultured cells has emerged as an alternative tool for toxicity assessment.
Among these, embryonic stem cells, and especially human embryonic stem cells (hESCs), which have the potential to differentiate into all fetal cell lineages represent an attractive cellular system for in vitro studies in developmental toxicology (Visan et al. 2012), and has been utilized to examine the embryotoxicity of MeHg Seiler and Spielmann 2011;Stummann et al. 2007). As the effects of MeHg on neuronal differentiation have been extensively studied Prince et al. 2021;Stummann et al. 2009;Tamm et al. 2008;Theunissen et al. 2011;Zimmer et al. 2011), our objective is to examine the developmental toxicity of MeHg on other non-neural-derived organs using hESCs as a model. Since the potentially affected organs such as the liver, pancreas, and the digestive tract (Eto 1997;Eto et al. 2010;Harada 1995) are endoderm derived, and the endoderm is the first germ layer to develop out of the three germ layers (Muhr and Ackerman 2020), we hypothesized that MeHg exposure affects the early differentiation of the endoderm in the hESCs. Using a transcriptomics approach, we aim to not only characterize the effects of MeHg during definitive endoderm formation but also explore the potential cellular mechanisms that lead to developmental toxicity due to MeHg exposure. We hypothesized that MeHg could disrupt the differentiation of endoderm formation in human stem cells by distortion of intracellular calcium and glutamate homeostasis, oxidative stress generation, as well as mitochondrial dysfunction.

Ethics review and approval
This project was reviewed and approved by the Health Canada and Public Health Agency of Canada's Research Ethics Board (File No. REB 2016-027H) and by the Office of Research Ethics and Integrity of the University of Ottawa .

Cell maintenance and induction
Human embryonic stem cells (H9; passage 30; WiCell Research Institute) were maintained according to the protocols established in our previous study . Briefly, hESCs were maintained and expanded in six-well plates, which were pre-coated with Matrigel (Catalog no. 354 230; BD Biosciences). Cells were cultured under 37 °C, 4% oxygen (O 2 ) and 10% carbon dioxide (CO 2 ) in a tri-gas incubator (Thermo Fisher Scientific) and fresh Essential 8 Flex Medium (Catalog no. A2858501; Thermo Fisher Scientific) was changed daily. On the first day of seeding, 10 µM Rho-associated kinase inhibitor (Catalog no. 72304; Stem-Cell Technologies) was added to the medium to enhance the survival of the cells. Cells were gently detached with Gibco StemPro Accutase cell dissociation reagent (Catalog no. A1110501; Thermo Fisher Scientific) and seeded onto new plates for subsequent experiments when the confluency reached around 85%. Only hESCs under passage 40 (within 10 passages from purchase) were used to avoid any serious genetic instability associated with prolonged passaging.
The hESCs were seeded in six-well plates at a density of 2 × 10 4 cells/mL with the growth medium (Essential 8 Flex Medium) and incubated at 37 °C, 4% oxygen (O 2 ) and 10% carbon dioxide (CO 2 ) with daily medium change. When the cell confluency reached around 20%-30%, the growth medium was replaced with the definitive endoderm induction medium A from Gibco ™ PSC Definitive Endoderm Induction Kit (A3062601, Thermo Fisher Scientific). 24 h later, the culture medium was disposed of, and the definitive endoderm induction medium B was used for another 24 h of incubation. The hESC-derived DE cells were verified by immunostaining for both self-renewal (OCT4) and endodermal protein markers (SOX17 and CXCR4) (Fig. S1). In brief, the cells were washed twice with phosphate-buffered saline (PBS) after disposing the induction medium B and fixing with 4% paraformaldehyde (Santa Cruz Biotechnology). The cells were then blocked and permeabilized with 5% bovine serum albumin (BSA) diluted in PBS containing 0.3% Triton X-100 (Sigma-Aldrich). The cells were, respectively, incubated with rabbit anti-human OCT4 antibody (Catalog no. C30A3; Cell Signaling Technology), goat antihuman SOX17 antibody (Catalog no. AF1924; R&D systems), and rabbit anti-human CXCR4 antibody (Catalog no. ab124824; Abcam) diluted at 1:200 in 5% BSA, and incubated overnight at 4 °C. The cells were then gently washed with PBS for three times and, respectively, incubated with Alexa Fluor 488 goat anti-rabbit IgG (Catalog no. A11008; Thermo Fisher Scientific), and Alexa Fluor 594 donkey antigoat IgG (Catalog no. A11058; Thermo Fisher Scientific) diluted at 1:500 in 5% BSA at 4 °C, overnight. Prolong Gold antifade reagent containing 4',6-diamidino-2-phenylindole (DAPI; Catalog no. S36939; Thermo Fisher Scientific) was added and the immunofluorescent images were taken with a Nikon A1RsiMP confocal microscope (Nikon Instruments Inc.) with Nikon's Imaging Software NIS-Elements.

Chemical preparation and exposure
MeHg chloride (purity: 99.9%) was obtained from Sigma-Aldrich and dissolved in dimethyl sulfoxide (DMSO, Sigma) to obtain 1000 × stock solutions (i.e., 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 1, 5 and 10 mM). hESCs were seeded in sixwell plates and allowed for growth until they reached a confluency of around 20%-30%. The cells were then exposed to 0 (DMSO as vehicle control), 0.001, 0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 1, 5 or 10 µM MeHg by 1000 × dilution of the stock solutions in induction medium for the initial cell viability assay. The experimental doses were narrowed down according to the results of cell viability and the cells were treated with 0 (DMSO as vehicle control), 10, 100 or 200 nM MeHg by 1000 × dilution of the stock solutions in induction medium for the rest of assays. The cells were repeatedly dosed each day throughout the induction process.

Cell viability and morphology
Cell viability was determined by a Cell Counting Kit 8 (WST-8/CCK8) (Catalog no. ab228554, Abcam Inc.) according to the manufacturer's instructions. In brief, hESCs were seeded, induced, and exposed to MeHg as described above. At the end of induction and exposure, 200 μL of WST-8 solution was added into each well on the six-well culture plate followed by 10 s of shaking and incubation at 37 °C for 1.5 h in the dark. The absorbance at 460 nm wavelength was read on a BioTek Cytation 3 Cell Imaging Multi-Mode Reader (Fisher Scientific). The effects of 0, 10, 100, and 200 nM MeHg on cell viability were also checked by cell counting. Briefly, hESCs were seeded, induced, and exposed to MeHg as described above. At the end of induction and exposure, the medium was disposed, and cells were gently lifted by StemPro Accutase cell dissociation reagent after rinsing by PBS. Then the cells were centrifuged and re-suspended in 1 ml of PBS and gently pipetted to generate a single-cell suspension. 10 µL of the cell suspension was taken for cell counting, and the cell numbers were counted with a Countess Automated Cell Counter (Invitrogen, Thermo Fisher Scientific) after 1:1 dilution of cell suspension in Trypan blue solution (Catalog no. T10282; ThermoFisher Scientific).
For morphology observations, hESCs were seeded, induced to DE, and treated with MeHg, as described above. The cell morphology was observed and imaged at each time point (before induction, during induction and after induction) with a 10 × objective on a Zeiss Axiovert 40 CFL inverted microscope (Carl Zeiss Microscopy LLC) with a SPOT RT3 digital camera and SPOT basic software (Diagnostic Instruments Inc.).

RNA sample preparation and next-generation sequencing
MeHg doses used for RNA sequencing were chosen based on cell viability and morphology results, as well as a literature review on the Hg concentrations found in maternal blood (Table 1). The final doses of 0, 10, 100, and 200 nM of MeHg exposure were chosen for RNA sequencing as these doses did not induce significant cell death and covered the range of Hg concentrations typically found within human blood. hESCs were cultured, induced, and treated as described above. RNA samples were extracted and purified with a RNeasy Plus Mini kit (Catalog no. 74134, Qiagen) following the manufacturer's protocol. In brief, the culture medium was completely removed, and cells were lysed directly with Buffer RLT Plus containing 10% β-mercaptoethanol (β-ME). The cell lysates were gently pipetted into microcentrifuge tubes and mixed by vortexing. Lysates were homogenized by centrifugation for 2 min at maximum speed (14,800 rpm) in QIAshredder spin columns placed in 2 ml collection tubes and were transferred to gDNA eliminator spin columns for gDNA removal. After centrifuging for 30 s at 10,000 rpm, 70% ethanol was added to the flow-through, and mixed well by pipetting. The mixed sample solutions were transferred to RNeasy spin columns and centrifuged for 15 s at 10,000 rpm. The flow-through was discarded, and spin column membranes were washed by centrifuging for 15 s at 10,000 rpm in Buffer RW1 and discarding the flow-through. The spin column membranes were then washed by centrifuging for 15 s at 10,000 rpm in Buffer RPE, followed by a second wash in Buffer RPE for 2 min at 10,000 rpm. The spin columns were carefully placed in new collection tubes and centrifuged for 1 min at maximum speed (14,800 rpm) to eliminate any possible residual Buffer RPE. To elute the RNA samples, the spin columns were placed in 1.5 mL collection tubes and 30 μL RNase-free water was added directly to the column membrane. The collection tubes were centrifuged for 1 min, and the eluted RNA samples were stored immediately at − 80 °C.
The RNA samples were sent to Genome Quebec for quality check and whole genome-scale mRNA sequencing by Illumina NovaSeq 6000. Sequencing was performed by paired end reading with a sequencing depth of 50 M reads for each sample and three biological replicates for each treatment group.

Sequencing data analysis and marker selection for verification
The raw data were obtained from Genome Quebec, and quality control was conducted with FastQC (https:// www. bioin forma tics. babra ham. ac. uk/ proje cts/ fastqc/). All analyses were conducted in R (https:// www.r-proje ct. org/). Gene sequences were aligned to the reference human genome database GRCh38 (https:// www. ncbi. nlm. nih. gov/ data-hub/ genome/ GCF_ 00000 1405. 26/), and reads counted using the "Rsubread" package (Liao et al. 2019). Data were then normalized, and lowly expressed genes were filtered out using the "edgeR" package (Robinson et al. 2010). Boxplots showed filtered data to reflect quantile normalization, and principal component analysis (PCA) was used to demonstrate the clustering of samples in different groups. The R package "edgeR" (Robinson et al. 2010) was used to identify differentially expressed genes (DEGs) with cutoff for DEGs set as adjusted p value < 0.05 and absolute log2 [fold change (FC)] > 1. Significant DEGs were further visualized in volcano plots and heat maps generated with the "ggplot2" package.
Gene set testing for DEGs was performed with the "clus-terProfiler" package (Yu et al. 2012) for Gene Ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Adjusted p value < 0.05 and count ≥ 30 was set as the cutoff for GO functional enrichment, whereas adjusted p value < 0.05 was used as the cutoff for KEGG pathway analysis. Data were visualized using the "ggplot2" package.
GO functional enrichment analysis identified 32 GO terms that were significant, of which 7 were related to early development. A protein-protein interaction network analysis of the genes within each of these seven early developmentrelated GO terms was conducted using Cytoscape (https:// cytos cape. org/) to identify hub genes for each of these GO 3.0 (1.7-11) 15.0 (8.5-54.8) Rothenberg et al. (2013) terms. Lineage genes classified based on the TaqMan hPSC scorecard assay (Tsankov et al. 2015) were also identified in the dataset, and a cross comparison of hub genes identified from GO terms and DEGs from lineage genes with genes driving endoderm differentiation in the literature (Aksoy et al. 2014;Chu et al. 2016;Faial et al. 2015;Teo et al. 2011;Xuan and Sussel 2016) narrowed down genes for further downstream qPCR and western blot verification. We explored the potential mechanisms of the developmental toxicity of MeHg on hESCs while differentiating toward DE cells by identifying the sequence of responses to the increasing dose of MeHg. We entered the DEGs into a web-based tool, FastBMD (https:// www. fastb md. ca/) (Ewald et al. 2021), to calculate the benchmark doses (BMD). The DEGs were then ranked by their BMDs, and the enriched KEGG pathways were ranked based on the average rank of DEGs involved in each pathway (cutoff of involved DEG number ≥ 5).

Reverse transcription quantitative polymerase chain reaction (RT-qPCR)
RNA samples were prepared as described above. RNA concentration was determined using NanoDrop ™ 2000 spectrophotometer (Thermo Scientific) and RNA quality was verified using gel electrophoresis and qualitative assessment of the relative intensities of the 28 and 18 s ribosomal RNA bands. RNA was treated with DNaseI (Catalog no. 18068015, Invitrogen, ThermoFisher Scientific) prior to cDNA synthesis from 900 ng total RNA using a high-capacity cDNA reverse transcription kit (Catalog no. 4368814; ThermoFisher Scientific). Real-time PCR was performed in 20 µl reactions containing 10 µl SsoFast EvaGreen supermix (Catalog no. 1725201, Bio-Rad), 500 nM forward and reverse primers (Table S1), and 0.5 µl cDNA template obtained from cDNA synthesis. Controls (no template added) were run in every assay, and non-reverse transcribed samples were run for every primer pair. Reaction efficiency for all primer pairs was between 90 and 110%. Final data were normalized to the geometric mean of the three reference genes (gapdh, hddc2 and znf324b) and the 0 nM group using the delta-delta CT method (Rao et al. 2013).

Western blotting
To determine whether the effects of MeHg on gene expressions also occur at the protein level, we measured the protein expression of the selected markers by western blots.
hESCs were cultured, induced, and exposed to MeHg as described above. Culture medium was discarded, and the cells were washed with PBS and lysed with Pierce radioimmunoprecipitation assay buffer (RIPA buffer; Thermo Fisher Scientific) containing proteinase and phosphatase inhibitor cocktails (Catalog no. 5872S; Cell Signaling Technology). The cell lysates were then sonicated with a Microson Ultrasonic Cell Disruptor XL2000 (Misonix Inc.). The whole process of protein extraction was conducted on ice. Protein concentrations of the lysates were quantified with a Pierce Bicinchoninic Acid Protein Assay kit (Thermo Fisher Scientific). The cell lysates were then denatured by heating at 100 °C for 4 min after mixing with 4 × Laemmli sample buffer (Catalog no. 1610747; Bio-Rad).
Cell lysate samples were electrophoresed on 4-20% Mini-Protein TGX Stain-Free Precast Gels (Catalog no. 17000436, Bio-Rad). The same amount of protein was loaded for each batch, and the amount of total protein loaded ranged from 20 to 24 μg between different batches. Gels were then activated under ultraviolet light using the Chemi-Doc XRS + imaging system (Bio-Rad), and proteins were electro-transferred from gels to polyvinylidene difluoride (PVDF) membranes with the Mini Trans-Blot Cell (Catalog no. 1703930; Bio-Rad). Membranes were then blocked with 5% non-fat dry milk (NFDM, Catalog no. 1706404; Bio-Rad) or 3% bovine serum albumin (BSA, Catalog no. BSA-50; Rockland Immunochemicals, Inc.) depending on the protein targets for 1 h at room temperature and incubated with antibodies of interest (Table S2) diluted in their corresponding blocking buffer at 4 °C overnight. Then the membranes were washed three times with Tris-buffered saline containing 0.1% Tween 20 (TBST) on a shaker for 5 min and incubated with 1:5,000 dilution of horseradish peroxidase (HRP)-linked secondary antibodies (Table S2) diluted in corresponding blocking buffer at 4 °C overnight. The membranes were rewashed for three times with TBST and incubated with enhanced chemiluminescence substrate (ECL, Catalog no. 1705062; Bio-Rad) in the dark and imaged with the ChemiDoc XRS + imaging system. The Image Lab software (Bio-Rad) was used to quantify the band intensities, and all target band intensities were normalized against the total protein band intensities.
As autophagy plays a critical role in the differentiation of stem cells (Sharma et al. 2022), the expression of autophagic protein markers (LC3B and SQSTM1, antibody information is shown in Table S2) was also measured by western blots as described above.

Statistical analysis
All data were displayed as the means plus or minus the standard errors of the mean (means ± SEMs) and were tested for statistical significance with one-way ANOVA, followed by Dunnett's post hoc tests. The differences were considered significant if p < 0.05. Except for the RNA sequencing data, which was analyzed with R as described above, all statistical analyses were performed using GraphPad Prism (version 8).

Effects on morphology and cell viability of MeHg exposure during definitive endoderm (DE) differentiation
To determine sub-lethal doses of MeHg exposure during hESC to DE differentiation, cell viability was measured following 0-10 μM of MeHg exposure. Differentiating cells exposed to 300 nM of MeHg showed a significant reduction in cell viability, and exposure to higher than 1 μM of MeHg resulted in almost complete cell death (Fig. 1). Exposure to concentrations at or below 200 nM of MeHg did not have a significant impact on cell viability, and differentiated DE cells showed no morphological difference when compared to the controls (Fig. 2). Therefore, the doses of 0, 10, 100, and 200 nM MeHg were identified as the sub-lethal doses for further transcriptomic analysis.

Identification of differentially expressed genes (DEGs)
Differentially expressed genes were identified with the "edgeR" package following data cleanup and quantile normalization of the samples (Fig. 3a, b). Based on the cutoff for DEGs (adjusted p value < 0.05 and absolute log2 FC > 1), a total of 1163 DEGs were identified (Excel Table S1), with 26 DEGs in the 10 nM treatment group, 393 DEGs in the 100 nM treatment group, and 1136 DEGs in the 200 nM treatment group. A total of 16 genes were differentially expressed in all 3 treatment groups. These DEGs were further visualized via a heat map (Fig. 3c) and a volcano plot (Fig. 3d-f). The heat map showed a clear clustering of each MeHg treatment group.

Functional analysis of DEGs
The 1163 DEGs were analyzed for GO functional enrichment with a cutoff of p adj < 0.05. A total of 87 GO terms were identified to be significant, and among those, 32 GO terms with greater than 30 count number were kept for further analysis (Fig. 4a). The clustering of these 32 significant GO terms was organized into an enrichment map (Fig. 4b), showing that many of the significant GO terms were related to RNA metabolic processes including ribonucleoprotein complex biogenesis (GO:0022613), ribosome biogenesis (GO:0042254), ncRNA metabolic process (GO:0034660), ncRNA processing (GO:0034470), rRNA processing (GO:0006364), rRNA metabolic process (GO:0016070), and nuclear-transcribed mRNA catabolic process, nonsense-mediated decay (GO:0000184) or protein synthesis/ transport such as mRNA catabolic process (GO:0006402), protein targeting (GO:0006605), establishment of protein localization to membrane (GO:0090150), translation initiation (GO:0006413), and protein targeting to membrane (GO:0006612). In addition, seven of these GO terms, including negative regulation of cell differentiation (GO:0045596), gland development (GO: 0048732), urogenital system development (GO: 0001655), stem cell differentiation (GO: 0048863), renal system development (GO:0072001), mesenchyme development (GO:0060485), and kidney development (GO:0001822) were GO terms that were linked to the GO term developmental process (GO:0032502), and the genes that are linked to these seven GO terms are shown in Fig. 4c. To select gene markers of interest for verification, a protein-protein interaction (PPI) network analysis of the genes within each of these seven development-related GO terms was conducted, and the top ten hub genes were identified from each of these seven GO terms. A total of 27 hub genes were identified after accounting for repetitions, with BMP4 being a hub gene in all 7 GO terms, and BMP2, WNT4, and SMAD3 appearing as hub genes in six out of the seven GO terms (Excel Table S2).
The KEGG pathway enrichment analysis (Fig. 5) also yielded similar results as the GO functional enrichment, with multiple ribosome-related pathways such as ribosome (hsa03010), ribosome biogenesis in eukaryotes (hsa03008), spliceosome (hsa03040), and RNA polymerase (hsa03020) being upregulated in all three MeHg treatment groups (Fig. 5a-c). In addition, pathways that are known to be perturbed by MeHg (Farina et al. 2011;Grotto et al. 2009;Novo et al. 2021;Roos et al. 2012), such as calcium signaling pathway (hsa04020), nucleotide excision repair (hsa03420), Fig. 1 Cell viability (as percentage of control) of definitive endoderm (DE) cells differentiated from human embryonic stem cells (hESCs) exposed to 1 nM to 10 µM MeHg during DE differentiation. One-way ANOVA and Dunnett's post hoc test were used at each time point to compare differences between treatment and 0 nM group with p < 0.05 being considered as significant. Data are presented as means ± SEMs with N = 6-9 for each group and chemical carcinogenesis-reactive oxygen species (hsa05208), are also reflected in the KEGG pathway enrichment analysis (Fig. 5e). In addition, according to the average rank of DEGs based on their BMDs (Excel Table S3), 21 KEGG pathways were successfully ranked with protein digestion and absorption (hsa04974) and calcium signaling pathway (hsa04020) considered as the early responses triggered by increasing doses of MeHg (Excel Table S4).
Lineage marker genes within DEGs were also identified based on previously published lineage genes (Bock et al. 2011;Tsankov et al. 2015), and ten lineage genes were found to be differentially expressed (Fig. 6). These included mesodermal/endodermal genes PDGFRA, FOXA2, GATA6, CLDN1, EOMES, and LEFTY2, ectodermal genes COL2A1 and MAP2, and self-renewal genes DNMT3B and SOX2. Of these genes, CLDN1, LEFTY2, and COL2A1 were downregulated following 200 nM of MeHg exposure during DE differentiation, whereas all the other lineage genes were upregulated.
As BMP4, SMAD3, FOXA2, GATA6, EOMES, and LEFTY2 were also major genes involved in the endoderm differentiation (Aksoy et al. 2014;Chu et al. 2016;Faial et al. 2015;Teo et al. 2011;Xuan and Sussel 2016), these genes were chosen for further downstream qPCR and western blot verification. In addition to the above genes, two other genes BHMT and URB1 were also selected for downstream verification. BHMT was chosen based on the justification that it has the biggest fold change among all DEGs, and BHMT itself is highly expressed in the endoderm-derived organs such as the liver and pancreas (Pajares and Pérez-Sala 2006), and URB1 was chosen as it is a ribosome biogenesis protein that can modulate endoderm-derived organ development (He et al. 2017), and ribosomal pathways were top enriched in both GO functional analysis and KEGG pathway analysis. Fig. 2 Cell attachment, colony formation, and morphology of definitive endoderm (DE) cells differentiated from human embryonic stem cells (hESCs) exposed to MeHg during DE differentiation. a Representative images of hESCs before induction. b-e Representative images of DE cells exposed to 0 nM, 10 nM, 100 nM, and 200 nM MeHg during DE cell differentiation. All images were taken under a 10× objective

RT-qPCR and western blot verification
Gene expression of the above-mentioned genes was further verified using RT-qPCR. Significant downregulation of bmp4 (Fig. 7a) and smad3 (Fig. 7b) was observed in the 200 nM MeHg treatment group, showing the same trends as the transcriptomics dataset. However, the other four genes [eomes, foxa2, gata6, urb1 (Fig. 7c-f)] tested, which showed a significant but low fold change difference in the transcriptomics dataset, did not show any statistical differences when analyzed using RT-qPCR.

Discussion
The main goal of the present study was to assess the influence of MeHg exposure on the development of human embryonic stem cells (hESCs) during definitive endoderm (DE) cell differentiation using a transcriptomics approach. Our results showed that MeHg targeted mainly ribosomal processes to alter protein synthesis and translocation during hESCs to DE cell differentiation, resulting in an overall enhancement of cell differentiation.
The viability of differentiating cells was not affected when exposed to MeHg doses below or at 200 nM, whereas MeHg doses above 200 nM decreased its viability in a dosedependent manner. In contrast, undifferentiated hESCs showed close to 50% reduction in cell viability when exposed to 200 nM of MeHg for the same duration , suggesting that the hESCs which are differentiating toward DE cells have a higher tolerance for MeHg exposure in terms of cell viability when compared to undifferentiated hESCs. However, this does not indicate that MeHg exposure under 200 nM has no effect on these differentiating hESCs, Fig. 4 Significant enrichment GO terms of differentially expressed genes (DEGs). a GO functional enrichment dot plots of DEGs. The dots size represents the number of genes among the significant DEGs associated with the GO term. The color of the dots represents their significance in accordance with the adjusted p value. b Enrichment map of the 32 GO terms identified. c Gene-concept network plot of the seven GO terms related to early embryo development as sub-lethal doses of MeHg exposure have been shown to alter protein expression and induce cognitive deficits in animal models (Bittencourt et al. 2019). This is supported by findings in this study with respect to the effects of MeHg at sub-lethal doses on the altered gene and protein expressions of hESCs during the differentiation toward DE cells.
The transcriptomics dataset showed that MeHg exposure at 10 nM affected the least number of genes, with only 26 genes showing a significant difference in expression pattern compared to the vehicle control. The number of differentially expressed genes (DEGs) increased as MeHg concentrations increased, but the genes that were affected remained mostly consistent, as 85% of DEGs in the 10 nM MeHg-dosed group and 93% of DEGs in the 100 nM MeHg-dosed group were also DEGs in the 200 nM MeHg-dosed group, suggesting a generally similar manner in gene pattern alterations among the three doses. Functional analysis of DEGs via GO functional enrichment or KEGG pathway analysis also confirmed that MeHg disrupted hESCs differentiation into DE cells via similar pathways, with ribosomal and protein MeHg-dosed group compared to the controls. e Enrichment plot of KEGG pathways known to be perturbed by MeHg in 200 nM MeHgdosed group compared to the controls synthesis/translocation pathways consistently appearing in all three MeHg-dosed groups. Thus, it is likely that MeHg exerts its developmental toxicity mainly by disrupting ribosome biogenesis during early cell lineage differentiation.
It is well-recognized that ribosome biogenesis is important for the maintenance of stem cell pluripotency and its subsequent differentiation (Gabut et al. 2020;Saba et al. 2021;Tahmasebi et al. 2019). In undifferentiated stem cells, ribosome biogenesis is generally at a high level, though global protein synthesis rates tend to be low, which primes the undifferentiated stem cell for rapid and efficient reassembly of their proteome into the target cell types during differentiation. During differentiation, rates of protein synthesis and ribosome biogenesis are tightly and dynamically regulated in accordance with the needs of the differentiating cells, with protein synthesis rates increasing during differentiation. In contrast, ribosome biogenesis transiently decreases only during the early phases of differentiation and then rises back up to high levels (Baser et al. 2019;Corsini et al. 2018). Following MeHg exposure, upregulation of genes within ribosomal-related processes enriched from GO terms, upregulation of the ribosome biogenesis homolog Fig. 6 Lineage genes identified within the transcriptomics data set and the number of lineage genes that were differentially expressed following MeHg exposure during human embryonic stem cell to definitive endoderm differentiation Fig. 7 Expression of bmp4, smad3, eomes, foxa2, gata6, and urb1 genes in definitive endoderm (DE) cells exposed to MeHg during human embryonic stem cells to DE cell differentiation. Data are presented as means ± SEMs with N = 9 for each gene. The small circle symbols represent individual values of each experiment. One-way ANOVA and Dunnett's post hoc test were used at each time point to compare differences between treatment and control groups with p < 0.05 being considered as significant marked by asterisk URB1, and the upregulation of the KEGG pathway ribosome biogenesis in eukaryotes (hsa03008) all suggest that MeHg upregulates ribosomal biogenesis during hESC to DE cell differentiation. Despite undifferentiated stem cells having high ribosome biogenesis, and MeHg exposure on undifferentiated hESCs showing enhanced pluripotency , it is not likely that MeHg exposure during hESC to DE cell differentiation inhibited differentiation and enhanced hESCs' pluripotency. Firstl, despite "negative regulation of cell differentiation" being enriched through GO functional enrichment, most genes within the GO term were downregulated following MeHg exposure (39 downregulated vs 27 upregulated), suggesting inhibition of negative regulation, and thus enhancement of cell differentiation. Second, multiple pieces of evidence within the dataset also point to an upregulation of protein synthesis, translocation, and turnover. This includes genes within the GO terms "translational inhibition" (seven upregulated and none downregulated) and "protein targeting" (nine upregulated and six downregulated) being mostly upregulated, and core enriched genes within the KEGG pathways for "spliceosome", "RNA polymerase", and "aminoacyl-tRNA biosynthesis" being upregulated. In addition, MeHg is capable of inducing protein synthesis in animal models, with increases in levels of liver ribosomes, ribosomal subunits, and polyribosomes, along with a threefold increase in the incorporation of labeled amino acids in the livers of rats exposed to MeHg (Brubaker et al. 1971). Thus, the combination of ribosome biosynthesis and protein synthesis upregulation suggests that MeHg upregulates ribosome biogenesis leading to increased global protein translation and enhancing differentiation which could further lead to the overgrowth of specific organs. For example, in zebrafish, MeHg caused upregulation of URB1 during DE differentiation, leading to alterations in ribosome biogenesis and overgrowth of digestive organs (He et al. 2017).
What is seemingly contradictory to this conclusion is that nine proteasome subunits (PSMA5, PSMA7, PSMB3, PSMB6, PSMD1, PSMD3, PSMD8, PSMD11, PSME3) and the KEGG pathway proteasome (hsa03050) were also all upregulated following MeHg exposure during hESC to DE cell differentiation. High levels of proteasome activity have been associated with the maintenance of pluripotency (Saez et al. 2018;Schröter and Adjaye 2014;Yan et al. 2020). However, studies have also shown that the ubiquitin proteasome pathway plays a role in modulating the toxicity of MeHg. Specifically, overexpression of the ubiquitin-protein ligase CDC34 significantly increased cellular-ubiquitinated protein levels in yeast and human cell lines, increasing their resistance to MeHg exposure (Hwang 2007;Hwang et al. 2002). Thus, the upregulation of the proteasome pathway Fig. 8 Expression of BHMT, p-SMAD3, URB1, FOXA2, BMP4, EOMES, GATA6, LEFTY2, and SMAD3 protein in definitive endoderm (DE) cells exposed to MeHg during human embryonic stem cells to DE cell differentiation. Data are presented as means ± SEMs with N = 4-6 for each protein. Each batch was first normalized to 0 nM exposure group before comparison. The small circle symbols represent individual values of each experiment. One-way ANOVA and Dunnett's post hoc test were used at each time point to compare differences between treatment and control groups with p < 0.05 being considered as significant marked by asterisk may be due to cells resisting MeHg toxicity and not maintaining pluripotency.
Our results in the ranked KEGG pathways based on BMDs of the gene set involved showed that the calcium signaling pathway (hsa04020) is the most sensitive responses to MeHg exposure. This indicates that MeHg toxicity began with the disruption of calcium homeostasis which can result in the generation of excessive amounts of ROS (Farina and Aschner 2017). This is nicely supported in our data with KEGG pathways chemical carcinogenesisreactive oxygen species (hsa05208) and oxidative phosphorylation (hsa00190) ranking behind the KEGG pathway calcium signaling pathway. Proteasome (hsa03050), which ranked behind chemical carcinogenesis-reactive oxygen species and oxidative phosphorylation, would then need to be upregulated to clean up the damaged proteins caused by the elevated oxidative stress. Meanwhile, the overproduction of ROS has also been linked to hyperactive ribosome biogenesis, and protein synthesis in animal embryos and human cell lines (Huang et al. 2021;Oliveira et al. 2019), and is considered to promote stem cell differentiation (Bigarella et al. 2014), including toward the DE fate (Lv et al. 2022).
Under normal circumstances, antioxidants such as glutathione can prevent oxidative stress through the removal of ROS. However, glutathione is also involved in Hg detoxification by specifically binding with MeHg to form a complex that prevents Hg from binding to cellular proteins and causing cellular damage (Kromidas et al. 1990) or by forming a glutathione-mercury complex in the liver for the excretion of mercury (Zalups 2000). Therefore, depleting existing glutathione concentrations will result in a decreased capacity for the removal of ROS and subsequent cell death. To respond to the depletion of glutathione by Hg, upregulation of the trans-sulfuration pathway that synthesizes glutathione from homocysteine (Hcy) has been observed (Woods and Ellis 1995). From our results, this upregulation of glutathione synthesis is reflected as a downregulation of betaine-homocysteine S-methyltransferase (BHMT), one of the genes with the biggest fold change following MeHg exposure. BHMT normally plays an important role in the regulation of Hcy metabolism by catalyzing the resynthesis of methionine from Hcy (Pajares and Pérez-Sala 2006). Thus, the downregulation of BHMT would facilitate Hcy in entering the trans-sulfuration pathway for glutathione synthesis. Interestingly, folate biosynthesis (hsa00790), which produces folate, that also catalyzes the resynthesis of methionine from Hcy was significantly upregulated following MeHg exposure. This could partly reflect the lack of methionine within the cell due to Hcy being directed for glutathione synthesis to protect against oxidative stress induced by MeHg exposure.
As MeHg resulted in enhanced differentiation during hESC to DE cell differentiation, it is likely that endodermal fate was over-promoted. Under normal conditions, three major extracellular signals, Wnt, Nodal/Activin, and BMPs, regulate definitive endoderm differentiation. These signals all, in turn, phosphorylate SMAD2/3 proteins (Faial et al. 2015;Funa et al. 2015;Teo et al. 2011), which depending on the transcriptional regulator, it binds to either activate endodermal fate genes when bound to SMAD4 (Yang and Jiang 2020) or activate mesodermal fate genes when bound to FOXH1 (Charney et al. 2017). In addition, BMP4 can act synergistically with activin to further upregulate the expression of transcription factors involved in endoderm differentiation and promote the endodermal fate (Teo et al. 2012). BMP4 itself can also cooperate with FGF2 via ERK to further induce mesoderm and inhibit the endoderm differentiation (Bernardo et al. 2011). To understand how MeHg may have disrupted stem cell specification during the differentiation toward DE cells, we examined lineage gene markers within our dataset. We found that the most significantly affected genes by MeHg were mesodermal or endodermal genes, including pdgfra, foxa2, gata6, cldn1, eomes and lefty2. Most notably, the expression of endodermal genes foxa2, eomes, and gata6 were upregulated, and FOXA2 was confirmed to be upregulated at the protein level as well. FOXA2 is an important transcription factor that is expressed during early endoderm differentiation and, when overexpressed, can lead to the expression of genes associated with endodermal lineage, but does not induce the expression of genes involved in late endoderm differentiation (Levinson-Dushnik and Benvenisty 1997). EOMES, on the other hand, interacts with SMAD2/3 to initiate the transcriptional network governing endoderm formation and, when overexpressed, leads to the expression of DE markers such as FOXA2 (Faial et al. 2015;Teo et al. 2011). These changes are likely initiated by the upregulation of phosphorylated SMAD3 (p-SMAD3) though smad3 gene expression was significantly downregulated in this study, suggestive of a negative feedback loop between p-SMAD3 and smad3 gene expression. Like smad3 gene and p-SMAD3 protein, there might be a negative feedback loop for BMP4 regulation as well, as bmp4 gene was downregulated, yet BMP4 protein showed an increasing though non-significant trend following MeHg exposure, which could further promote endoderm differentiation. Thus, MeHg exposure during hESC to DE cell differentiation leads to overexpression of endodermal genes, which may ultimately result in the over-promotion of endodermal fate (Fig. 9).
In summary, our study showed that MeHg exposure at environmentally relevant doses during hESC to DE differentiation has sub-lethal effects that could cause global gene expression changes, leading to calcium dyshomeostasis and ROS overproduction, resulting in an upregulation of ribosome biogenesis and global protein synthesis/ translocation to promote hESC differentiation. Considering that the key transcription factors such as EOMES and FOXA2, which are critical in endodermal development, are overexpressed following MeHg exposure, hESC specification could be affected with endodermal fate being overpromoted, abnormal organ development, and subsequently congenital diseases later in life. A summary of this effect is presented in the form of an adverse outcome pathway in Fig. 10. It is not possible to extrapolate the dose-response relationship observed in this in vitro study to humans. However, the treatment doses covered a range of total Hg typically found in human maternal blood (0.1-284.7 nM) ( Table 1). Our results suggest that the risk of MeHg during pregnancy, especially the newfound effects on endodermal development, needs more consideration for public health promotion and disease prevention.