Neuronal differentiation pathways and compound-induced developmental neurotoxicity in the human neural progenitor cell test (hNPT) revealed by RNA-seq.

There is an increased awareness that the use of animals for compound-induced developmental neurotoxicity (DNT) testing has limitations. Animal-free innovations, especially the ones based on human stem cell-based models are pivotal in studying DNT since they can mimic processes relevant to human brain development. Here we present the human neural progenitor test (hNPT), a 10-day protocol in which neural progenitor cells differentiate into a neuron-astrocyte co-culture. The study aimed to characterise differentiation over time and to find neurodevelopmental processes sensitive to compound exposure using transcriptomics. 3992 genes regulated in unexposed control cultures (p ≤ 0.001, log2FC ≥ 1) showed Gene Ontology (GO-) term enrichment for neuronal and glial differentiation, neurite extension, synaptogenesis, and synaptic transmission. Exposure to known or suspected DNT compounds (acrylamide, chlorpyrifos, fluoxetine, methyl mercury, or valproic acid) at concentrations resulting in 95% cell viability each regulated unique combinations of GO-terms relating to neural progenitor proliferation, neuronal and glial differentiation, axon development, synaptogenesis, synaptic transmission, and apoptosis. Investigation of the GO-terms 'neuron apoptotic process' and 'axon development' revealed common genes that were responsive across compounds, and might be used as biomarkers for DNT. The GO-term 'synaptic signalling', on the contrary, whilst also responsive to all compounds tested, showed little overlap in gene expression regulation patterns between the conditions. This GO-term may articulate compound-specific effects that may be relevant for revealing differences in mechanism of toxicity. Given its focus on neural progenitor cell to mature multilineage neuronal cell maturation and its detailed molecular readout based on gene expression analysis, hNPT might have added value as a tool for neurodevelopmental toxicity testing in vitro. Further assessment of DNT-specific biomarkers that represent these processes needs further studies.

acrylamide

adverse outcome pathway

complement DNA

chlorpyrifos

developmental neurotoxicity

4′,6-diamidino-2-phenylindole

differentially expressed gene

dulbecco's phosphate-buffered saline

gene ontology term

fluoxetine

human embryonic stem cell

human neural progenitor test

integrated approaches to testing and assessment

inhibitory concentration, x = percentage decrease

log2 fold changes

methyl mercury

neural progenitor cell

organisation for economic co-operation and development

principal component analysis

quantitative polymerase chain reaction

RNA sequencing

valproic acid

Introduction

The prevalence of neurologic diseases with a developmental origin, like autism spectrum disorder, bipolar disorder, attention deficit hyperactivity disorder, and schizophrenia seems to be rising (EPA, 2019; Grandjean and Landrigan, 2014, 2006; Hertz-Picciotto and Delwiche, 2009; Schettler, 2001). The aetiology of these disorders is multifactorial, even for those with a strong genetic background (De Felice et al., 2015; Hallmayer et al., 2011). It is suggested that compound exposures may play an important role in these disorders. In recent decades, for example, a sharp upward trend has been noticed in autism cases that could not be fully explained by other environmental factors (EPA, 2019). Information on the neurotoxic potential of many compounds is scarce, but epidemiological data and animal experiments suggest that there is indeed a link between compound exposure and adverse outcomes (Carlson et al., 2020; Miodovnik, 2019; Zeliger, 2008). Strong causative links between compound exposures and specific human health outcomes, however, are lacking (De Felice et al., 2015; Grandjean and Landrigan, 2006).

While proven useful in the past, animal studies to test the safety of chemicals or pharmaceuticals may not always be sufficient in representing human brain development (Paparella et al., 2020). Neurodevelopment becomes increasingly complex and species-specific over the course of development (Cardoso-Moreira et al., 2019; Rice and Barone, 2000; Silbereis et al., 2016), and failure to account for human-specific effects is reflected in the high attrition rates when drugs enter clinical trials (Authier et al., 2016; Weaver and Valentin, 2019). In the context of both chemical and pharmaceutical safety assessment, DNT studies are carried out in specific cases only, and the available test guidelines (OECD TG 426 or 443 DNT cohort) in place only assess gross morphological, physiological and behavioural alterations (Beekhuijzen, 2017; Fritsche et al., 2018b; Tsuji and Crofton, 2012). Furthermore, animal-testing is low-throughput, costly, and ethically debatable (Crofton et al., 2012; Meigs et al., 2018; Paparella et al., 2020; Tsuji and Crofton, 2012). A different, human-relevant, and mechanism-based approach is therefore highly needed (Bal-Price et al., 2018b; Crofton et al., 2014; Fritsche et al., 2018a).

Animal-free safety testing of compounds requires a rethinking of the current testing paradigm (National Research Council, 2007). The intact organism may in principle be replaced by a series of in vitro assays that together represent the key processes in human neurodevelopment in combination with in silico models to integrate experimental outcomes and extrapolate them to human hazard assessment (Bal-Price et al., 2018a; Hessel et al., 2018). Current initiatives focus on building frameworks in which these key neurodevelopmental processes are defined, for example in ontologies, adverse outcome pathways (AOPs), and integrated approaches to testing and assessment (IATA) (Bal-Price et al., 2018b; Fritsche et al., 2017; Hessel et al., 2018; Masjosthusmann et al., 2020; Sachana et al., 2019, 2021b). An ontology is defined as a framework to qualitatively and quantitatively integrate and structure relevant data from various sources (Desprez et al., 2019). Interpretation of the data and data integration of in vitro data in these frameworks are major challenges that need to be overcome, but some first steps are taken towards addressing these issues in non-DNT and DNT contexts (Dobreniecki et al., 2022; Hernández-Jerez et al., 2021; Masjosthusmann et al., 2020; Spînu et al., 2022).

The use of human stem cell-based models in the abovementioned frameworks is especially promising, because they naturally transition through various developmental stages upon differentiation (E Fritsche et al., 2018). However, due to the inherent reductionist nature of these and other in vitro assays, their technical applicability domain (i.e. which kind of compounds are feasible to test in cell-based systems) as well as their biological domain (i.e. which biological processes does the assay mimic) are limited and need to be defined (Carlson et al., 2020; Hessel et al., 2018).

A first step to define this biological domain can be achieved using a transcriptomics approach. Whole transcriptome analysis allows for an unbiased assessment of the expression of genes that are involved in (neurodevelopmental) biological processes, thereby identifying which of these processes may be represented in a given model. These neurodevelopmental processes may include neural progenitor commitment, neuronal and glial cell differentiation and maturation, synaptogenesis, and network formation (Fritsche et al., 2017; Hessel et al., 2018; Wegner et al., 2020). Moreover, transcriptomics can provide mechanistic insights into the mode-of-action of compounds and offer biomarkers for DNT in vitro.

This study builds on previous research in which we have shown successful differentiation of hESC-derived neural progenitor cells (NPCs) into a neuron-astrocyte co-culture. Importantly, spontaneous electrical activity was already present after three days of differentiation (de Leeuw et al., 2020c). The current study aimed to develop an exposure protocol using this differentiation method that is suitable for DNT testing: the human neural progenitor test (hNPT). To this end, the biological domain of hNPT and potential endpoints to measure DNT in this system were explored over time and in response to five commonly used compounds in in vitro assays for DNT with diverse mechanisms of action using a whole transcriptome approach with RNA sequencing (RNA-seq).

Materials and methods

Neuron-astrocyte differentiation from embryonic stem cell-derived neural progenitor cells: the human neural progenitor test (hNPT)

Neural progenitor cells (NPCs) were generated from H9 human embryonic stem cells (WA09, passage 26, WiCell, Madison, WI, USA) as described previously (de Leeuw et al., 2020c). Cells are regularly screened for mycoplasma infections every year. For each experiment, NPCs were thawed and maintained according to the neuronal induction protocol from Stemcell (Document #28782, Stemcell Technologies, Vancouver, Canada). Briefly, cells were kept in a humidified chamber (37 °C, 5% CO2, 3% O2) in STEMdiff™ Neural Progenitor (NP) medium (Stemcell) on Poly-l-Ornithine (PLO, 15 μg/mL, Sigma-Aldrich, Saint Louis, MO, USA) - laminin (10 μg/mL, Sigma-Aldrich) coated 6-well plates (Corning, New York, NY, USA). Culture medium was refreshed daily for 7 days until the NPCs were 100% confluent. For neuron-astrocyte differentiation, NPCs were dissociated and seeded at 2.56 × 105 cells/cm2 on PLO-laminin-coated 8-well micro-slides (Ibidi, Gräfelfing, Germany), 12/24/48-well plates (Corning) or 96-well plates (Greiner Bio-One, Kremsmünster, Austria), depending on the application. Cells were kept for one day in Neural Progenitor medium (NP medium), after which they were switched to differentiation medium adapted from Gunhanlar et al. (2018). This medium consisted of neurobasal medium, 10 μl/mL N2 supplement (100x), 20 μl/mL B-27 without retinoic acid supplement (50x), 10 μl/mL 5000 IU/mL Penicillin/5000 μg/mL Streptomycin, 10 μl/mL nonessential amino acids (100x), 20 ng/mL recombinant glial cell line-derived neurotrophic factor (GDNF) and 20 ng/mL recombinant brain-derived neurotrophic factor (BDNF; all from Gibco, Waltham, MA, USA), 1 μM dibutyryl cyclic adenosine monophosphate (Sigma-Aldrich), 200 μM ascorbic acid (Sigma-Aldrich) and 2 μg/mL laminin (Sigma-Aldrich). Recombinant ciliary neurotrophic factor (CNTF; Gibco) was added during the first week of neuron-astrocyte differentiation at a final concentration of 10 ng/mL. GDNF, BDNF, ascorbic acid, dibutyryl cyclic adenosine monophosphate and CNTF were all dissolved in sterile solutions. Half of the medium was refreshed every three to four days for up to 28 days.

Compound exposure and analysis

Compounds were all purchased from Sigma-Aldrich freshly prepared for each exposure. Acrylamide (ACR, CAS# 79-06-1, ≥99% purity) and valproic acid sodium (VPA, CAS# 1069-66-5, ≥98% purity) were dissolved in medium, and chlorpyrifos (CPF, CAS# 2921-88-2, ≥98% purity), fluoxetine hydrochloride (FLX, CAS# 59333-67-4, ≥98% purity) and methyl mercury (II) chloride (MeHg, CAS# 115-09-3, Cl ∼13%) were dissolved in dimethyl sulfoxide (DMSO, CAS# 67-68-5). Stock solutions and final solutions were prepared in polypropylene vials. Final concentration of DMSO in the medium of all experimental conditions and controls was 0.25% (v/v).

NPC viability

Cell viability of compounds to NPCs was measured in a 5-day protocol (see also Fig. 1G). Briefly, when NPCs reached 100% confluency, cells were dissociated and seeded in a PLO-laminin coated 96-wells plate at a density of 1.25 × 105 cells/cm2 in NP medium. Cells were allowed to adhere to the plate for 2 h before exposure to either of the five compounds. NP medium was used as solvent control for ACR and VPA, NP medium with 0.25% DMSO was used as solvent control for CPF, FLX, and MeHg, 0.3 μg/mL 5-Fluorouracil (Sigma-Aldrich) was used as positive control (average cell viability 66%; 95% confidence interval 60–74%), and 0.5 mg/mL Penicillin G (Sigma-Aldrich) as negative control (average cell viability 92%; 95% confidence interval 85–97%) (both dissolved in medium). Three independent experiments, i.e. different vials of the same NPC bank started on different days, with six technical replicates per condition were performed. Exposure medium was refreshed on day 3. On day 5, cell viability was determined using the CellTiter-Blue® kit (Promega, Leiden, the Netherlands). Cells were incubated for 2 h and measured in a Spectramax® M2 spectrofluorometer (Molecular Devices, Berkshire, United Kingdom) at 544Ex/590Em.

Fig. 1

Neuronal differentiation of hESC-derived NPCs and establishment of hNPT. (A-D) Protein expression over time of a selection of markers for cell types: neurons (TUBB3, MAP2), astrocytes (GFAP), synapses (PSD95), and inhibitory (VGAT) neurotransmitter vesicles. (E) Gene expression relative to NPC culture of the same markers plus NES (neural progenitor), SYNPR (pre-synapse) and SLC17A6 (excitatory neurotransmitter vesicle). DLG4 encodes PSD-95, SLC32A1 encodes VGAT. N = 2 and n = 5-6. Filled and open circles indicate the two independent experiments. (F) Differentiation protocol of hNPT, and endpoints assessed at given time points. (G) Proliferation protocol of NPC, and endpoints assessed at given time points. Scale bar: (A,B) 100 μm, (C,D) 20 μm. White arrows in (A) indicate examples of GFAP + cells. Inserts (C, D) show details of PSD95 and VGAT. Statistical significance of gene expression changes (adjusted p-value) only relative to the NPC culture is indicated as *: p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001. Full statistics can be found in Supplementary Table 5.

hNPT viability

hNPT was prepared as mentioned above in PLO-laminin coated 48-well plates. Cells were exposed for ten days from the moment they were switched to differentiation medium. For the half medium refreshments on day 3 and 7, the same concentration of compound was added as on day 0. At least three experiments with three technical replicates per condition were performed. To determine cell viability, cells were incubated for 1 h with Cell Proliferation Reagent WST-1 (Roche, Mannheim, Germany) instead of the CellTiter-Blue® kit to avoid auto-fluorescence during fluorescence imaging.

hNPT neurite area

Immunocytochemistry was performed as described before (de Leeuw et al., 2020c). All washing steps were performed in 300 μL, all other steps in 150 μL. Briefly, after WST-1 measurements, hNPT was rinsed two times with warm DPBS (no calcium, no magnesium; Gibco) and fixed with 4% pre-warmed paraformaldehyde (Electron Microscopy Sciences, Hatfield, PA, USA) in DPBS for 30 min. Cells were then permeabilised with 0.2% Triton X-100 (Sigma-Aldrich) for in DPBS 5 min, washed two times 5 min with DPBS, and blocked for 1 h in blocking solution (1% bovine serum albumin (BSA; w/v; Sigma-Aldrich), 0.5% Tween-20 (Sigma-Aldrich) in DPBS). After washing, primary antibodies were applied overnight at 4 °C in 0.5% BSA/0.5% Tween-20 in DPBS (Supplementary Table 1). Samples were washed and incubated with secondary antibodies for 1 h in the same antibody incubation mixture (Supplementary Table 1). Nuclei were stained using DAPI (20 ng/mL, Sigma-Aldrich). Imaging for neurite area measurements was performed on a Leica DMi8 microscope system with 10x objective (Leica, Wetzlar, Germany) using the appropriate Leica Software (LAS X). Two independent experiments were performed with three wells per experimental condition, in which at least three images from each well were taken and used for image processing in Fiji/ImageJ (version 1.51n) (Schindelin et al., 2012) and CellProfiler (version 3.1.9) (McQuin et al., 2018). To measure neurite area, a mask was produced from the DAPI and MAP2 channel (fixed manual cut-off per experiment), which was then inverted to calculate the neurite area, normalised for the DAPI area. Additional high-resolution images to show synaptic markers were taken using a Nikon A1Rsi microscope (Nikon, Tokyo, Japan) using NIS Elements software.

Concentration-response curves of cell viability and neurite area were created using PROAST software (version 67.0) (Slob, 2002) in R (version 4.0.0) (R Core Team, 2020) using the exponential model. Effect concentrations with a 95% confidence interval were derived by PROAST from the resulting curves. Statistical analysis of differences between IC20 and IC5 values for each compound was performed using a two-way ANOVA test and post-hoc Sidak's multiple comparisons test using GraphPad Prism (version 9.1.0).

RNA isolation, qPCR and RNA-seq

RNA isolation

Samples of NPCs, control and exposed day 10 hNPT were fixed in QIAzol (Qiagen, Hilden, Germany). Medium was pipetted out of eight wells and QIAzol was applied onto the cells. Cells were resuspended by vigorous pipetting and transferred to vials on ice. Vials were then transferred and stored at −80 °C until further processing. Whole RNA extraction was performed following the manufacturer's protocol using the RNeasy® mini kit (Qiagen) with a DNase digestion step (Qiagen). RNA concentration was determined using the Qubit3 (Invitrogen, Carlsbad, CA, USA) and RNA quality was analysed using the 2100 Bioanalyzer (Agilent Technologies, Amstelveen, the Netherlands). For construction of the time line using qPCR, five or six samples per time point were analysed from each of two experiments over the course of differentiation on day 0, 1, 4, 7, 10, 14 and 28. For RNA-seq analysis, eight replicates from a single separate experiment were taken from the compound exposure on day 10, plus eight replicate samples of unexposed NPCs and hNPT. Throughout the manuscript, “N” reflects the number of independent experiments and “n” reflects the number of biological replicates.

qPCR

Synthesis of cDNA was performed using the high capacity cDNA reverse transcription kit containing random hexamer primers (Applied Biosystems, Foster City, CA, USA). Ingredients from the kit were prepared according to manufacturer's instructions and mixed with 400 ng mRNA per sample in 25 μL demineralised water (dH2O), resulting in a total volume of 50 μL. 3 μL diluted cDNA (3x in dH2O) was combined with 0.5 μL primer (Supplementary Table 2), 5 μL TaqMan Fast Universal PCR Master Mix (no AmpErase™ UNG; Applied Biosystems) and 1.5 μL dH2O. Gene quantification was done on a 7500 Fast Real-Time PCR system (Applied Biosystems) in optical 96-well fast thermal cycling plates (Applied Biosystems) with the following thermal cycling conditions: 95 °C for 20 s, 40 cycles of 95 °C for 3 s, 60 °C for 30 s. Relative gene expression differences were calculated using the 2−ΔΔCt -method (Applied Biosystems, 2001), normalised against the housekeeping genes Hypoxanthine phosphoribosyltransferase 1 (HPRT1), and Glucuronidase beta (GUSB). Statistical analysis was performed using a nested one-way ANOVA test and post-hoc Tukey's multiple comparisons test for each gene in GraphPad Prism (version 9.1.0).

RNA-seq

All methodological information of the RNA-seq analysis is summarised in Supplementary Table 3. Samples were processed according to the TruSeq Stranded mRNA protocol (Illumina, San Diego, CA, USA). Enrichment of mRNA was done using polyA-affinity purification converting it into a multiplex library of 16 samples per flow cell. Sequencing was performed on the NextSeq 500 sequencer (Illumina) using the NextSeq 500/550 High Output Kit version 2.5 (75 cycles). Raw bcl files were base called and demultiplexed into FASTQ files with Bcl2fastq (Illumina, version 2.20.0.422). Quality control reports (FastQC, version 0.11.8) were generated using MultiQC (version 1.0) (Ewels et al., 2016) to identify potential anomalies regarding library size, read length distribution, mean read quality distribution, mean quality for each position in the read, and base frequency for each position in the read. This revealed a GC- and a fragment length bias between the multiplex runs. Expression of the transcripts was quantified against the human transcriptome (GRCh38, version 30, Ensembl 96) using Salmon (version 1.3.0) (Patro et al., 2017). The estimated transcript abundances were imported into RStudio (version 1.3.959) (RStudio Team, 2020) and summarised into genes using tximport (version 1.16.1) (Soneson et al., 2016). Gene expression data have been deposited in the NCBI GEO database under accession number GSE166297.

Principal component analysis (PCA) was performed on variance stabilizing transformation (vst)-normalised data to identify potential outliers. The PCA analysis revealed one outlier (MeHg IC20/100 sample 8), which was due to the high sequence duplication levels seen in the raw fastq data (59.3%), and was therefore left out of further analysis. Differential gene expression between conditions was calculated using the DESeq2 package (version 1.30.0) (Love et al., 2014), adjusting for batch and using the lfcShrink ashr function to extract final log2 fold changes (log2FC) (Stephens, 2017). Differentially expressed genes (DEGs) were obtained by filtering results for an absolute p-value ≤0.001 and log2FC ≥ 1 for the comparison between hNPT and NPCs, and an absolute p-value ≤0.001 and log2FC ≥ 0.5 for the comparison between compounds. These cut-off values were chosen to keep statistically robust changes in gene expression without omitting more subtle gene expression changes that participate in biological pathways. Functional enrichment of genes in Gene Ontology (GO-) terms (Biological processes, GO FAT) was performed using DAVID (consulted on 8 November 2020, version 6.8) (Huang et al., 2009) against a background of genes with at least one count in the analysis (38,078 genes for time, 41,485 for compounds). GO-terms with false discovery rate (FDR) ≤0.05 were considered to be enriched. Categorisation of GO-terms was performed based on their GO classification (http://amigo.geneontology.org/amigo) for producing the pie charts summarizing all GO-terms. The terms relating to neurons were further categorised under one of the biological processes predefined in Hessel et al. (2018). Interactions between genes were obtained with STRING (consulted on 15 November 2020, version 11.0) (Szklarczyk et al., 2019).

Pearson's correlation coefficients comparing log2FC values from PCR and RNA-seq data for a selection of genes showed a strong to very strong correlation, which demonstrates the validity of the RNA-seq data (Supplementary Table 4).

Data visualisation

Graphs were visualised using GraphPad Prism (version 9.1.0), heatmaps (ward.D clustering) and PCA plots were built in RStudio (version 4.0.0) (R Core Team, 2020), the cell lineage map was made in Pathvisio (version 3.3.0) (Kutmon et al., 2015), Venn diagrams were produced using Venny (version 2.1.0) (Oliveros, 2007) and gene interactions from STRING (full network, medium confidence) were visualised in Cytoscape (version 3.7.4) (Shannon et al., 2003), using Perfuse force-directed layout.

Results

hNPT consisted of excitatory neurons, inhibitory neurons, and astrocytes after ten days of differentiation

In a previous study, we have characterised a neuronal-astroglial differentiation protocol in which we found excitatory neurons, inhibitory neurons, and astrocytes after two weeks of culturing (de Leeuw et al., 2020c). Spontaneous electrical activity was measurable already after a few days of differentiation, suggesting that a network was already forming within these two weeks. For the development of the human neural progenitor test (hNPT), we therefore followed the culture in more detail in the first two weeks to find the moment at which these three cell types were present.

Immunocytochemistry between day 0 and 14 revealed that neuronal markers TUBB3 and MAP2 were already present in some of the cells in the neural progenitor cells (NPC) phase (Fig. 1A–D). The same was true for protein expression of synaptic markers SYNPR and PSD95 (encoded by DLG4), and vesicle marker VGLUT2 (encoded by SLC17A6) (PSD95: Fig. 1C, SYNPR, VGLUT2: de Leeuw et al., 2020c). Astrocyte marker GFAP emerged from day 7 (Fig. 1A). Vesicle marker VGAT (encoded by SLC32A1) had a similar delay in protein expression as GFAP, only being present from day 7 (Fig. 1D). Gene expression of the same set of markers showed that all markers related to functional neurons and astrocytes were upregulated between day 4 and day 14 and most of them plateaued from that time point (Fig. 1E). Neural progenitor marker NES was the only downregulated gene, consistent with the differentiation occurring in hNPT. Together, these results led to the establishment of hNPT as a 10-day differentiation protocol for exposure experiments (Fig. 1F).

RNA-seq revealed neurodevelopmental processes and cell types present in hNPT

To define which cell types and neurodevelopmental processes were regulated during differentiation, RNA-seq was performed on the NPCs and hNPT. In total, 3992 differentially expressed genes (DEGs) were regulated over time (2155 up, 1837 down), of which 2343 genes were linked to GO-terms by DAVID. Analysis of selective markers for neuronal cell types using our previously described cell lineage map (de Leeuw et al., 2020a) showed an upregulation of glutamatergic, GABA-ergic, and glycinergic neuron markers, mixed regulation of NPC, neuronal, astroglial, and synaptic component markers, and downregulation of oligodendrocyte markers (Fig. 2A). Less and mixed regulation was observed among cholinergic, dopaminergic, and serotonergic markers.

Fig. 2

Cell types present and GO-terms regulated in hNPT. (A) Cell lineage map (de Leeuw et al., 2020a) showing gene expression regulation (log2 fold change (log2FC)) of selective markers for a range of cell types along the ectodermal and neuronal lineage in hNPT relative to the NPC culture. Light green boxes indicate a p-value of ≤0.001, dark green boxes indicate when a gene was additionally regulated more than 1 log2FC. (B) Pie charts list GO-terms enriched by regulated genes summarised per category according to their GO classification (http://amigo.geneontology.org/amigo). Synaptic communication included terms from Neurotransmission in Table 1, Neurodevelopment included terms from Cell differentiation and development, and Neuronal differentiation included terms from Neurite outgrowth and synapse formation. Samples for these graphs were hNPT day 10 and NPC (“day 0”) control samples (N = 1 and n = 8). (C) Top-10 GO-terms with highest fold change enrichment scores calculated by DAVID. Numbers between brackets present the number of genes regulated in the GO-term and the total number of genes in the GO-term, respectively. All individual GO-terms with classification and statistics can be found in Supplementary Table 6. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1

Number of neuronal-related GO-terms and number of regulated genes involved in hNPT development relative to the NPCs, grouped to neurodevelopmental processes as defined by Wegner et al. (2020) and Hessel et al. (2018). All individual GO-terms with classification and statistics are listed in Supplementary Table 6. Samples were hNPT day 10 and NPC (“day 0”) control samples (N = 1 and n = 8). Percentages in brackets indicate which part of the GO-terms or genes of the total were regulated. “-” indicates that there were no GO-terms regulated.

Neurodevelopmental processes
Wegner et al. (2020)	Hessel et al. (2018)	GO-terms regulated genes (471)	No. of regulated genes (2343)
Proliferation & stem cell maintenance	Neural progenitor cell formation	–	–
	Neural progenitor cell proliferation	3 (0.6%)	38 (1.6%)
Cell differentiation and development	Neurogenesis	17 (3.6%)	589 (25.1%)
	Neurons phenotypic specification	8 (1.6%)	275 (11.7%)
Gliogenesis and myelination	Gliogenesis	2 (0.4%)	56 (2.4%)
	Myelination	–	–
Neural migration	Neuron migration	–	–
Neurite outgrowth and synapse formation	Neurite outgrowth	5 (1.1%)	188 (8.0%)
	Dendritogenesis/dendrite formation	–	–
	Axon formation	2 (0.4%)	101 (4.3%)
	Axon guidance/axon path finding	1 (0.2%)	56 (2.4%)
	Synaptogenesis/synapse formation	5 (1.1%)	57 (2.4%)
Neurotransmission	Synaptic connectivity	16 (3.4%)	159 (6.8%)
	Synaptic pruning	–	–
	Neural network functioning/synaptic plasticity	1 (0.2%)	42 (8.9%)
Apoptosis
Regional specificity	Anatomical	11 (2.3%)	202 (8.6%)

Enrichment analysis of the regulated genes resulted in 471 regulated Gene Ontology (GO-) terms (Fig. 2B, Supplementary Table 6). 14% of GO-terms enriched were of a neuronal nature, while 18% of the GO-terms enriched were related to other differentiation terms. The other 68% were related to general cell processes. Still, among the most regulated GO-terms (i.e. with the lowest fold change enrichment score) were specific terms relating to glutamatergic and GABA-ergic synaptic signalling (Fig. 2C). Interestingly, the top regulated GO-term was “regulation of endodermal cell differentiation”, of which all genes but one were downregulated. These neuronal GO-terms were further categorised in a predefined list of important neurodevelopmental processes. We used a categorisation published by Wegner et al. (2020) that was previously applied to a neuronal differentiation protocol and linked these to a more detailed list of important neurodevelopmental processes defined by Hessel et al. (2018) (Table 1). DEGs and related GO-terms clustered in neural progenitor cell proliferation (3 GO-terms, e.g. ‘positive regulation of neural precursor cell proliferation’), neurogenesis (17 GO-terms, e.g. ‘nervous system development’, ‘neuron differentiation’), gliogenesis (2 GO-terms, e.g. ‘glial cell differentiation’), axonal outgrowth, and synapse formation (13 GO-terms, e.g. ‘axon development’, ‘axon guidance’, ‘synapse organisation’), neurotransmission (17 GO-terms, e.g. ‘glutamatergic synaptic transmission’, ‘GABA-ergic synaptic transmission’, ‘regulation of synaptic plasticity’), and regional specificity (11 GO-terms, e.g. ‘forebrain development’, ‘hindbrain development’). There were no GO-terms related to neuronal migration, synaptic pruning or neuronal apoptosis regulated in hNPT.

In summary, in hNPT, NPC differentiated in a mixed culture of excitatory neurons, inhibitory neurons, and astrocytes in which GO-terms for essential neurodevelopmental processes were regulated such as neuronal and glial differentiation, axon guidance, synaptogenesis, and neurotransmission.

Exposure to five DNT compounds affected cell viability and neurite area at similar concentrations

We next exposed hNPT to either of five known or suspected DNT compounds: acrylamide (ACR), chlorpyrifos (CPF), fluoxetine (FLX), methyl mercury (MeHg), or valproic acid (VPA). Both cell viability and neurite area (Fig. 3A) were measured, and NPC cell viability was measured in a separate, shorter 5-day exposure. For each compound an IC5 (just below effects on viability) and IC20 divided by 100 (IC20/100, no effects on viability) concentration was derived (Pistollato et al., 2020). The methods to determine cell viability were all comparable since no statistically different inhibitory concentrations (IC) were found between any of the three applied methods, except for a lack of NPC effects on cell viability by VPA exposure (Fig. 3B, Fig. S1). Due to the ease of the readout relative to ‘neurite area’, we chose hNPT cell viability measures to derive the IC5 and IC20/100 concentrations for RNA-seq experiments (Table 2). Immunocytochemistry showed that using these concentrations resulted in similar morphology of all experimental conditions relative to the control condition, except for ACR IC5, which showed signs of cytotoxicity (Fig. 3C).

Fig. 3

Compound effects on NPCs and hNPT. (A) Summary of toxicity outcomes for all compounds. Error bars represent 95% confidence interval for NPC viability (N = 3 and n = 6), hNPT viability (N = 3 and n = 3) and neurite area (N = 2 and n = 3). Red triangles represent mean concentration in umbilical cord blood as measured in Hogervorst et al. (2021) (ACR), Silver et al. (2018) (CPF), Laine et al. (2003) (FLX), HBM4EU dashboard ((VITO), n.d.) (MeHg), Bank et al. (2017) (VPA). (B) Graphical flow of CellProfiler image processing pipeline for calculating the neurite area. (C) Protein expression of cell-type markers for NPCs (PAX6), neurons (MAP2), and astrocytes (GFAP) in experimental conditions used for RNA-seq. Significance level: * = adjusted p-value ≤0.05. Scale bar: 200 μm. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 2

IC5 and IC20/100 values, and selected compound concentrations for RNA-seq experiments. All values are given in μM. N = 3 and n = 3 for all experimental conditions.

	IC5 (95% confidence interval)	Chosen concentration	IC20/100 (95% confidence interval)	Chosen concentration
ACR	320 (130–500)	450	5.4 (3.0–7.3)	6.5
CPF	21 (15–33)	15	0.36 (0.27–0.48)	0.3
FLX	9.5 (4.1–6.4)	10	0.14 (0.077–0.11)	0.12
MeHg	0.11 (0.040–0.14)	0.1	0.0016 (0.00083–0.0024)	0.0015
VPA	180 (57–530)	250	5.3 (3.3–8.0)	6

Exposure to five DNT compounds caused global gene expression changes and regulated common and distinct GO-terms in hNPT

RNA-seq analysis of exposed hNPT revealed that compounds regulated between 21 and 2764 DEGs (p < 0.001 and log2FC ≥ 0.5, Fig. 4A). These large differences in DEG numbers were also reflected in a principal component analysis (PCA) and a heatmap based on the 3737 uniquely regulated DEGs (Fig. 4B and C). For example, IC5 of ACR and VPA showed distinct gene expression regulation patterns relative to each other and the other compounds, while CPF, FLX, and MeHg clustered together, showing more similar gene expression regulation patterns.

Fig. 4

Global gene expression changes in hNPT in response to compound exposure. (A) Number of genes up and downregulated by each of the compounds and compound concentrations (p < 0.001 and log2FC ≥ 0.5). (B) PCA plot of gene expression after exposure to compounds and the control condition based on 3737 genes that are differentially regulated in any of the IC5 exposure conditions. (C) Same data as in (B) presented in a heatmap. ACR: acrylamide, CPF: chlorpyrifos, FLX: fluoxetine, MeHg: methyl mercury, VPA: valproic acid. PC: principal component. N = 1 and n = 8 for all experimental conditions.

Enrichment of regulated genes revealed a distinct pattern of regulated GO-terms for each of the compounds. ACR and VPA IC20/100 did not enrich any GO-term. The other IC20/100 concentrations regulated mostly general cell processes and few differentiation GO-terms, while IC5 concentrations elicited specific neuronal GO-terms as well (Fig. 5). In the case of CPF and FLX, the IC5 enriched additional GO-terms relative to the IC20/100, while in the case of MeHg the type rather than the amount of GO-terms changed to a neuronal nature after exposure to the IC5. Interestingly, while ACR IC5 regulated most genes, these genes converged in the lowest number of GO-terms.

Fig. 5

Summary of GO-terms regulated by each experimental condition. Pie chart list GO-terms summarised per category according to their GO classification (http://amigo.geneontology.org/amigo). Categories as in Fig. 2. N = 1 and n = 8 for all experimental conditions. All individual GO-terms with classification and statistics can be found in Supplementary Table 7.

Categorisation of the neuronal terms in neurodevelopmental processes as done in Table 1 resulted in distinct but overlapping combinations of GO-terms regulated by each of the compounds. FLX IC20/100 did not enrich any neuronal-specific GO-terms. At least one of the concentrations of each compound regulated GO-terms related to nervous system development, indicating that all compounds had the potential to cause DNT in hNPT. While neuronal apoptosis was not among the GO-terms regulated over time, it was regulated by CPF and by FLX IC5.

Regulated neurodevelopmental processes revealed common biomarker genes and specificity of compounds

Three of the enriched GO-terms related to specific neurodevelopmental processes were shared by at least three compounds: ‘axon development’ (GO:0061564, 443 genes), ‘neuron apoptosis’ (GO:0051402, 202 genes), and ‘synaptic signalling’ (GO:0099536, 587 genes, the other three synaptic terms contain the same genes) (Table 3). ‘Axon development’ and ‘neuron apoptosis’ had respectively eight and seven genes in common among the compounds that regulated these GO-terms (Fig. 6A, C). The regulation of these genes appeared highly similar across all compounds, except for ACR IC20/100 and VPA IC20/100 (Fig. 6B, D). Within the GO-term ‘neuron apoptosis’, MeHg IC5 also presented a different gene expression regulation pattern relative to the other compounds. UNC5B was regulated in both GO-terms (Fig. 6B, D).

Table 3

Denomination of GO-terms relating to neurodevelopmental processes regulated by DNT compounds. All individual GO-terms with classification and statistics are listed in Supplementary Table 7. N = 1 and n = 8 for all experimental conditions. GO-terms were summarised under the broadest denominator when they directly related to each other in Gene Ontology. Numbers in brackets indicate the total number of GO-terms per category and number of genes identified by DAVID, respectively. “-“ indicates that there were no GO-terms regulated.

	ACR IC5	CPF IC20/100	CPF IC5	FLX IC5	MeHg IC20/100	MeHg IC5	VPA IC5
Total no. of regulated GO terms	33	50	113	136	55	68	122
Total no. of unique genes in regulated GO-terms	949	191	281	230	128	165	559
Proliferation & stem cell maintenance	–	–	–	–	Neural progenitor proliferation (2, 9)	–	–
Cell differentiation and development	Nervous system development (4, 298)	–	Nervous system development (3, 74)	Nervous system development (6, 65)	–	Nervous system development (9, 64)	Nervous system development (4, 144)
Gliogenesis and myelination	–	–	–	Glial cell differentiation (2, 11)	–	–	–
Neural migration	–	–	–	–	–	–	–
Neurite outgrowth and synapse formation	–	–	Neuron projection/axon development and guidance (1, 24)	Neuron projection/axon development and guidance (3, 19)	–	Neuron projection/axon development and guidance (8, 32)	–
Neurotransmission	Synaptic signalling (4, 97)	–	–	Synaptic signalling (4, 27)Glutamate receptor signalling (1, 8)Regulation of neurotransmitter levels (2, 18)	–	Synaptic signalling (5, 29)Glutamate receptor signalling (1, 8)	Synaptic signalling (5, 52)Neurotransmitter transport (1, 23)Regulation of neurotransmitter levels (1, 23)
Apoptosis	–	Neuron apoptosis (1, 12)	Neuron apoptosis (2, 15)	Neuron apoptosis (2, 15)	–	–	–
Regional specificity

Fig. 6

Comparison of common genes affected by different compound exposures. (A) Common genes among CPF IC5, FLX IC5, and MeHg IC5 for the GO-term ‘Axon development’. (B) Expression of the seven common genes (log2FC) of (A) in all experimental conditions. (C) Common genes among CPF IC20/100, CPF IC5, and FLX IC5 for the GO-term ‘Neuron apoptotic process’. (D) Expression of the eight common genes (log2FC) of (C) in all experimental conditions. N = 1 and n = 8 for all experimental conditions.

‘Synaptic signalling’, regulated in four experimental conditions, showed overlap between four test conditions of only two genes (HRH2 and SLC6A9; Fig. 7A), which were also statistically significantly regulated in CPF IC5. Further investigation was performed on all regulated genes in this GO-term by analysing their interaction in STRING. This revealed a network of gene clusters that were specific for eight major neuron subtypes (Fig. 7B). Projecting gene expression regulation of each compound showed their diverse action on this gene set. ACR IC5 regulated the majority of the genes, mostly downwards, across almost all neuron subtypes (Fig. 7C). VPA IC5, in contrast, mostly upregulated genes and specifically regulated the serotonergic and dopaminergic markers in the network. MeHg IC5 showed a selective downregulation of genes, among others five of the nineteen glutamatergic markers. FLX IC5 showed a similar pattern, although the direction of regulation was more mixed. Analysis of CPF IC5, which did not show regulation of this GO-term, altered gene expression of only nineteen genes which did not focus on a specific neuron subtype.

Fig. 7

Gene expression regulation patterns by compound exposures within the GO-term ‘synaptic signalling’. (A) Venn diagram of four experimental conditions showing the overlap of genes in the GO-term. (B) Gene interaction network of 128 of the 138 genes in (A) that had a connection with each other in STRING, visualised in Cytoscape. Genes that were not connected are not shown in the image. Colours indicate neuron subtype specificity of genes. (C) Expression of genes in (A) projected on the network in (B), per compound, and CPF IC5 expression of 19 DEGs that were present among the 128 genes. N = 1 and n = 8 for all experimental conditions. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Together, these results suggest that hNPT may be able to both detect DNT compounds by monitoring a limited number of genes, as well as to show the specificity of DNT compounds on one of the dominant and most complex processes in hNPT, which is synaptic signalling.

Discussion

hNPT is a 10-day protocol that has been shown to present functional electrical activity between neurons (de Leeuw et al., 2020c). This current study indicates using RNA-seq that various processes essential for the development of the central nervous system over time are represented in hNPT. These include neuronal and glial differentiation, neurite extension, synaptogenesis, and synaptic transmission. Exposure to known DNT compounds regulated GO-terms related to neural progenitor cell proliferation, neuronal and glial differentiation, axon development, synaptogenesis, synaptic transmission, and apoptosis. These results infer that hNPT might be a possible model to study neuron differentiation and maturation, but also to study more complex processes such as synaptogenesis and neuronal network functioning.

The field of developmental neurotoxicity (DNT) needs fast, reliable in vitro models that truthfully resemble human neurodevelopmental processes. OECD is currently developing a guidance document for testing and assessment of an in vitro battery for DNT (Hernández-Jerez et al., 2021; Masjosthusmann et al., 2020; Sachana et al., 2021a). While there are already many cell models available representing earlier time windows in neurodevelopment, such as NPC proliferation, migration, differentiation in various cell types and neurite outgrowth (reviewed in Bal-Price et al., 2018a; Fritsche et al., 2018), models including markers for synaptogenesis and neural network formation are currently sparse and still mainly based on rodent models (Sachana et al., 2019, 2021a). The current results show that this hNPT model may be an interesting model to add to the available in vitro battery due to the short protocol and the complex processes represented. As this was a study based on a single RNA-seq experiment, Further studies will need to refine the model, specifically regarding readout parameters for more detailed analysis of synaptogenesis and neuronal function e.g. using the microelectrode array (mwMEA), in addition to related gene and protein biomarkers.

There are multiple promising human embryonic and induced pluripotent stem cell (hiPSC)-based models that have recently been developed for the combined study of neurite outgrowth, synaptogenesis, and neurotransmission (Barenys et al., 2016; Baumann et al., 2016; Bu et al., 2020; Kobolak et al., 2020; Nimtz et al., 2020; Pistollato et al., 2020; Wegner et al., 2020; Zhong et al., 2020). The advantage of hNPT over other currently available protocols is the short duration for these processes to occur and the concurrent development of synaptic activity in the same period (de Leeuw et al., 2020c). There is one protocol for 3D neurospheres that showed a similar development of key genes as hNPT (Sandström et al., 2017). While 3D cultures have the obvious advantage that they more reliably mimic the in vivo architecture, 2D cultures benefit from the ease of readout and by extension throughput possibilities (Jensen and Teng, 2020; Pacitti et al., 2019; Schmidt et al., 2017). A potential reason why hNPT is developing fast may be because the starting material consists of both NPC and young neurons, as shown by the immunostainings. With regard to neurite outgrowth, multiple robust protocols already exist for this endpoint (Harrill et al., 2011; Li et al., 2021; Ryan et al., 2016; Stiegler et al., 2011). The hNPT model might be useful to study this endpoint as well, but the current study was focussed on detailed molecular readouts based on gene expression analysis. Follow-up studies using clear endpoint-specific toxicants for neurite outgrowth (such as rotenone and colchicine) and optimisation of the image analysis are needed to further explore this endpoint.

Transcriptomics represents an unbiased, relatively high-throughput and sensitive measure of molecular changes, by which effects may be detected at concentrations below those causing viability effects. In the present study, the concentration at which effects on neuronal development became apparent was just below those inducing overt decrease in cell viability (Pistollato et al., 2020). This is in line with some studies in which molecular effects of compounds were apparent just below effects on cell viability were present (Chen et al., 2019; Krug et al., 2013). Other studies, however, have shown robust transcriptomic effects already at concentrations considerably lower than those at which cytotoxicity occurred (Theunissen et al., 2012; Waldmann et al., 2014). This may be due to the resilience of this specific cell culture the nature of the compounds, mainly early teratogens, that cause widespread differentiation changes as opposed to quite selective synaptic changes by most of these compounds.

RNA-seq analysis over time revealed gene expression regulation of expected processes such as progenitor proliferation, neurite development (both dendrites and axons), glial differentiation, synaptogenesis and neurotransmission, a broad array of important neurodevelopmental processes as defined by many authors (Bal-Price et al., 2018b; Hessel et al., 2018; Li et al., 2018; Rice and Barone, 2000; Silbereis et al., 2016; Wegner et al., 2020). We did not observe gene expression changes related to neuron migration, which may suggest that the neurons grew within colonies of cells instead of migrating from these clusters. Glial differentiation-related GO-terms were regulated, but not as many as for neuronal differentiation terms. This is in line with the immunostainings where glia were present at day 10 at far fewer numbers than neurons. Longer culturing may increase the strength of the astroglial differentiation GO-terms in hNPT.

Using the same framework as for neurodevelopmental processes over time, compounds showed a unique combination of affected GO-terms. Low compound concentrations affected either none (ACR, VPA) or mainly general cell function GO-terms (CPF, FLX, MeHg), while high compound concentrations elicited specific neuronal-related effects reminiscent of published findings by other groups. ACR, for example, is a neurotoxicant that affects the developing cerebellum and peripheral motor system in vivo (Lindeman et al., 2021; LoPachin et al., 2003), which may explain the lack of effect in the low concentration of ACR on hNPT and other CNS models (Hoelting et al., 2016; Ryan et al., 2016). ACR IC5, specifically affected genes relating to general neurogenesis and synaptic signalling, which is comparable to impaired synaptic function by ACR exposure in the CNS (Barber et al., 2007; LoPachin and Gavin, 2012). However, it should be noted here that the tested ACR IC5 concentration appeared to reduce cell viability more than expected, therefore gene expression changes may show potential interference of cytotoxicity. CPF only had a specific effect on a neurite outgrowth GO-term, an effect also found in other in vitro studies (Pistollato et al., 2020; Qiao et al., 2001; Wu et al., 2017). The DNT effects of FLX are still under debate (Campagne, 2019), but this compound was previously shown to affect synaptic signalling in rodent in vitro and in vivo models (Chang et al., 2015; Frank et al., 2017; Mack et al., 2014), which was also the case in hNPT. MeHg is a classic neurotoxicant that affects neuron extension and glutamatergic neurotransmission (Farina et al., 2011; Fujimura and Usuki, 2015; Ryan et al., 2016), in line with enriched GO-terms in hNPT (i.e. neuron projection/axon development and guidance, synaptic signalling and glutamate receptor signalling). VPA is known for its widespread effects on differentiation, as was also apparent from the relatively large proportion of general development and differentiation GO-terms affected, as well as neuronal-specific terms. These latter terms were mainly focused around synaptic signalling in hNPT, which is in line with effects on the excitatory/inhibitory balance in rodents (reviewed in Tartaglione et al. (2019)). In short, hNPT mimicked some of the key effects of the tested DNT compounds in a dose dependent manner. This offers opportunities in hNPT to unravel physiological responses from adverse compound effects and to further explore the biological domain of the model.

Some of the expected effects of compounds were not detected by gene expression analysis in hNPT, which can mostly be traced back to the lack of the specific biological processes or cell types involved. For example, the effect of FLX on serotonergic markers, dopaminergic neuron differentiation, and suspected effects on oligodendrocytes (de Leeuw et al., 2020a; Jha et al., 2016; Kroeze et al., 2016; Lupu et al., 2018; Sommi et al., 1987) were not detected, probably due to the absence of these specific cell types. Similarly, CPF did not have expected effects on dopaminergic, serotonergic, and cholinergic neurotransmission (Slotkin et al., 2012; Todd et al., 2020). MeHg is known to affect neuron migration in the developing cortex (Antunes dos Santos et al., 2016), but migration was not affected in hNPT, in line with its absence in this model. There were also a number of processes present in hNPT that were not disrupted by the tested compounds. For example, FLX decreases cortex thickness in rodents (Swerts et al., 2010), which was not among the upregulated GO-terms in this study. Also, the effect of CPF on synapse formation found in similar in vitro models (Pistollato et al., 2020; Sandoval et al., 2019) was not observed in hNPT. Therefore it is essential to take into account the biological domain of the model and to understand that an in vitro model can only partly represent the in vivo brain development.

Upon further investigation of a selection of common GO-terms among compounds, we found a set of common genes for axon development and neuron apoptosis, indicating conserved pathways. The genes involved in axon development encode a chemorepellant (SEMA3D), receptors for diverse chemoattractors or repellants (EPHA8, UNC5B), transcription factors and binding protein for motor axon outgrowth (CRABP2, ETV1, ETV4), a stabiliser of myelin (PLP1) (Bashaw and Klein, 2010; Diehl et al., 1986; Han et al., 2018). UNC5B is also involved in external apoptosis, as are DDIT3, CEBPB, BBC3, PMAIP1, CLN3, and PRKCG (Persaud-Sawin et al., 2002; Thiébault et al., 2003; Yang et al., 2017). AARS1 is involved in tRNA synthesis (Shiba et al., 1995). These lists of genes might serve as a biomarker set for these specific processes or even as a general gene set for DNT detection. It is interesting to note that UNC5B is upregulated in virtually all exposure conditions and is involved in both axon development and apoptosis. This paradoxical function can be explained by the fact that UNC5B encodes a dependence receptor of netrin-1, which means that the receptor can have both pro- or anti-apoptosis effects depending on the presence of netrin-1 (Hofmann and Tschopp, 1995; Ko et al., 2012). UNC5B is suggested to be involved in microtubule reorganisation, which is needed for axon guidance but can also be an early sign of apoptosis (Kellermeyer et al., 2018).

The genes involved in synaptic signalling, on the other hand, showed little overlap both in occurrence and direction of effect. This exemplifies the complexity of synaptic transmission as opposed to the other two processes (Grant, 2019) and may provide further mechanistic insight into how compounds may affect these processes. Follow-up experiments focussing on functional parameters and protein expression using a wider range of compounds should be performed to confirm these results.

In conclusion, this first RNA-seq study aimed to set the boundaries of the biological domain of the hNPT by investigating which processes are represented in the model and which processes can be disrupted by chemical exposure. Future efforts in hNPT may focus on further characterisation of the model and benchmarking it against in vivo data (Wegner et al., 2020) for which various tools are available (Stein et al., 2014; van de Leemput et al., 2014) to further define its biological domain. This is important in the light of studying the molecular basis of consequences of compound exposure on cognitive development, of which effects can be very time and region-dependent (Carlson et al., 2020; Prem et al., 2020). Applying the hNPT protocol on hiPSC is also an important future step due to the ethical issues surrounding the use of hESC. Previous research from our research lab has shown that this should in principle be feasible (de Leeuw et al., 2020b). Furthermore, a quantitative readout on protein level, such as with flow cytometry, can enhance the interpretation of the consequences of transcriptomics findings from this study, and to investigate cell types specifically affected by compounds. Testing more compounds with diverse modes of action, both neurotoxic and non-neurotoxic, can also aid in mapping the boundaries and potential of hNPT to detect effects on various neurodevelopmental processes (Aschner et al., 2017; Harrill et al., 2018; Sund and Deceuninck, 2021). Knowing its biological domain, hNPT might have added value as a tool for neurodevelopmental toxicity testing in vitro.

Funding statement

This research is funded by the Dutch NGO Stichting Proefdiervrij, the Dutch Ministry of Agriculture, Nature and Food Quality, the Dutch Ministry of Health, Welfare and Sports, and European Union’s Horizon 2020 research and innovation programme under Grant Agreement 874583 (ATHLETE Project). This publication reflects only the authors' view and the European Commission is not responsible for any use that may be made of the information it contains.

Credit author statement

Victoria de Leeuw: Investigation, Formal analysis, Visualisation, Writing – original draft. Conny van Oostrom: Investigation. Paul Wackers: Formal analysis, Data Curation. Jeroen Pennings: Investigation. Hennie Hodemaekers: Investigation. Aldert Piersma: Supervision, Conceptualisation, Funding acquisition. Ellen Hessel: Supervision, Conceptualisation, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.