Genome-wide expression screening in the cardiac embryonic stem cell test shows additional differentiation routes that are regulated by morpholines and piperidines.

The cardiac embryonic stem cell test (ESTc) is a well-studied non-animal alternative test method based on cardiac cell differentiation inhibition as a measure for developmental toxicity of tested chemicals. In the ESTc, a heterogenic cell population is generated besides cardiomyocytes. Using the full biological domain of ESTc may improve the sensitivity of the test system, possibly broadening the range of chemicals for which developmental effects can be detected in the test. In order to improve our knowledge of the biological and chemical applicability domains of the ESTc, we applied a hypothesis-generating data-driven approach on control samples as follows. A genome-wide expression screening was performed, using Next Generation Sequencing (NGS), to map the range of developmental pathways in the ESTc and to search for a predictive embryotoxicity biomarker profile, instead of the conventional read-out of beating cardiomyocytes. The detected developmental pathways included circulatory system development, skeletal system development, heart development, muscle and organ tissue development, and nervous system and cell development. Two pesticidal chemical classes, the morpholines and piperidines, were assessed for perturbation of differentiation in the ESTc using NGS. In addition to the anticipated impact on cardiomyocyte differentiation, the other developmental pathways were also regulated, in a concentration-response fashion. Despite the structural differences between the morpholine and piperidine pairs, their gene expression effect patterns were largely comparable. In addition, some chemical-specific gene regulation was also observed, which may help with future mechanistic understanding of specific effects with individual test compounds. These similar and unique regulations of gene expression profiles by the test compounds, adds to our knowledge of the chemical applicability domain, specificity and sensitivity of the ESTc. Knowledge of both the biological and chemical applicability domain contributes to the optimal placement of ESTc in test batteries and in Integrated Approaches to Testing and Assessment (IATA).

Introduction

The extensive use of animals under current international regulations for human chemical safety testing is increasingly conflicting with ethical and scientific principles. This is especially relevant in the field of developmental and reproductive toxicology (DART), as it uses relatively large numbers of animals partially because effects on multiple generations are assessed (van der Jagt et al., 2004, Rovida and Hartung, 2009, Beekhuijzen, 2017). To move away from animal-model based ‘black box’ testing and to focus more directly on human biology, alternative methods should preferably be mechanism based as this would increase comparability between perturbation of developmental pathways between species and would facilitate extrapolation of laboratory test methods to human individuals (Gibb, 2008).

The cardiac embryonic stem cell test (ESTc) is a well-studied in vitro assay for developmental toxicity testing. It determines chemical-induced perturbations of the differentiation of pluripotent stem cells into beating cardiomyocytes (Genschow et al., 2004). Within the ESTc, a heterogenic cell population is generated with mixed cell types besides cardiomyocytes. For example, neural crest cells and neurons are also present within the ESTc (Mennen et al., 2019, Mennen et al., 2021a). However, a complete inventory of the cell type composition generated by stem cell differentiation within the ESTc has not been mapped so far. Such an understanding of the complete biological domain of the ESTc could improve our mechanistic understanding of cell differentiation in the ESTc and provide more information on which mechanisms within the ESTc can be perturbed, aiding the specificity and the sensitivity of this in vitro tool for the assessment of chemically induced developmental toxicity.

It is already known that not all developmental toxicants that show an in vivo response in laboratory mammals also show a response in the ESTc. While this difference in assay sensitivity is not yet well understood, it is logical given that the stem cells cannot fully mimic the biological complexity and diversity in whole organisms. The reason for these differences is becoming better understood using the growing knowledge of toxicity mechanisms at the sub-cellular level e.g. using gene expression profiling. This would explain why the original ESTc method by scoring beating cardiomyocytes has a limit to its biological applicability domain which is at the cellular level. By changing the endpoint of the assay from subjective observation of beating cardiomyocyte inhibition to gene expression profiles of the differentiation route, it may become possible to improve the predictability of this assay. This will result in a better understanding of its biological applicability domain and therefore the knowledge of the limits and scope of its chemical applicability domain will increase. For example, the developmental toxicant 2,3,7,8-Tetrachlorodibenzo-p-dioxin (TCDD) is not detected by the cardiomyocyte readout of the ESTc, but can be predicted a developmental toxicant by an additional EST test in which osteogenesis is stimulated (de Jong et al., 2014). Defining the biological and chemical applicability domains of in vitro test systems such as the ESTc could facilitate test assay selection for chemical screening strategies or Integrated Approaches to Testing and Assessment (IATA) and with that improve toxicity predictions.

The ESTc has been shown to be an appropriate screening tool for triazoles by their interference with beating cardiomyocyte differentiation (Dimopoulou et al., 2018). Triazoles are designed to interfere with fungal cholesterol biosynthesis by inhibiting sterol 14α-demethylase cytochrome P450 (CYP51), which demethylates lanosterol in the cholesterol biosynthesis pathway (He et al., 2011, Pan et al., 2018). Like triazoles, morpholines and piperidines are classes of fungicides but which are less well studied including within the ESTc. They are also designed and shown to interfere with fungal cholesterol (=ergosterol) biosynthesis, but are structurally different from triazoles. Within the ergosterol biosynthesis pathway the morpholines and piperidines inhibit sterol Δ14-reductase and sterol Δ8,Δ7-isomerase which are important in the formation of 4,4-dimethylzymosterol or ergosterol, respectively (Fig. 1) (Pan et al., 2018, FRAC, 2021). Morpholines and piperidines can cause foetal malformations in rats, such as cleft palate formation after oral exposure to tridemorph at doses not toxic to the dams (Act, F.a.E.P., 1999).

Fig. 1

Ergosterol biosynthesis interference by azoles, morpholines and piperidines. Adapted and modified from (Pan et al., 2018). HMG-CoA = β-Hydroxy-β-methylglutaryl-CoA.

In order to apply a mechanism based ESTc readout, previous studies have involved gene transcript analysis using a hypothesis-driven targeted approach by preselecting gene transcript biomarkers based on existing literature (Mennen et al., 2019, Mennen et al., 2020, Mennen et al., 2021b, Mennen et al., 2021c). This approach has been successful for chemicals with known adverse effects for which it is possible to generate such hypotheses. However, unknown effects not supported by existing literature can be missed. Therefore, a comprehensive genome-wide expression screening could help improve our mechanistic understanding of chemical perturbations. Such an approach can be used to generate reasonable hypotheses by linking regulated pathways to phenotypic changes (Currie, 2012), and may ultimately improve hazard and risk assessment (Liu et al., 2019, Merrick, 2019). This may potentially avoid future confirmatory in vivo testing.

The objective of the present investigations was to derive an inclusive predictive biomarker profile for embryotoxicity, using a hypothesis-generating data-driven approach using genome-wide gene expression screening ‘NGS’, that would be able to eventually distinguish compounds within and between chemical classes. This approach was used to describe which differentiation routes appear during embryonic stem cell differentiation in the ESTc, comparing immature early stage (day four) and mature late stage (day ten) differentiation timepoints. The effects of the morpholines and piperidines on the ESTc biological domain were studied, by investigating differences in gene expression level changes between the structurally similar compounds within compound groups.

Methods

Stem cell culture

Murine embryonic stem cells (ES-D3 (D3), ATCC® (Manassas, VA, USA)) were maintained according to the previously described protocol (Mennen et al., 2019, Spielmann et al., 1997). The embryonic stem cells (ESCs) were maintained in 35 mm culture dishes (Corning, New York, NY, USA) in a humidified atmosphere at 37 °C with 5 % CO2 for stimulation of cell proliferation.. ESCs were replated in fresh medium every 2–3 days. The culture medium (CM) consisted of Dulbecco’s Modified Eagle’s Medium (DMEM; Gibco, Waltham, MA, USA), 20 % Foetal Bovine Serum (FBS; Greiner Bio-One, Kremsmünster, Austria); 2 mM l-Glutamine (Gibco); 1 % Non-Essential Amino Acids (NEAA; Gibco); 1 % 5000 IU/ml Penicillin/5000 µg/ml Streptomycin (Gibco); and 0.1 mM β-mercaptoethanol (Gibco). In order to preserve pluripotency, the ESCs in CM were supplemented with 1000 units/ml leukemia inhibitory factor (LIF; ESGRO®, Millipore, Burlington, MA, USA). These pluripotent ESCs were used in the differentiation assays.

Cell differentiation assay

A previously described protocol was used to differentiate ESCs into cardiomyocytes (Genschow et al., 2004, Spielmann et al., 1997). To enable differentiation, the CM as described for stem cell-culture was used without the addition of LIF. The differentiation protocol started with the hanging-drop method to form cell-aggregates called embryoid bodies (EBs) at differentiation day 0 (DD0). For this method, a 3.75·104 cells/ml cell suspension was added in droplets to the inside of the lid of a 100/20 mm CELLSTAR® cell culture dish (Greiner Bio-One). 5 ml of ice-cold phosphate buffered saline (PBS; Ca2+, Mg2+ free; Gibco) was added within the base of the culture dish and then lids were added after which the complete dishes were incubated for 3 days at 37 °C and 5 % CO2. At differentiation day 3 (DD3) exposure started by collecting the EBs in 5 ml of exposure medium containing test compound (see section test compounds) and were added to a 60 mm bacterial petri dish (Greiner Bio-One) to prevent attachment to the bottom of the dish. At differentiation day 5 (DD5), one EBs per well of a 24-wells plate (Greiner Bio-One) was transferred, each containing 1 ml of exposure medium. These EBs were cultured without further medium changes until differentiation day 10 (DD10). Four to five independent experiments were performed and EB samples were collected for gene expression analysis.

Test compounds

Two morpholines, tridemorph (TDM, CAS# 24602–86-6) and fenpropimorph (FPM, CAS# 67564–91-4), and two piperidines, fenpropidin (FPD, CAS# 67306–00-7) and spiroxamine (SPX, CAS# 118134–30-8) were tested. The triazole flusilazole (FLU, CAS# 85509–19-9), was included as a known positive control. The compounds were obtained from Sigma-Aldrich (Zwijndrecht, The Netherlands) and were tested in concentrations previously determined from concentration-response curves for inhibition in beating cardiomyocyte development (Mennen et al., 2021a). The tested ID10 (=inhibitory concentration at which 10 % beating inhibition occurs) and ID50 (=inhibitory concentration at which 50 % beating inhibition occurs) values are depicted in Table 1. All experimental conditions contained 0.25 % dimethyl sulfoxide (DMSO, CAS# 67–68-5, Sigma-Aldrich).

Table 1

ID10 and ID50 values of tested compounds (obtained from: (Mennen et al., 2021a).

	FLU	TDM	FPM	FPD	SPX
Differentiation ID10a	26 µM	57 µM	1.3 µM	0.1 µM	0.54 µM
Differentiation ID50b	42 µM	230 µM	5.5 µM	0.21 µM	1.2 µM

ID10: Concentration at which 10% beating inhibition occurs. bID50: Concentration at which 50% beating inhibition occurs.

Next Generation sequencing (NGS)

RNA collection and quality control

The samples comprised of collected EBs exposed to ID10 or ID50 compound concentrations or the vehicle control (0.25 % DMSO) at differentiation day 4 (DD4) and DD10. These days correspond to exposure periods of 24 h and 7 days in the ESTc, respectively. DD4 samples consisted of ∼ 56 EBs from one 60 mm plate per sample. The larger DD10 samples consisted of 24 EBs (one from each well) from one 24-wells plate per sample. Four to five samples per condition (five for all controls and the ID50 samples of DD10, four for all other samples) were collected from independent experiments and were transferred into Qiazol (Qiagen, Cat # 79306). The collected samples were stored at −80 °C prior to RNA isolation (RNeasy Mini-kit (Qiagen, Cat. # 74104) according to manufacturer’s protocol). Two additional steps were added to the RNA isolation protocol: a homogenisation step using QIAshredder columns (Qiagen, Cat. # 79654) and a DNase step using a RNase-Free DNase set (Qiagen, Cat # 79254). RNA quantity and quality were assessed using the Qubit3 (Invitrogen, Carlsbad, CA, USA) and the 2100 Bioanalyzer (Agilent Technologies, Amstelveen, The Netherlands).

RNA sequencing

The collected RNA samples were processed and sequenced by GenomeScan (Leiden, The Netherlands) using NGS which can sequence millions of fragments simultaneously per run and therefore can sequence hundreds to thousands of genes at one time. The NEBNext Ultra II Directional RNA Library Prep Kit for Illumina was used to process the samples. The sample preparation was performed according to the protocol “NEBNext Ultra II Directional RNA Library Prep Kit for Illumina” (NEB #E7760S/L). Briefly, mRNA was isolated from total RNA by polyA affinity purification using oligo-dT magnetic beads. After fragmentation of the mRNA, a cDNA synthesis was performed. This was used for ligation with the sequencing adapters and PCR amplification of the resulting product. The quality and yield after sample preparation was measured with the Fragment Analyzer. The size of the resulting products was consistent with the expected size distribution (a broad peak between 300 and 500 bp). Clustering and DNA sequencing using the NovaSeq6000 v1.5 was performed according to the manufacturer's protocols including a concentration of 1.1 nM of DNA, two flow cells, and NovaSeq control software NCS v1.7. Image analysis, base calling, and quality check was performed with the Illumina data analysis pipeline RTA3.4.4 and Bcl2fastq v2.20. Expression levels of the transcripts (20 million paired end reads per sample) were quantified against the mouse reference genome (Ensembl GRCm38.p6, containing 22.519 coding genes) using TopHat version 2.0.14.

RNA sequencing analysis

The obtained count-tables were used for further differential gene expression analysis using Rstudio statistical software (version 4.1.0). Control samples for DD4 and DD10 were compared by differential gene expression using DESeq2 (version 1.30.0), adjusting for the sampling day per individual independent experiment to extract log2FC values of genes with at least one count in the analysis (=20.335). Differentially Expressed Genes (DEGs) were obtained by filtering results for p < 0.001 and log2FC > 1.5. Functional enrichment analysis of genes in Gene Ontologies (GO) terms (biological processes GO BP5) was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/ consulted in September 2021, version 6.8) as a gene list with genes of at least one count in the analysis. GO-terms were clustered using the ‘functional annotation clustering’ tool and summary names (see supplementary material) were obtained with help of http://amigo/geneontology.org/amigo based on their tree-view.

The compound treated samples were analysed in a similar manner. First, for visualisation, principal component analysis (PCA) was performed on count data after filtering out the 0 values and applying a Variance Stabilizing Transformation (VST). Additionally, data was corrected for differences between experiments. Determination of DEGs was performed using DESeq2 after filtering results for p < 0.001 and log2FC > 0.5 and the number of up and downregulated genes per condition were determined compared to the control (DMSO) values. PCA plots and heatmaps were generated in RStudio (version 4.0.0), Venn diagrams were assembled using Venny (version 2.1.0), and the remaining graphs were visualised using GraphPad Prism (version 8.1.2, https://www.graphpad.com).

Results

Gene expression changes in unexposed control cultures

To investigate cell types developing within the ESTc, we compared DD4 and DD10 gene expression in controls relative to each other by DESeq analysis. Results showed that 1255 / 20,335 genes were upregulated on DD4 and 2987 / 20,335 genes were upregulated on DD10 (Fig. 2A). These DEGs were organised into GO-clusters and the enrichment scores were given (Fig. 2B, 3C). Six GO-term clusters were regulated on DD4 were mainly related to general cell processes and mechanisms. These clusters included RNA metabolic process, small molecule metabolic process, primary metabolic process, meiotic cell cycle process, embryonic morphogenesis, and synaptic signalling (Fig. 2B). At DD10, cell differentiation related GO-term clusters were significantly higher regulated as compared to DD4. The top five differentiation related clusters, according to enrichment score, were circulatory system development, (skeletal) system development, heart development, muscle and organ tissue development, and nervous system and cell development. In line with the experimental purpose of understanding the differentiation occurring in this test system, these top five clusters (in blue, Fig. 2C) were used for the selection of related GO-terms which were examined on DD10 for their perturbation by the test compounds.

Fig. 2

Genome-wide specific gene expression changes in controls per differentiation day. A) The number of Differentially Expressed Genes (DEGs) after comparison between day 4 and day 10 as part of the total number of assessed genes. B) Enrichment scores per enriched cluster of GO-terms related to the 1255 genes specific to day 4, or C) 2987 genes specific to day 10. The top 5 enriched clusters related to specific differentiation routes are depicted in blue and were selected for further analysis. Enrichment score = -log(p-value) for which the p-value is the mean/median of GO-terms belonging to each cluster. P < 0.001, log2FC > 1.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Gene expression changes after compound exposures

To study the extent that the morpholines and piperidines perturbed normal cell differentiation in the ESTc, effects of exposure to ID10 and ID50 concentrations of test compounds were assessed on DD4 or DD10, corresponding to 24 hrs or 7 days exposure duration, respectively. Using a PCA plot of DD4 data, individual exposure and control samples per experiment showed a scattered pattern with low PC1 and PC2 values of 14.9 % and 12 %, respectively (Fig. 3A). DEGs count on DD4 was highest after exposure to FLU, resulting in > 20 regulated DEGs per exposure condition (Fig. 3B). On DD10, the PCA plot showed more clear distinctions between samples exposed to the vehicle, FLU or the test compounds with higher PC values: PC1 53.6 % and PC2 8.9 % (Fig. 3C). There was a clear difference in the PC2 related direction of FLU versus the PC1 related shift of the other test compounds relative to the position of the vehicle control. As anticipated, concentration responses were apparent for all compounds with ID10 samples closer to the vehicle controls as compared to the ID50 sample responses. The TDM ID50 showed the largest distance from the vehicle controls, followed by FPD. Also, TDM ID10 levels showed a relatively large distance from the vehicle controls as compared to the ID10 concentrations of the other test compounds. These differences between test groups in the PCA plots were also reflected in the DESeq analysis and DEGs showed most regulated genes in the samples exposed with TDM ID50 or FPD ID50 (Fig. 3D).

Fig. 3

Genome-wide gene expression regulations in exposure groups on days 4 and 10. A) Principal Component Analysis (PCA) of all experimental groups on day 4. B) Number of DEGs per exposure group with P < 0.001 and log2FC > 0.5; solid coloured bars indicating upregulated genes and the white bars indicate downregulated genes. C) PCA and D) DEGs per exposure group on day 10. ID10: Concentration at which 10 % beating inhibition occurs. ID50: Concentration at which 50 % beating inhibition occurs. Positive control - FLU: flusilazole. Morpholines - FPM: fenpropimorph, TDM: tridemorph. Piperidines – FPD: fenpropidin, SPX: spiroxamine.

Differentiation day 4 gene expression analysis

The expression of regulated genes on DD4 was compared between exposures using Venny and a heatmap (Fig. 4). Commonly regulated genes by at least two compounds at ID10 (Fig. 4A) and ID50 levels (Fig. 4B) were extracted from the Venny diagrams. The commonly regulated genes exhibited some overlap between compounds at ID10 and ID50 levels and resulted in the 6 genes that were regulated by two or more compounds: methylsterol monooxygenase 1 (Msmo1), actin alpha cardiac muscle 1 (Actc1), and α-smooth muscle actin (Acta2), structural maintenance of chromosome 1b (Smc1b), transforming growth factor beta receptor 3 (Tgfbr3), and solute carrier family 38 member 3 (Slc38a3). All conditions resulted in upregulations of these genes with log2FCs up to 1.1 (Fig. 4C). The conditions clustered pairwise per compound, except for the SPX conditions. Smc1b showed most diversity in expression levels among the test conditions. The uniquely regulated genes in Fig. 4A and 4B were listed in supplementary table 1.

Fig. 4

Common and unique DEGs per test group at day 4. A) Venny diagram for ID10 concentrations and B) ID50 concentrations. C) Heatmap of genes commonly regulated by at least two test compounds. Colour key indicates the log2FC. ID10: Concentration at which 10% beating inhibition occurs. ID50: Concentration at which 50% beating inhibition occurs. Morpholines - FPM: fenpropimorph, TDM: tridemorph. Piperidines – FPD: fenpropidin, SPX: spiroxamine.

Differentiation day 10 gene expression analysis

Test conditions were examined for their interferences with the top five differentiation routes regulated in control cultures, based on GO-terms and on individual gene expression related to these GO-terms.

Analysis of effects on GO-terms

Within each of the five GO-terms, the enrichment value and number of regulated genes by each condition were assessed (Fig. 5). ID10 and ID50 showed clear concentration-responses for all compounds and GO-terms analysed. The enrichment values of GO-terms differed in magnitude between test compounds. TDM and FPD ID50 conditions showed highest enrichment for the GO-term circulatory system development. FLU conditions showed highest enrichment for skeletal system development and nervous system development. Smaller differences in enrichment values were found for the GO-terms heart development and muscle organ development. Although the ESTc was designed to detect cardiomyocyte differentiation effects, the largest effects of the test compounds were found on nervous system development.

Fig. 5

GO-term enrichment per exposure group. The -log(p-value) displays the significance of enrichment. Circle size indicates the number of regulated genes (as indicated) by each condition. The red dotted line indicates a P-value of 0.001. All values to the right of this line are P < 0.001. ID10: Concentration at which 10 % beating inhibition occurs. ID50: Concentration at which 50 % beating inhibition occurs. Positive control - FLU: flusilazole. Morpholines - FPM: fenpropimorph, TDM: tridemorph. Piperidines – FPD: fenpropidin, SPX: spiroxamine. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Analysis of effects on individual gene expression

Shared and unique genes per GO-term regulated between test compounds were visualised in the Venn-diagrams (Fig. 6). To identify the exclusivity of the GO-terms by means of unique genes per differentiation route and to rule out potential overlap of DEGs between GO-terms, the common DEGs per GO-term identified in Fig. 6 were visualised in Fig. 7. The genes responsive in the GO-term heart development were all shared by the GO-term circulatory system development. The GO-terms circulatory system development, skeletal system development, muscle organ development, and nervous system development also showed overlap between commonly regulated DEGs, but also contained uniquely regulated genes. Chemical regulation of the common DEGs per GO-term obtained from Fig. 6 were evaluated and visualised in heatmaps (supplementary fig. S1-S5). The assessed genes per GO-term are listed in supplementary table 3. Generally, the heatmaps separated FLU and ID10 conditions of FPM, SPX, and FPD from the remaining ID50 test compounds and TDM ID10. For all the test compounds except FLU, genes were always regulated in the same direction of up- or down-regulation. Therefore, no distinction as to specific gene up- or down-regulation could be made between the morpholine and piperidine group, nor within these structural groups. Xin Actin Binding Repeat Containing 2 (Xirp2) was most upregulated after chemical exposure compared to controls within the GO-term circulatory system development, whereas C-X3-C Motif Chemokine Receptor 1 (Cx3cr1) and C—C Motif Chemokine Receptor 5 (Ccr5) were the most downregulated genes (fig S1). The latter two genes were also most downregulated within the GO-term nervous system development (fig S5). Additionally, G-protein coupled receptor 183 (Gpr183) was one of the genes that was most downregulated within the nervous system development GO-term (fig S5). Lumican (Lum) involved in collagen fibril organisation was downregulated within the GO-term skeletal system development, especially after exposure to TDM (fig S2). Xin Actin Binding Repeat Containing 2 (Xirp2) was most upregulated within both GO-terms heart development and muscle organ development (fig. S3 and S4). In some cases FLU regulated genes in an opposite direction (fig S1-S5, Table 2). These genes were related to multiple GO-terms and were not unique to one of the GO-terms.

Fig. 6

Venn-diagrams of common and specific Differentially Expressed Genes (DEGs) per GO-term of ID50 test groups at day 10. ID10: Concentration at which 10% beating inhibition occurs. ID50: Concentration at which 50% beating inhibition occurs. Morpholines - FPM: fenpropimorph, TDM: tridemorph. Piperidines – FPD: fenpropidin, SPX: spiroxamine.

Fig. 7

Venn-diagram of common chemical regulated DEGs per development related GO-term of ID50 test groups at day 10.

Table 2

List of genes that were upregulated by FLU, but downregulated by the other test compounds relative to the DMSO control.

	Circulatory system development	Skeletal system development	Heart development	Muscle organ development	Nervous system development
Epha3Ephrin type-A receptor 3	X		X		X
Igf1insulin growth factor 1	X	X	X	X	X
Kcnab1potassium voltage-gated channel subfamily A regulatory beta subunit 1	X		X	X	X
Tcf21transcription factor 21	X
Sfrp1Secreted frizzled related protein 1		X			X
Col13a1collagen type XIII alpha 1 chain		X
HlxH2.0 like homeobox				X	X

Crosses indicate the genes and the GO-terms in which they appear.

Uniquely regulated genes by the test compounds within ID50 levels are listed per GO-term in supplementary table 2. TDM generally showed the most uniquely regulated genes, whereas FPM showed the least uniquely regulated genes. Also, in this case, the listed genes were not always unique to each GO-term.

Discussion

This study has advanced the understanding of the biological domain within the ESTc. This has been possible through mapping the biological domain in terms of emerging cell differentiation routes within the ESTc by applying a hypothesis-generating data-driven approach using a genome-wide gene expression screen (NGS). This tool has the advantage to sequence millions of fragments simultaneously per run and therefore can sequence hundreds to thousands of genes at one time. A comparison between an early and a late time-point within the ESTc protocol confirmed cardiomyocyte differentiation and also revealed the presence of additional differentiation routes in terms of developmental systems or even cell type specific GO-term analysis. These GO-terms were regulated by all compounds and thus gave mechanistic insight into the perturbations by morpholines and piperidines on multiple stem cell differentiation routes.

Compared to existing literature describing the use of transcriptomics with the ESTc (Pennings et al., 2011, Schulpen et al., 2014, van Dartel et al., 2009a, van Dartel et al., 2010, van Dartel and Piersma, 2011, van Dartel et al., 2014), the current study used a more comprehensive approach for gene expression analysis as there was no ‘a priori’ selection of a subset of genes to be measured and instead the expression level of all expressed genes (>20 k) was determined and then organised into GO-terms. This approach revealed that the GO-cluster for nervous system development in particular represented the highest number of regulated genes that were unique to this GO-term as illustrated in Fig. 7. The presence of a sub-population of neuronal cells within the ESTc has been previously investigated. Using a different study approach, the presence of tubulin beta 3 class III (TUBB3) positive cells, indicative of neuronal differentiation, has been shown in the ESTc (Mennen et al., 2021a). Also, neural crest cells have been observed within the ESTc (Mennen et al., 2019, Mennen et al., 2021b), but this differentiation route is relatively less represented in the current study (Fig. 2; supplementary material). Within the employed ESTc protocol it was previously shown that embryoid body formation (day 0–3) is driven by cell proliferation and cardiomyocyte differentiation starts after plating the embryoid bodies on tissue culture plates from day 3 onwards (van Dartel et al., 2009b). Therefore it was of interest to start exposure from day 3. However, early commitment of ESCs to lineage fates also occurs before day 3 after the omission of LIF in the culture medium (Fehling et al., 2003, Anton et al., 2007). This is exemplified by the expression of Brachyury as a marker for mesoderm commitment by differentiation day 3 (van Dartel et al., 2009b). Therefore effects of compounds on early ESC commitment are not measured using this protocol, which focuses on subsequent differentiation from committed lineage fates.

On DD4, genes specific to basic cell processes and mechanisms were mainly expressed and chemical treatments showed regulation of six common genes between the test compounds. These genes had mixed functions with half of them being related to the intended mode of action of the test compounds or to the original readout of the ESTc. These genes concerned methylsterol monooxygenase 1 (Msmo1), actin alpha cardiac muscle 1 (Actc1), and α-smooth muscle actin (Acta2). Interestingly, Msmo1 has previously been found to be commonly regulated by the test compounds within the ESTc (Mennen et al., 2021a). The selection of this gene was based on existing literature, whereas in this genome-wide analysis it was one of the significant genome-wide regulated genes on DD4. FLU also induces Msmo1 in the rat whole embryo culture (Robinson et al., 2012). The test compounds affect the sterol biosynthesis pathway by accumulation of sterol intermediates upon FPM and SPX exposure in hiPSCs (Wages et al., 2020) and by inhibition of Δ8-Δ7 isomerase by exposure to FPM in rat liver homogenates (Ruan et al., 2000). Our results confirm Msmo1 is a useful gene for early developmental toxicity screening for compounds with similar modes of action. The genes Actc1 and Acta2 are both markers for early heart development and have been studied in relation to perturbations of murine embryonic stem cell differentiation before (KalantarMotamedi et al., 2016, Potta et al., 2010, Smirnova et al., 2014). The other three genes, structural maintenance of chromosome 1b (Smc1b), transforming growth factor beta receptor 3 (Tgfbr3), and solute carrier family 38 member 3 (Slc38a3), were studied in cell differentiation in different contexts. Smc1b has been studied in relation to the formation of germ cells from embryonic stem cells and meiosis and was influenced by BMP4 (bone morphogenetic protein 4) (Esfandiari et al., 2017). Tgfbr3 has been studied in relation to osteogenic differentiation (Atala, 2020, Cook et al., 2019, Shojaee et al., 2019), but also in relation to hypoxia, epithelial to mesenchymal transition, liver development, heart morphogenesis and muscular septum morphogenesis (www.genecards.org; assessed April 2022). Lastly, the amino acid transporter Slc38a3 has been studied in relation to lung cancer development (Wang et al., 2017, Person et al., 2015). Although the differences in regulation of these six genes by the test compounds were small, especially the genes Msmo1, Actc1, and Acta2 related to cardiac muscle cell differentiation and chemical mode of action (MOA) may be relevant biomarker genes when studying differentiation perturbations on DD4.

On DD10, the test compounds also regulated other differentiation routes in addition to the original readout of the ESTc. The morpholines and piperidines have been tested before in the ESTc in relation to neuron differentiation and regulated gene expression of Tubb3 (Mennen et al., 2021a). The effects of FPM and SPX have also been tested in several human cell lines and developing human neural tissue (Wages et al., 2020). In these cells, FPM and SPX increased levels of 7- and 8-dehydrocholesterol and reduced levels of desmosterol and cholesterol. Therefore, the authors concluded these compounds may be developmental neurotoxicants as cholesterol is an essential lipid in the central nervous system and its metabolism is affected in many neurodevelopmental disorders (Wages et al., 2020). FPM has also been assessed with regard to craniofacial malformations in zebrafish embryos where a set of marker genes were selected that gave insight into the mode of action associated with skeletal mal-development (Heusinkveld et al., 2020).

At equipotent concentrations, as measured by inhibition of cardiac differentiation, TDM and FPD regulated ∼ 5000 DEGs while the other compounds regulated ∼ 2000 DEGs (Fig. 2D). Large differences in gene expression regulation by equipotent concentrations of compounds have been observed before, e.g. in ESTn, in which carbamazepine regulated far fewer genes as compared to valproic acid, which are both anti-epileptic drugs (Schulpen et al., 2015). This indicates that gene expression analysis offers a very different perspective on compound effects and potency, in particular providing additional information on molecular regulation, that can inform about mechanism of action.

The test compounds regulated genes within all prioritized GO-terms. The commonly affected genes were all regulated in the same direction when comparing test compounds and therefore discriminate between compounds. However, uniquely regulated genes were also found for each test compound within each GO-term that may help in the identification of compound specific effects when comparing the morpholines and piperidines (supplementary table 2). These dozens of genes should be verified for their uniqueness and functional properties in additional experiments. These findings are indicative of common and unique mechanisms of toxicity induced by the selection of the test compounds in these experiments. The presence of such potential differences in mechanisms were not observed or studied in in vivo studies (Act, F.a.E.P., 1999, EFSA, 2008, Pfeil, 2004, Report, 2007, ECHA, 2015, EFSA, 2017, Adcock et al., 2007). Apart from mechanistic comparisons of dynamic effects, differences in potency, kinetics and metabolism in vivo may affect embryotoxicity. Studies in adult rats show clearance of FPM, TDM, and FLU from the body into urine and faeces (Adcock et al., 2007, EC, Directive, 2009, Hawkins et al., 1974), while information on kinetics in pregnancy is often not available for non-pharmaceutical chemicals. This hampers the estimation of relative potency.

The expression changes of commonly affected genes occurred in some cases in the opposite direction for the test compounds as compared to FLU, which is indicative of the perturbations of different mechanisms. Seven genes were differently regulated: Ephrin type-A receptor 3 (Epha3), insulin growth factor 1 (Igf1), potassium voltage-gated channel subfamily A regulatory beta subunit 1 (Kcnab1), transcription factor 21 (Tcf21), Secreted frizzled related protein 1 (Sfrp1), collagen type XIII alpha 1 chain (Col13a1), and H2.0 like homeobox (Hlx). These genes were all upregulated by FLU and were not unique to the related GO-terms (Table 2), but are involved in embryo development. Genes Igf1, Epha3, and Tcf21 all play a role in heart development. Igf1 is involved in expanding the developing mesoderm and promoting cardiac differentiation (Engels et al., 2014), which would hold true for FLU since it upregulated Igf1 expression, but not for the test compounds. A repression of Igf1 by bisphenol A in human ESC (hESC) EBs is correlated with a decreased neural cell differentiation (Huang et al., 2017). hESC maintenance by the addition of IGF to Activin containing medium supported pluripotency through PI3K/mTOR signalling (Wamaitha et al., 2020). Epha3 is a receptor kinase necessary for the fusion of the ventricular septum and atrioventricular cushions during heart development (Dilg et al., 2016). Tcf21 is an epicardial marker in heart development and a progenitor of ventricular cardiomyocyte and pharyngeal muscle (Lupu et al., 2020, Braitsch et al., 2013, Dohn et al., 2019). Tcf21 is also involved in the mesenchyme of developing organs like the kidney, lung and gut (Cui et al., 2003). Col13a1 and Hlx are also involved in mesenchyme cells and organ development. Col13a1 is a collagen involved in the mesenchymal subtype in the lungs and causes congenital myasthenic syndrome type 19 (Yuan et al., 2020, Logan et al., 2015). This collagen is also important in the basal lamina of neuromuscular junctions and mice lacking Col13a1 show immature nerve terminals and reduced neurotransmission (Maselli et al., 2012). Hlx enhanced the appearance of premature reprogramming cells in hiPSCs and interfered with pluripotency (Yamakawa et al., 2016). During embryogenesis, Hlx is prominently expressed in visceral mesenchyme of the developing liver, gall bladder and gut (Hentsch et al., 1996). Sfrp1 is the counter-acting molecule of Wnt and seems to have a role in rostral- and caudal regulation of ESC-derived neuroectoderm (Takata et al., 2017), and in differentiation of stem cells to dopaminergic neurons (Kwon et al., 2014). Kncab1 has not been studied in relation to stem cells or embryo development, except for its association with elevated birth weight, which may have been attributed to gestational duration (Beaumont et al., 2018). In summary, these differently regulated genes between FLU and the test compounds are involved in multiple differentiation routes.

Overall, this hypothesis-generating data-driven approach provided a valuable and additional perspective on the biological domain of the ESTc, through the novel mechanistic information from the large quantity of gene expression analysis compared to methods with a priori gene selection. Given the conserved nature of the developmental mechanisms of vertebrate embryonic cell differentiation represented in the mouse-derived ESTc, these mechanisms are likely to be relevant for human safety prediction and protection as well. The overlapping and unique gene regulations of the tested compounds advances our knowledge of the chemical applicability domain of the ESTc. This progressive understanding and knowledge of both the biological and chemical applicability domains for this assay could contribute to future toxicity predictions by facilitating selection of reliable and relevant test assays for effective chemical screening strategies, themselves part of Integrated Approaches to Testing and Assessment (IATA). This refined ESTc method, together with other in vitro and in silico test systems in combination with kinetic modelling, can be instrumental for the contextualisation of such test batteries for improved prediction and protection of human development, enabling hazard and risk assessment with reduced dependence on in vivo animal models.

CRediT authorship contribution statement

R.H. Mennen: Investigation, Formal analysis, Visualization. N. Hallmark: Supervision. M. Pallardy: Supervision. R. Bars: Supervision, Funding acquisition. H. Tinwell: Supervision. A.H. Piersma: Supervision, Conceptualization.

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.