Multi-omics approach reveals dysregulated genes during hESCs neuronal differentiation exposure to paracetamol

Summary Prenatal paracetamol exposure has been associated with neurodevelopmental outcomes in childhood. Pharmacoepigenetic studies show differences in cord blood DNA methylation between unexposed and paracetamol-exposed neonates, however, causality and impact of long-term prenatal paracetamol exposure on brain development remain unclear. Using a multi-omics approach, we investigated the effects of paracetamol on an in vitro model of early human neurodevelopment. We exposed human embryonic stem cells undergoing neuronal differentiation with paracetamol concentrations corresponding to maternal therapeutic doses. Single-cell RNA-seq and ATAC-seq integration identified paracetamol-induced chromatin opening changes linked to gene expression. Differentially methylated and/or expressed genes were involved in neurotransmission and cell fate determination trajectories. Some genes involved in neuronal injury and development-specific pathways, such as KCNE3, overlapped with differentially methylated genes previously identified in cord blood associated with prenatal paracetamol exposure. Our data suggest that paracetamol may play a causal role in impaired neurodevelopment.


INTRODUCTION
][7][8][9][10][11][12][13] The association is reported to be stronger with long-term exposure and higher dose. 14In 2019, the EUs pharmacovigilance safety committee (PRAC) reviewed all the available evidence, including nonclinical and epidemiological studies, regarding the impact of prenatal paracetamol exposure on impaired neurodevelopment in offspring.PRAC concluded that the available evidence is inconclusive, and recommended that the summary of product characteristics (SmPC) of paracetamol containing products should be updated to reflect the current state of scientific knowledge; ''Epidemiological studies on neurodevelopment in children exposed to paracetamol in utero show inconclusive results'' (PRAC, 2019). 15 published a consensus statement and literature review in Nature Reviews Endocrinology concluding that there is growing evidence supporting the hypothesis that in utero exposure to paracetamol can impair fetal development. 16This conclusion is still highly debated and contested by others (ENTIS, 2021) [17][18][19] reflecting a need for further research.
Using samples from the Norwegian Mother, Father and Child cohort (MoBa) biobank, we have previously shown differences in DNA methylation (DNAm) in cord blood from children diagnosed with ADHD who have been exposed to paracetamol (>20 days) during prenatal development compared to unexposed children. 20These findings suggest that DNAm might be involved in the pathogenesis of ADHD, but the causality and effect on neuronal differentiation and brain development is not known.It is well established that normal prenatal neurodevelopment involves cellular differentiation and establishment of cell-type specific epigenetic patterns, and that these events are prone to influences by environmental factors.2][23] If and how paracetamol modulates the apparent increased risk of NDDs is currently unknown.
Although paracetamol has been used for more than a hundred years, the mechanisms of action still remain unclear and appear to involve numerous physiological pathways, that vary among in vitro and in vivo studies. 24,25Following administration of therapeutic doses, paracetamol is primarily metabolized into pharmacologically inactive glucuronide and sulfate conjugates, while a small portion is oxidized to form the highly reactive metabolite, NAPQI. 26Although the precise neurotoxic mechanism of action for paracetamol and its metabolites remains unclear, it is known that deacetylation of paracetamol by N-deacetylase yields p-aminophenol, later conjugated with arachidonic acid by FAAH to form the active metabolite N-arachidonoylphenolamine (AM404). 27Paracetamol, once it crosses the blood-brain barrier, 28 may provide neuroprotection at low doses, whereas high doses may induce neurotoxicity, 29 possibly through oxidative stress.Furthermore, to our knowledge no studies have addressed whether and how the different cell populations derived from hESCs during early neuronal differentiation could potentially metabolize paracetamol and form bioactive compounds.
Recently, we established a protocol for neuronal differentiation of human embryonic stem cells (hESCs), which can be used in neuropharmacological studies. 30In the present study, we have used this model system of early human brain development and investigated the effects of paracetamol exposures on transcriptional and epigenetic regulation.2][33] By integrating multiple omics methods (bulk RNA-seq, bulk DNAm, single-cell RNA-seq and ATAC-seq, Table 1) we observed time and dose effects of paracetamol exposure during neuronal differentiation.

Neuronal differentiation exposure, timeline, and morphology
We investigated epigenetic and transcriptomic effects of exposure to paracetamol using an in vitro neuronal differentiation protocol that drives hESCs toward anterior neuroectoderm. 30,34The neuronal differentiation is divided into three stages: the neural induction phase (Stage I) ends at Day 7, the self-patterning phase (Stage II) ends at Day 13, and the FGF2/EGF2-induced maturation phase (Stage III) ends at Day 20 (Figure 1A).We replaced culture media daily, and the cells were exposed to 100 or 200 mM paracetamol during differentiation from Day 1 and onwards.2][33] Unexposed (control) and paracetamol-exposed cells were harvested for downstream analyses on Day 7 and 20.We also harvested control cells at the onset of differentiation (hESCs; Day 0), and on the intermediate timepoint that cells were passaged (Day 13) to assess whether paracetamol exposed cells had mRNA abundance changes related to proliferation or delayed differentiation.
The timeline of brightfield images of the control cells versus cells exposed to 100 (P100) or 200 mM paracetamol (P200) documents the morphological changes and cell culture density at differentiation Days 2, 4, 7, 8, 10, 13, 14, 17 and 20 (Figure 1B).Tightly packed neuroepithelial cells form the neural rosettes by Day 7 reassembled at the next stage under high density cell passaging and proceeded to maturation.We did not observe any distinct morphological changes in the differentiating cultures exposed to 100 (P100) or 200 mM paracetamol (P200).However, preliminary paracetamol titration experiments showed that exposing the differentiating cells to 400 mM paracetamol increased cell death and unpatterned morphology in cells and were thus discontinued.A set of representative images following the 20-day timeline in control cultures and 100, 200 and 400 mM paracetamol-exposed cells is presented in supplemental information (Figures S1A-S1C).

Validation of differentiation markers
The effect of paracetamol on gene expression at Days 7 and 20 was assayed using digital droplet PCR (ddPCR) (Figure 1C).Expression of the pluripotency transcription factors (TFs) POU5F1 and NANOG decreased significantly after neural induction.As anticipated, we observed an increase in expression of the neural markers SOX2, OTX2, FOXG1, and MAP2 on Day 7 and VIM, TUBB3, and NEUROD1 on Day 13.Notably, the expression of FOXG1, MAP2, and NES was significantly different between exposed and control cells on Day 7. Differential expression upon paracetamol exposure was also documented for filament NES, PAX6, and NEUROD1 on Day 20. (C) ddPCR results from 4 to 6 replicates of mRNA expression of selected marker genes from Days 0, 7, 13 and 20.Significant comparisons are marked with asterisks (Student's t test, *: p % 0.05, ***: p % 0.001).

Paracetamol-induced gene expression changes in genes involved in neural development
To delineate the effect on gene expression, we performed bulk gene expression analysis using RNA-seq in controls and paracetamol exposed cells (P100 and P200; Figures 2 and S2; Table S1).Principal component analysis showed that samples clustered according to differentiation day (Figure S2A; STAR Methods).Overall, when we compared gene expression from P100 versus control we identified 121 differentially expressed genes (DEGs, Figure 2A) and 1 433 DEGs between Day 7 and Day 20 P200 versus control (Figure 2B) (Table S1).Of these, 67 DEGs overlapped between P100 and P200 compared to control (Figure 2C).Pairwise comparisons between paracetamol (P100 and P200) and control at each day are shown in Figures S2B and S2C.The bulk RNA-seq analysis of the previously selected marker genes (Figure 1B) correlated well with the ddPCR results of selected marker genes (Figure S1D).
Gene set enrichment analyses (GSEA) identified downregulated biological processes (BPs) involved in synaptic organization, transmission, and regulation in the P100 time-response analysis (Figure 2D).Notably, one upregulated BP in the P200 time-response analysis were enriched for forebrain regionalization, whereas downregulated BPs were enriched for synaptic signaling, regulation of transsynaptic signaling, regulation of synaptic plasticity and regulation of neurotransmitter levels (Figure 2E).Moreover, the DEGs associated with P100 and P200 reflected BP enrichment of transmitter transport and regulation, synaptogenesis, synaptic organization, and plasticity.We also performed g:profiler 35 analysis of BPs on up and downregulated DEGs for P100 time-response and P200 time-response (Tables S1H and S1I).Top P100 time-response downregulated GO terms were generation of neurons, nervous system development, neuron differentiation and for P200 time-response downregulated GO terms were nervous system development, transsynaptic signaling and synaptic signaling (Tables S1J and S1K).Of the 46 BPs for downregulated DEGs P100 time-response, 40 overlapped with downregulated P200 time-response correlating with many overlapping GO terms identified above.There were no upregulated BPs for P100 time-response DEGs, and top identified BPs for upregulated P200 time-response GO terms were anatomical structure morphogenesis, animal organ morphogenesis and animal organ development (Table S1J and S1K).

Single-cell RNA sequencing reveals dose-specific changes in several major cellular processes after paracetamol exposure
To explore cell-type specific gene expression and maturation signatures over time, we performed single-cell RNA sequencing (scRNA-seq) analysis of control and paracetamol-exposed cells (P100 and P200) at Days 7 and 20 (Table S2).Specifically, our aim was to determine whether paracetamol exposure caused deviations in neuronal differentiation compared to control cells.A total of 15 201 cells (n = 6 924 Day 7 cells, n = 8 277 Day 20 cells) from two time-course experiments were aggregated and projected in Uniform Manifold Approximation and Projections (UMAPs, Figure 3; STAR Methods).The scRNA-seq data may also be visualized in the open access webtool (hescneuroparacet), where expression of genes can be explored per cell, cluster, and time point.
The cells clustered according to differentiation day and not exposure to paracetamol (Figure 3A).Consistent with previous results, 34 Seurat-predicted cell cycle phase showed an increase in G1 cells at Day 20 (Figure 3B).However, the existing cell cycle analysis tools are unable to decipher the proportion of neurons that are in G0 phase.Using CytoTRACE pseudotime differentiation trajectory analysis, the most differentiated cells were found at Day 20 (Figure 3C).The composite neuronal differentiation P1-P13 clusters were manually annotated with genes CRABP1, PAX6, TYMS, KIF20A, FGF17, HES5, WNT5A, ASCL1, GNG8, DLX1, SRSF9, VEPH1 and TAGLN2, respectively (Figure 3D).At Day 7, we observed a shift in cell composition from P1 to P2 in the cells exposed to paracetamol (Figure 3E).The effect was more prominent at Day 20, albeit not statistically significant (Tables S2I and S2J and STAR Methods), 53 with a higher proportion of paracetamol-exposed cells compared to controls annotated to cluster P5 and a more prominent effect at P200 exposure, whereas a lower proportion were annotated to P6, P9 and P10 for both concentrations (Figure 3E).

Paracetamol exposure is associated with differential expression of neural lineage markers
Performing GO analyses, we identified enrichment of BPs after paracetamol exposure in the scRNA-seq datasets.First, we identified the 10 most upregulated and downregulated BPs between control cells and P100 or control cells and P200 cells at Day 7 (Figures S3A and S3B, respectively).Notably, BPs involved in DNA replication and cell-cycle regulation were upregulated in cells exposed to both paracetamol doses, whereas we identified a downregulation of a BP which involves generation of neurons.DEGs in P100 cells compared to control cells were enriched for more neuronal specific annotations.Furthermore, P200 compared to control cells identified upregulated GO terms involved in cellular responses to toxic insults and DNA checkpoint activation.At Day 7, the relative gene expression of top 20 DEGs between P100 (Figure S3C) and P200 (Figure S3D) is shown.Furthermore, the gene expression of the top overlapping DEGs between the P100 and P200 cells compared to control cells at Day 7 (Figure S3E) identified gene specific similarities.We found changes in crucial genes, such as the master gene of forebrain development FOXG1 55 and genes of the HES and ID gene families involved in differentiation and neurogenesis. 56ext, we extracted the top 10 upregulated and downregulated BPs among the DEGs at Day 20 between control and P100 cells (Figure S3F) or control and P200 cells (Figure S3G).In both comparisons, downregulated GO terms included neuron/nervous system development and microtubule polymerization or depolymerization.Of the top 20 DEGs at Day 20 between paracetamol-exposed and control cells, we found major patterning TFs, such as NKX2-1 and EMX2 (Figures S3H and S3I).Moreover, genes that belong to the ZIC family, among other genes involved in Notch and Wnt signaling, were also identified (Figures S3H and S3I).The gene expression of the top overlapping genes between P100 and P200 cells were compared to control cells at Day 20 (Figure S3J).The analysis of GO terms delineated gene specific changes, such as upregulation of SELENOW (Figure S3D), previously associated with neuroprotection from oxidative stress. 57Paracetamol exposure induced differentiation lag as evidenced by the PAX6 expression in some P200 cells (Figure 1C).We also observed a downregulation of tissueand stage-specific genes, such as FABP7, ISL1, STMN2 and INA.In addition, TUBB1A, the isotype associated with postmitotic neurons, 58 and TUBB2B, that constitutes 30% of all brain beta tubulins, 59 were also found among the downregulated genes.Interestingly, PEG10 and C1orf61 and MIR9-1HG appear in the Day 20 P200 comparisons, genes which have recently been linked to cortical migration and intercellular communication, 60 further documenting how the P200 dose of paracetamol could affect proper network formation and cell-to-cell signaling.

Integration of scATAC-and scRNA-seq link paracetamol-induced changes in chromatin opening to transcriptional activity at Day 20
To understand whether paracetamol exposure during differentiation influenced the chromatin landscapes at Day 20 we performed single-cell assay for transposase-accessible chromatin using sequencing (scATAC-seq). 61We obtained 3 042 nuclei for Day 20 control and 3 480 and 4 282 nuclei for Day 20 exposed to 100 mM (P100) and 200 mM paracetamol (P200), respectively.First, we reanalysed scRNA-seq data from Day 20 controls, P100 and P200, and remapped the P3-P13 clusters (Figures 4A and 4B and webtool hescneurodiffparacet).The maturation trajectory cohered with the initial Day 7-Day 20 time point analysis (Figures 3C and 4C).Next, we mapped the scATAC-seq data (Figure 4D) to 15 scA-TAC-seq clusters (C1-C15; Figure S6A) that we integrated with the annotated scRNA-seq P3, P5, P9, P10 and P13 clusters (Figure 4E).The quality of the combined scATAC-seq datasets was documented with an even distribution of integrated clusters (P3, P5, P9 and P10) over TSS, promoters, exons, introns, and distal genomic regions (Figures S6B and S6C).The P13 cluster is represented by very few nuclei and displayed lower enrichment across all genomic regions.To further explore the differential chromatin accessibility across the genomes, we correlated distal accessible regions to gene expression.Cis-regulatory interactions with active genes were predicted by the integrative single-cell analysis using ArchR. 62This allows for identification of putative cis regulatory element (CREs) -gene pairs (peak-to-gene links (P2GLs) for biological and functional comparisons of Day 20 control cells with P100 and P200 exposed cells.We grouped P2GLs into five clusters and plotted heatmaps of gene expression and gene scores for all three datasets (P2GLs represent actively expressed genes linked to chromatin opening, Figures 4G-4I).Heatmaps show differential gene expression to the left, and gene scores for chromatin opening to the right, for CRE-gene pairs across cell clusters P3, P9 and P10 for control, and P3, P5, P9 and P10 for P100 and P200.The rows represent k-means 1-5 and are based on z-score-scaled associated gene expression levels.The scores for scRNA-seq and scATAC-seq were nicely correlating within the k-mean clusters showing good integration of the two modalities.When we integrated the scATAC-seq and the scRNA-seq datasets for Day 20, the clusters P3, P5, P9 and P10 and P13 were mapped (Figure 4E).Notably, the heatmaps of control and P100 and P200 cells were different with an absence or presence of the integrated cell cluster P5, also described as having differential prominence in the scRNA-seq data (Figure 3E).
We observed relatively similar numbers of putative CREs in controls (n = 30,771), P100 (n = 26,216) and P200 (n = 28,883) (Figures 4G-4I).Only 6,960 of the putative CREs overlapped between the three datasets, whereas 4,732, 3,435, and 6,534 putative CREs specifically overlapped between control and P100, control and P200, and P100 and P200, respectively (Figure 4J; Tables S3A-S3C).Moreover, we observed that many of the putative CREs were only detected in individual datasets (Day 20 control; n = 15 644, P100; n = 10 460 and P200; n = 12 304 CREs).These variations in putative CREs suggests that paracetamol exposure results in changes in chromatin accessibility.For example, in the locus of the neuronal marker gene STMN2, which showed higher gene expression in unexposed cells compared to exposed cells, showed changes in chromatin opening peaks in paracetamol exposed cells (Figures S3J, S5, and S6D).Furthermore, accessibility peaks in the ELAVL4 locus displayed differential chromatin opening in clusters P3, P5, P9, and P10.We also identified putative CREs in the exposed cells that correlated with higher ELAVL4 expression in control cells compared to P100 or P200 cells (Figure 4K; Tables S3A-S3C).

Paracetamol affects region-specific chromatin accessibility
The putative enhancer-gene interactions identified likely represent chromatin regulomes in the control and the paracetamol-exposed cells.We therefore compared the level of overlap of linked genes to better understand the effect of paracetamol exposure on chromatin regulation.Linked genes are defined as actively expressed genes by scRNA-seq that have linked chromatin accessibility regions in proximity to the gene locus.A larger proportion of linked genes overlapped between the Day 20 control and paracetamol exposed cells (n = 4,047), and the P100 (n = 1 224) and P200 (n = 1,048) cells (Figure 4L; Tables S3A-S3C and S3V).GO analyses of the linked genes revealed a common enrichment of BPs, such as nervous system development, biological regulation, neurogenesis, and neuron differentiation, varying in the different k-means (km) clusters (Figure S6D; Tables S3D-S3F).Linked genes identified in paracetamol exposed cells (P100, n = 652; P200, n = 915 and 1424 for both P100 and P200) indicate exposure-induced changes in gene expression.Interestingly, linked genes including neuronal lineage transcription factors (e.g., PAX6, NEUROD4, NEUROG1, SOX9, and SOX2) and genes involved in chromatin modification (e.g., histone H3K27 acetyltransferase EP300, 63 Histone H3K27me3 demethylase KDM6A 64 and chromatin remodeller SNF2 subfamily member SMARCAD1 65 ), suggest that paracetamol may contribute to modulation of transcriptional regulation and chromatin structure (Table S3V).In addition, GO analysis of the putative CREs identified enrichment of BPs per k-means group and sample are shown in Tables S3G-S3U.
Enriched TF motifs in the CREs may infer binding events that regulate gene expression programs.We identified TFs (n = 29) that were common for control, P100 and P200 cells, including NEUROD1, NEUROG1, SOX9, HMGA1 and members of the ONECUT and NHLH families, which have previously been described in the neuronal differentiation protocol. 34TFs were found to be common between P100 and P200 cells (n = 18) or dependent on paracetamol dose (Figures 4M, 4N, S3I, and S6F; Table S4).The TF footprinting predicts binding events due to the protection of the TF from Tn5 transposition in accessible chromatin. 62Based on the differential TF enrichment in P100 and P200 cells (Figure 4N), we generated TF footprints for OTX2, TBR1, and EMX2 as aggregate of genome-wide binding sequences adjusted for Tn 5 bias in the different integrated clusters (Figure 4O).Moreover, the UMAP plots for Motif Matrix show that these factors have binding events that are enriched in open chromatin in different integrated clusters.The OTX2 motifs were more enriched in clusters P3 and P10, whereas TBR1 motifs were mainly detected in cluster P9 and EMX2 motifs were found in clusters P9 and P10 (Figure 4O) suggesting that these TFs may potentially have cluster specific gene expression regulation.The computed sequence logos identified for these factors across enriched clusters are shown below each TF (Figure 4O).

Paracetamol induced changes in DNAm during neuronal differentiation
To assess whether exposure of hESCs undergoing neuronal differentiation to paracetamol induced DNAm changes, we analyzed control cells and cells exposed to 100 (P100) and 200 mM paracetamol (P200) at Day 7 and 20.The experimental set-up included analysis of control cells harvested at Day 0 and 13 as reference points of possible dysregulation (Figures 5 and S7; Table S1).Assessment of principal component analysis showed clustering of samples according to differentiation day (Figure S7C; STAR Methods).Overall, the distribution of average global DNAm was indistinguishable across exposed and control cells at all time points regardless of paracetamol dose (Figure S7A).In contrast, non-CpG DNAm levels decreased during differentiation and were lower in cells exposed to 200 mM paracetamol compared to control cells at Day 20 (Figure S7B).
To investigate whether the cells exposed to different doses of paracetamol (P100 and P200) respond differently than control cells from Day 7 to Day 20, we performed a DNAm time-response analysis.We observed no significant DNAm changes in P100 compared to control between Day 7 and Day 20 (not shown).In contrast, 3 113 CpGs responded differently to P200 compared to control cells over time (Figure 5A).CpGs showing an increase in DNAm are annotated to genes with key functions in dynamic cellular redox changes in the developing brain, such as the neural specification gene PRDM16 66 enriched for GO terms such as synaptic signaling and chemical synaptic transmission (Figure 5B).CpGs showing a decrease in DNAm are annotated to genes enriched for gene ontology (GO) terms involved in synaptic regulation, GABAergic signaling, and cell morphogenesis (Figure 5B).Overall, when we assessed DNAm levels at all CpGs (Figure 5C) compared to all significant CpGs (Figure 5D), we found a general dose-dependent increase in DNAm levels at significant sites in paracetamol-exposed cells compared to controls at both Day 7 and Day 20.The annotation of differentially methylated CpGs (DMCs) in relation to genes and CpG islands was similar across the different comparisons (Figures S7D and S7E).
Next, we performed comparisons of DNAm levels in paracetamol-exposed cells to controls at Day 7 (Figures S7F and S7H; Table S1) and Day 20 (Figures S3G and S3I; Table S1) for the different doses.As expected, we observed more significant DNAm changes after longer exposure (Day 20) and at higher concentration of paracetamol (P200) (Table S1).We found no dose-dependent DMCs at Day 7 whereas at Day 20 a larger number of CpGs (n = 8 940) were differentially methylated between P100 and P200.Moreover, there was some overlap between DMCs at 7 and Day 20 (Figure S7J).To define a regulatory role of paracetamol induced DNAm changes on gene expression, we assessed the overlap between DEGs and differentially methylated genes (DMGs) for the pairwise comparisons and time-response analyses (Table S1).The percentage of DMGs overlapping with DEGs varied between 3 and 63% for the different comparisons.The P200 time-response analysis revealed 180 overlapping genes, and some selected genes are visualized in Figure 6.DNAm levels for GRIK3, CACNA1D, ABAT, MAPT, and ANKRD6 were inversely correlated with gene expression, whereas DNAm levels for PAX7, CDH2 and WNT7B were positively correlated with gene expression at Day 20.GLI3 had DNAm levels that correlated with both negative and positive regulation (Figure 6).

Overlap of dysregulated genes in differentiating hESCs and cord blood from children exposed to paracetamol during pregnancy
We have previously identified an association between differential DNAm in cord blood and long-term paracetamol exposure during pregnancy in children with ADHD. 20To assess the translational potential and causality of these findings, we compared the dysregulated genes identified in the present model to the DMGs associated with paracetamol exposure in cord blood.In brief, the DMCs and the DEGs in paracetamol-exposed differentiating cells at Day 20 were correlated with DMCs/DMGs in cord blood. 20Interestingly, we identified 20 genes that were both differentially methylated and differentially expressed in Day 20 paracetamol-exposed cells which overlapped with DMGs identified in cord blood between paracetamol-exposed children with ADHD versus controls (Figure 7A).Furthermore, and only one gene (KCNE3) overlapped between Day 20 DMCs, Day 20 DMGs and paracetamol-exposed children with ADHD versus ADHD controls in cord blood (Figure 7B).
We assessed the expression of these 20 overlapping genes in scRNA-seq data (Figure 7C) and RNA-seq data (Figure 7D).Notably, genes involved in differentiation, such as GAB2 67 , or Notch and Hedgehog signaling pathways (for example, NOTCH4, 68 PTCH2, 69 SHISA2 70 ) show a dose-dependent downregulation of expression in paracetamol-exposed cells compared to controls.This high variation in expression levels after paracetamol exposure compared to controls, indicates a dose effect.Among genes identified in Day 7 P200 cells, we observed upregulation of ZDHHC14 71 and SGK1, 72 that are associated with control of neuronal excitability and neuronal response to injury, respectively.Several toxic insult response genes were also upregulated at Day 20 in the P200 treated cells, such as IER3, 73 SPRED2, 74 GPR130, 75 and SLC22A23. 76

DISCUSSION
Stem cell based systems for toxicological studies of medications allow for more realistic testing of concentrations in vitro and are often monitored using multi-omics approaches. 77Due to the inaccessibility of early human neurodevelopment for studies of the effects of paracetamol we have utilized a well characterized neuronally differentiating hESC protocol 30 as an easily accessible and 3R's compliant source.To our knowledge, this is the first multi-omics study investigating how paracetamol exposure affects epigenetic and transcriptional programmes using an in vitro hESC model of neuronal differentiation. 34e identified differential expression of genes enriched for BPs involved in transmitter transport and regulation, synaptic organization and synaptogenesis and synaptic plasticity.Interestingly, in cells exposed to the high dose of paracetamol (P200), we also identified enrichment of GO categories reflecting possible patterning deviations from forebrain differentiation.Thus, we observe an effect of paracetamol exposure on transcriptional dysregulation and possible developmental delays during neuronal differentiation.Additionally, among the common DEGs, irrespective of paracetamol dose, we identified SORCS3, which encodes a brain-expressed transmembrane receptor associated with neuronal development and plasticity that has been previously identified in a GWAS meta-analysis of ADHD significant risk loci. 78The scRNA-seq analysis identified dose-dependent DEGs involved in several major processes during neuronal differentiation.Exposure to paracetamol resulted in a subtle shift in cell populations from P1 to P2 at Day 7. At Day 20 the P6, P9 and P10 clusters decreased whereas P5 markedly increased in P200.The DEGs link paracetamol exposure to genes involved in cell-cycle length and phase transition and genes important for neuronal maturation, neurite outgrowth, cortical neurogenesis, expression of neurotransmitter transporters and WNT and FGF signaling.Further, the data also documented dose-dependent differential expression of crucial spatiotemporally regulated TFs associated with neuronal differentiation.These findings provide evidence of transcriptomic dose-dependent effects of paracetamol exposure related to cellular response to toxic insults and fate-determination deviation queues at Day 20.[81]  To further characterize and investigate epigenetic mechanisms involved in altered cellular heterogeneity and differentiation, we performed scATAC-seq to map CREs and compare the chromatin landscape changes between control and cells exposed to P100 and P200 at Day 20 (hescneurodiffparacet).Interestingly, the results from these analyses identified differences in several chromatin accessible regions in genes previously associated with fate mapping and neurogenesis, which have been identified in risk determination studies focusing on cognition and autism spectrum disorders (e.g., DLK1, DIO3, IPW 82,83 ) and ADHD (e.g., MEG3 84 ).Genome-wide association studies (GWASs) suggest that cognitive disorders might also result from the cumulative effect of gene variant regulation and parent-of-origin effect, 84 which in ADHD has been associated with genes such as DDC, MAOA, PDLIM1, and TTR.We also identified chromatin openness differences related to genes identified in the paracetamol exposed cells.For example, PCDH7, which acts as a cue for axonal guidance via roles in cell adhesion, was among the genes that were both differentially expressed and had different chromatin openness in Day 20 P200 cells.Dysregulation of PCDH7 could be relevant to the semaphorin-plexin signaling documented by the enriched BPs in the P100 cells, and it is noteworthy as it has also been associated to ADHD. 78e also identified dose-dependent differences in chromatin accessibility of the neural specification gene locus PRDM16, which correlates well with DNAm analyses where exposure to paracetamol was associated with increased DNAm at CpGs annotated to PRDM16.This gene has been linked to regulation of transcriptional enhancers activating genes involved in intermediate progenitor cell production and repressing genes involved in cell migration. 66The molecular mechanism of PRDM16 remains unknown but has been associated with reactive oxygen species regulation in the embryonic cortex. 85ntegrative analyses of scRNA-seq and scATAC-seq to explore the links between chromatin accessibility and gene expression documented high correlation and showed an overlap between scATAC-seq and scRNA-seq cell clusters.This analysis identified many putative CREs, and linked genes affected by paracetamol exposure.Although a portion of linked genes overlapped between the control and exposed cells, there was only a small overlap between the putative CREs.This suggest that many local changes in chromatin opening were affected by paracetamol exposure and that these CREs may regulate different gene expression programmes in control and paracetamol exposed cells.How paracetamol regulates chromatin opening is not known, and we have limited knowledge on the effect of paracetamol on the TFs recognition of DNA motifs.However, our data suggest that genes linked to putative CREs that were affected by paracetamol exposure have a role in regulation of transcription and chromatin.In particular, a linked gene specific for paracetamol exposure was EP300, which codes for a histone acetyltransferase responsible for the active enhancer mark H3K27ac. 86EP300 expression decreased subtly upon paracetamol exposure both at Day 7 and 20, and these changes in EP300 levels may change the number of active enhancers in the paracetamol exposed cells.Furthermore, we identified TFs such as EMX2, TBR1 and OTX2 whose regulatory activity is affected by paracetamol exposure at the differentiation endpoint.This is suggestive of a mechanism for how paracetamol exposure may change gene expression programmes in early brain development.More work is needed to follow up the functional role of paracetamol exposure on the enhancer regulation and gene expression.
We found a general dose-dependent increase in DNAm in paracetamol exposed cells compared to controls at both Day 7 and Day 20.The genes were enriched for BPs involved in cell morphogenesis and adhesion, but also synaptic transmission and signaling with a focus on catecholamines and GABAergic transmission.8][89][90] It is biologically plausible that medications that cross the placenta and the blood-brain barrier, may interfere with normal fetal brain neurodevelopment, previously shown for valproic acid, and more recently for topiramate. 91This is especially relevant for substances that cross the blood-brain barrier, which is considered functional by the 8th week of gestation. 92However, invasive testing of the fetal compartment in human is rarely indicated, which makes data on safety during pregnancy for most substances scant.][95] We previously found DMGs in cord blood from children with ADHD exposed to long-term maternal paracetamol use. 20However, the causality, and relevance of such findings to brain development are not known.The hESCs used in our experiments represent the inner cell mass (ICM) of Day 6 preimplantation-blastocysts. 96The ICM is known to differentiate to form the primitive ectoderm around Day 7.5 post fertilization, 97 but placental circulation is not established before Day 17-20. 98Prior to this, the fraction of paracetamol that reaches the developing embryo does so via passive diffusion. 99Interestingly, irrespective of the temporal differences of the two models, 100 comparing the differentially expressed and methylated genes identified in our study revealed overlap of several genes identified in cord blood with potential compelling relevance to early brain development.2][103] To our knowledge this is the first study that has identified an effect of paracetamol on differential DNAm of KCNE3 in a neuronal differentiation model of hESCs.Differential DNAm at KCNE3 was also identified in our MoBa study in cord blood. 20Furthermore, KCNE3 was also identified in a DNAm analysis of extremely low gestational age new born (ELGAN) cohort. 93CNE3 is an interesting candidate, which encodes a voltage-gated ion channel with important functions in regulating release of neurotransmitters and neuronal excitability. 104,105Among the overlapping genes, there was also IER3, whose differential DNA methylation at specific CpG sites has been previously associated with in utero exposure to bisphenol A (BPA). 106The same study had identified TSPAN15, which is also significantly differentially methylated between control and P200 cells in our neuronal differentiation model.
We also compared our results to previously published studies by exploring expression of relevant genes in the webtool (hescneuroparacet).In a rat model of fetal paracetamol exposure, 107 dysregulation in specific ABC efflux transporters and related enzymes was found after chronic treatment of the mothers in E19 fetal brain and choroid plexus.In our model of paracetamol exposure, we also identified differential expression of the relevant genes ABCA1, ABCC5, ABCD3, and GSTM3.Moreover, the expression of ABCC4, which encodes the multidrug resistance-associated protein MPRP4, known to play a role in paracetamol efflux was also upregulated at Day 20. 108 In another model assessing chronic paracetamol exposure effects in E15-19 rat brains, 109 similar to our findings, an increased expression of genes related to proliferation (e.g., MKI67, MTDH), cytoskeletal structural genes (e.g., NEFL at Day 20), metabolic stress alleviation (e.g., PFKP at Day 7) or decrease in genes related to differentiation (SOX12 at Day 7) migration and dendrite orientation (such as CRMP1 at Day 20) was found.Notably, a significant increase in APOE expression at P200 cells at Day 7 was observed indicating stress or neuronal damage, 110 and a similar response by the increased expression for glutamate transporter EAAT4. 111A DNA repair gene, RAD54L, was upregulated in P200 treated cells at Day 7 (Table S2C), and interestingly this gene was also found differentially expressed in the prefrontal cortex of paracetamol-exposed mice. 112Also, a stress response transcription factor ATF4 was found significantly downregulated in P200 cells at Day 7 (Table S2C) and has been shown to be differentially expressed in a single-cell study of mice hepatocytes upon paracetamol overdose. 113Moreover, the significant IL6ST (GP130) and MEF2A 114 upregulation and the similarities of our findings to the other in vivo models, validate both the trajectories identified by our approach, and the chosen paracetamol neurotoxicity platform for early development studies.
Abiding by the Findability, Accessibility, Interoperability, and Reusability (FAIR) principles, we provide full access to the scRNA-seq and scATAC-seq datasets in open access shiny web platforms for scRNA-seq (hescneuroparacet) and integrative scATAC-seq/scRNA-seq (hescneurodiffparacet). 115,116 These webtools allow for data correlation with other published gene expression datasets, and enable plotting, exporting, and downloading high resolution figures of gene expression and gene co-expression analysis UMAPs of any gene.Furthermore, there are multiple options for dataset exploration and visualization, such as heatmaps, violin-, box-, proportion-and bubble plots in different tabs, where gene expression may be viewed per cell, cluster, treatment and timepoint.Chromatin opening can be explored for genomic loci of interest for input sample, ATAC clusters or integrated clusters in a genome browser view with peak-to-gene link annotations (instantly plotted for 500 random cells each time).In addition, the processed datasets and code are made available for customization and integration with other studies.To our knowledge this is the first in vitro study of paracetamol exposure of hESCs under neuronal differentiation, that provides access to scRNA-seq and scATAC-seq data via interactive webtools.

Limitations of the study
Our study has several limitations and strengths.First, paracetamol is metabolized into inactive glucuronide and sulfate conjugates and the highly reactive metabolite, NAPQI, in the liver, but the presence of the required enzymes, transporters and breakdown molecules, were not measured in this study.Future studies are needed to characterize paracetamol metabolism in hESCs and potential effect of paracetamol conjugates and metabolites on neuronal differentiation.Moreover, the protein expression changes, and enzyme activities were not explored, as this was beyond the scope of this study.This study aimed to delineate the epigenetic impact of paracetamol on early human cortical developmental events.We used both bulk-and scRNA-seq and comparing these two datasets identified a varying degree of overlapping DEGs from 1 to 27%, confirming that both datasets are complementary to detect transcriptomic changes induced by paracetamol.Using a multiomics approach, we unveiled the diversified transcriptional networks related to paracetamol exposure.We could not accurately classify the maturation properties of the cells, or the proportion of cells in G0 phase, highlighting the necessity of new tools to deconvolute the neuronal G0 phase.The results of scATAC-seq analysis presented here are inherently limited by the relatively low genome-per-cell coverage.That means that some open chromatin regions that could have proved relevant for the individual cells or cell populations, may have been missed.Finally, the in vitro results presented here need to be validated by in vivo models and targeted human dataset exploration in future studies.This could confirm whether the pathways identified can explain the paracetamol-induced adverse effects in the early human brain.
In conclusion, using an in vitro hESC model of neuronal development, we identified altered gene expression, DNAm and chromatin openness involved in key neuronal differentiation processes at bulk and single-cell RNA levels.An overlap of genes identified in this model and in the cord blood of neonates exposed long-term to paracetamol during pregnancy, points toward altered epigenetic regulation of early brain development.Identification of common DNAm modification sites and chromatin openness regions with a regulatory role on gene expression can identify loci and underlying mechanisms involved in neurotoxicity of drugs on the developing fetus with potential effect on long-term outcomes.Such findings could strengthen causal inference and clinical translation of altered DNAm, such as in cord blood, on brain development.However, these in vitro results need further validation for potential clinical translation.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: Neuronal differentiation of hESCs and exposure to paracetamol hESCs HS360 were differentiated, as described, 30 in two separate time-course experiments, to investigate the in vitro effects of paracetamol exposure in physiologically relevant concentrations to long-term exposure in vivo.2][33] Thus, cells undergoing neuronal differentiation were exposed to changes of medium only (control) or medium containing 100 or 200 mM paracetamol from Day 1 to Day 20.Control cells were harvested at day 0, 7, 13 and 20, while cells exposed to 100 or 200 mM paracetamol were harvested at Day 7 and 20 (Table 1).

METHOD DETAILS DNA/RNA isolation
Genomic DNA and total RNA were isolated by direct lysis in the culture vessels followed by column-based isolation using RNA/DNA purification kit (Norgen Biotek).RNase-Free DNase I Kit (Norgen Biotek) was applied for on-column removal of genomic DNA contamination from RNA isolates.Five RNA isolates were processed using the RNeasy Mini Kit (Qiagen) followed by DNase-treatment using the RNAse-Free DNase Set (Qiagen).All isolations were performed according to the manufacturer's instructions.Nucleic acid quantification was performed using Qubit (ThermoFisher Scientific), purity was measured using Nanodrop 2000 (ThermoFisher Scientific), while RNA and DNA integrity was assessed using 2100 Bioanalyzer (Agilent Technologies) and 4200 TapeStation (Agilent Technologies), respectively.

ddPCR and RNA expression analysis
Reverse transcription of total RNA was performed using QuantiTect Reverse Transcription Kit (Qiagen).Subsequent ddPCR reactions were set up using ddPCR Supermix for Probes (No dUTP) (BioRad) and Taqman assays (ThermoFisher) or Universal Probes (Roche) in combination with target primers (Eurofins) as outlined in key resources table.Droplets for ddPCR amplification were generated using the QX200 Droplet Generator (BioRad).Data acquisition and primary analysis was done using the QX200 Droplet Reader (BioRad) and QuantaSoft software (BioRad).All steps were performed according to the manufacturer's instructions.To calculate the number of target copies per ng RNA input, samples were normalized using RPL30 and RAF1 as normalization genes. 159Statistical comparisons were performed in R using t-test in ggpubr package v.0.4.0. 151Results were visualized in R using the tidyverse package. 150lk RNA-seq The sequencing library was prepared with TruSeq Stranded mRNA Library Prep (Illumina) according to manufacturer's instructions.The libraries (n=39) were pooled at equimolar concentrations and sequenced on an Illumina NovaSeq 6000 S1 flow cell (Illumina) with 100 bp paired end reads.The quality of sequencing reads was assessed using BBMap v.34.56 160 and adapter sequences and low-quality reads were removed.The sequencing reads were then mapped to the GRCh38.p5index (release 83) using HISAT2 v.2.1.0. 161Mapped paired end reads were counted to protein coding genes using featureCounts v.1.4.6-p1. 162Differential expression analysis was conducted in R version 3.5.1 163 using SARTools v.1.6.8 152 and the DESeq2 v.1.22.1. 128For the time-response analysis, the time effect of 100 or 200 mM paracetamol compared to control (no paracetamol) from Day 7 to Day 20 (Day 7 as reference), was modelled using the following design: '' Concentration + Day + Concentration:Day'' using the interaction term to find any treatment-specific differences over time.Genes were considered significantly differentially expressed with an FDR < 0.05.Normalized counts were visualized using the tidyverse package v.1.3. 150The heatmaps were generated using the pheatmap package version 1.0.12. 164The Wald-test was used to calculate p-values and Benjamini-Hochberg was used to correct for multiple testing.The Gene Set Enrichment Analysis (GSEA) analysis was performed using pre-ranked gene lists, based on p-values and direction of expression change by a competitive test, into GSEA software v.4.1.0 129to identify enrichment of BP terms.The size of the analysed gene sets was restricted to 15-1000 genes, and the chip annotation used was ''Human_ENSEMBL_Gene_ID_ MSigDB.v7.2.chip''.Moreover, we performed analysis of BPs on differentially expressed genes from cells exposed to 100 or 200 mM paracetamol compared to control cells over time (Day 7 -Day 20) (Tables S1H and S1I) analysed using a server based g:Profiler 53 showing results with correlation greater than 0.45_significant.

Bulk DNAm analysis
DNAm status of 43 samples were assessed using the Infinium MethylationEPIC BeadChip v.1.0_B3(Illumina).Quality control and pre-processing of the raw data was performed in R using Minfi v.1.36.0. 130No samples were removed due to poor quality (detection p values >0.05).Background correction was performed using NOOB method 165 and b values (ratio of methylated signal divided by the sum of the methylated and unmethylated signal) were normalized using functional normalization. 166Probes with unreliable measurements (detection p values >0.01) (n = 12,538) and cross-reactive probes (45) (n = 42,844) were then removed, resulting in a final dataset consisting of 810,477 probes and 43 samples.Probes were annotated with Illumina Human Methylation EPIC annotation 1.0 B5 (hg38).Differential DNAm analysis for pairwise and time-response comparisons was performed on the M values (log2 of the b values) using the limma

Figure 1 .
Figure 1.Neuronal differentiation of hESC to model neurodevelopmental effects of paracetamol (A) Schematic illustration of the different phases of the neuronal differentiation protocol from Day 0 to Day 20.Cells were exposed to 100 or 200 mM paracetamol from Day 1 and onwards.The effect of paracetamol on gene expression and epigenetic profiles was evaluated at Day 7 and Day 20.

Figure 2 .
Figure 2. Paracetamol exposure of differentiating hESCs modulates expression of genes involved in neural development (A and B) Volcano plots of longitudinal differences in gene expression difference between Day 7 and Day 20 (Day 20 minus Day 7) in A) P100 and B) P200 cells compared to control cells.Marked in red are significantly differentially expressed genes with an FDR <0.05.(C) Corresponding Venn diagram shows the number of overlapping DEGs between the P100 and P200 time-response (TR) comparisons (Day 20 minus Day 7 compared to control).(D and E) Top 20 enriched BPs based on GSEA in (D) P100 cells and (E) P200 cells compared to control cells over time (Day 20 minus Day 7).Normalized enrichment scores are represented as bars with FDR q-values.(F) Gene expression levels of selected DEGs in both P100 and P200 cells compared to control cells over time (TR) (overlapping area in C).

Figure 3 .
Figure 3. scRNA-seq analysis revealed shifts in cell type composition of differentiating cells exposed to paracetamol (A-E) Day 7 and Day 20 control cells and cells exposed to 100 mM (P100) or 200 mM (P200) paracetamol were visualized with UMAP and colored by (A) sample identity, (B) Seurat predicted cell cycle phase, (C) CytoTRACE pseudotime differentiation trajectory, (D) defined Seurat clusters at resolution 0.4 with corresponding gene annotations and (E) cell proportions per cluster.(F) UMAP plots of cells colored by SingleR cell annotation to Early Human Brain reference with corresponding cell annotations.Cell types starting with ''e'' are hESC derived cells and cell types starting with ''h'' are in vivo human embryo cell types.Nb1-4; neuroblasts, Prog1-2; neuronal progenitors, Rglf; radial

Figure 4 .
Figure 4. Effects of paracetamol on chromatin accessibility and integration with scRNA-seq at Day 20 (A-C) scRNA-seq UMAP plots colored by A) sample identity, B) scRNA-seq clusters or C) CytoTRACE pseudotime differentiation trajectory.(D and E) scATAC-seq UMAP plots of colored by D) sample identity and E) remapped clusters following constrained alignment of cell populations by scATAC-seq and scRNA-seq integration.(F) Supervised pseudotime trajectory of integrated clusters.(G-I) Heatmaps of scATAC-seq and scRNA-seq side-by-side representing peak to gene links (P2GLs) in G) control cells, H) P100 cells and I) P200 cells.Rows were clustered using k-means (k = 5).(J) Overlap of putative CREs from P2GLs analyses of Day 20 control versus P100 and P200.(K) ELAVL4 locus browser view (GRCh38.p13)from 5000 cells with ATAC-seq signals in integrated clusters.Differential chromatin opening (black triangles), ATAC peaks (red), P2GLs (arcs) and putative CREs enriched in P200 or P100 and P200 cells (colored triangles) are shown.The bar plot to the right shows the corresponding percentage of cells expressing ELAVL4, and the plot below represent the expression levels.(L) A Venn diagram representing overlap of linked genes.(M) Venn diagram of TF regulators identified from P2GLs integrative analysis.(N) A selection of positive TF regulators computed using gene integration scores with motifs in the putative CREs.(O) Footprints for selected TF regulators OTX2, TBR1 and EMX2 demonstrating preferential opening per cluster; below are the corresponding ArchR motif deviation scores, motif matrix and representative sequence logos identified in the scATAC-seq dataset.

Figure 5 .
Figure 5. Exposure to paracetamol induces changes in DNAm over time during neuronal differentiation (A) Volcano plot showing the effect of 200 mM paracetamol from Day 7 to Day 20.CpGs with adjusted p value <0.05 were considered significant.(B) Corresponding enriched GO-terms for CpGs showing a decrease in DNAm (top) or an increase in DNAm (bottom) over time (Day 20 -Day 7) in P200 cells compared to control cells.Adjusted p-values are indicated by color.(C) Average DNAm levels per sample for all CpGs across differentiation for Day 0, 7, 13, and 20.(D) Average DNAm levels per sample for all significant CpGs (paracetamol-exposed cells vs. control comparisons) at Day 0, 7, 13, and 20.(C and D) *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Figure 6 .
Figure 6.DNAm and gene expression for selected overlapping DMCs and DEGs (A-I) Change in DNAm and gene expression levels between control and P100 or P200 for (A) GRIK3, (B) PAX7, (C) CDH2, (D) CACNA1D, (E) WNT7B, (F) GLI3, (G) ABAT, (H) MAPT and (I) ANKRD6.Each point represents a matched replicate between RNA-seq and DNAm and black lines represents a linear regression line of the mean of values.

Figure 7 .
Figure 7. Overlap with differentially methylated genes in cord blood (A and B) Venn diagrams showing the overlap between Day 20 DEGs and DMCs for paracetamol comparisons and (A) paracetamol-exposed children with ADHD versus controls and (B) paracetamol-exposed children with ADHD versus ADHD controls in cord blood from MoBa.(C and D) Gene expression levels of overlapping genes identified in cord blood derived from (C) RNA-seq and (D) scRNA-seq data.

Table 1 .
Overview over presented datasets

TABLE
d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS B hESC culture and maintenance B Neuronal differentiation of hESCs and exposure to paracetamol d METHOD DETAILS B DNA/RNA isolation B ddPCR and RNA expression analysis B Bulk RNA-seq B Bulk DNAm analysis B Collection of cells and scRNA-seq B scRNA-seq data analysis B scATAC-seq library preparation and sequencing B scATAC sequencing analysis d QUANTIFICATION AND STATISTICAL ANALYSIS SUPPLEMENTAL INFORMATION Supplemental information can be found online at https://doi.org/10.1016/j.isci.2023.107755.(Continued on next page)