An Integrative Multi-Omics Workflow to Address Multifactorial Toxicology Experiments

Abstract

Toxicology studies can take advantage of omics approaches to better understand the phenomena underlying the phenotypic alterations induced by different types of exposure to certain toxicants. Nevertheless, in order to analyse the data generated from multifactorial omics studies, dedicated data analysis tools are needed. In this work, we propose a new workflow comprising both factor deconvolution and data integration from multiple analytical platforms. As a case study, 3D neural cell cultures were exposed to trimethyltin (TMT) and the relevance of the culture maturation state, the exposure duration, as well as the TMT concentration were simultaneously studied using a metabolomic approach combining four complementary analytical techniques (reversed-phase LC and hydrophilic interaction LC, hyphenated to mass spectrometry in positive and negative ionization modes). The ANOVA multiblock OPLS (AMOPLS) method allowed us to decompose and quantify the contribution of the different experimental factors on the outcome of the TMT exposure. Results showed that the most important contribution to the overall metabolic variability came from the maturation state and treatment duration. Even though the contribution of TMT effects represented the smallest observed modulation among the three factors, it was highly statistically significant. The MetaCore™ pathway analysis tool revealed TMT-induced alterations in biosynthetic pathways and in neuronal differentiation and signaling processes, with a predominant deleterious effect on GABAergic and glutamatergic neurons. This was confirmed by combining proteomic data, increasing the confidence on the mechanistic understanding of such a toxicant exposure.

1. Introduction

Metabolomics is now a well-established component in the set of omics approaches applied to describe processes of neurological diseases. It constitutes an accepted evolution upon the development of genomics, epigenomics, transcriptomics and proteomics [1]. Being the final step in the cascade of events from genes to living organisms, metabolomics provides a holistic view of the current state of a system as a result of the interaction between the environment and its own genetic, transcript, and protein profiles and perturbations. Recently, efforts were made [2,3,4] to harmonize the attempts to include omics in the field of regulatory toxicology, highlighting the relevance of systems biology approaches to support decision-making based on a better mechanistic understanding of chemical hazards. Untargeted omics workflows are intrinsically well-suited to the collection of biological information without the need for previous knowledge, thus potentially unveiling unforeseen changes in the systems under study. Indeed, they have the potential to become a valuable tool when it comes to the exploration of new mechanisms of toxicity [5,6]. Proteins need a longer period to be either synthetized or degraded [7,8], while metabolomics can give a picture of both long-term and instantaneous reactions of the system induced by the stressor under study [9,10,11].

Toxicology studies aiming at revealing mechanisms of action of different substances often require a multitude of investigative techniques. Therefore, there is an increasing need to not only understand the mechanisms through which these compounds may affect biological systems, but also to identify key events shared among different toxicity pathways that can be used for faster assessment of potential hazard of compounds. Such key events are at the core of the Organization for Economic Cooperation and Development (OECD) framework Adverse Outcome Pathways (AOPs), and aim at bringing mechanistic toxicological data into a context relevant for risk assessment application [12]. This concept adds another challenge to the analytical strategy, requiring approaches making possible to deal with multifactorial experimental designs while enabling a joint biological interpretation of data issued from different analytical platforms.

Data generated from untargeted omics analyses typically contains information comprising thousands of variables. In order to summarize, explore and discover information within such large number of signals, dimensionality reduction methods have been widely adopted. Both unsupervised and supervised methods are employed for this purpose, with supervised ones presenting the advantage of being able to take into account the experimental design (DOE) in the process of maximizing the variance between samples coming from different experimental groups [13]. Nevertheless, classical dimensionality reduction methods fall short when several experimental factors are involved in the design, and although a few alternatives have been proposed (ASCA, ANOVA-PCA, etc.) to circumvent this limitation, they are not able to provide a single model for interpretation. Recently, ANOVA multiblock OPLS (AMOPLS) was proposed as an appropriate tool, allowing the study of complex omics datasets by taking into account the DOE information and offering a single model to understand the behavior of the system [14]. The model allows the assessment of the cumulative contribution of each metabolite according to its explained variance. Therefore, it is specifically possible to relate each compound to the way each different experimental factor has been exerted its influence. In addition, AMOPLS can be joined to a consensus approach, thus making possible to combine multiple tables obtained from separate analytical methods, and to evaluate the contribution of each technique to the variability observed in the model. In this study, consensus AMOPLS was used to combine data obtained from the different LC-MS platforms in metabolomics (Figure S1).

The present work focuses on a strategy comprising factor deconvolution and multi-omics integration to maximize information retrieval from toxicology studies. As an example, phenotypic changes taking place in primary 3D rat neural cell cultures upon exposure to the heavy metal trimethyltin (TMT) were studied. TMT has been shown to induce neurotoxicity and neuroinflammation in numerous studies, both in vitro and in vivo [15,16,17,18]. Taking previous in-house experience with this substance in the 3D rat model [19,20], it was considered an appropriate model substance for studying neuroinflammatory and neurotoxic effects. Three different experimental factors were considered: Maturation state of the cell culture, exposure duration and concentration of toxicant. A multiplatform metabolomic strategy was applied for the analysis of the cell contents, comprising four different and complementary LC-MS setups followed by a confirmatory analysis in proteomics. The contribution of each experimental factor was decomposed and independently quantified by means of the AMOPLS multivariate analysis and the information arising from metabolomic and proteomic analyses combined with a pathway enrichment tool. The merged omics signatures were thus used to identify putative toxicity pathways connecting the TMT exposure to the cellular response.

2. Materials and Methods

2.1. Solvents and Additives

UPLC-MS-grade methanol, acetonitrile and water were obtained from Fisher (Fisher Scientific, Loughborough, UK). Ammonium formate and ammonium hydroxide for MS were supplied by Sigma-Aldrich (Buchs, Switzerland). UPLC-MS-grade formic acid was obtained from Biosolve (Valkenswaard, The Netherlands).

2.2. Culture of 3D Rat Brain Cell Cultures and TMT Exposure

Serum-free 3D brain cell cultures were prepared from the whole brain of 16-day embryonic rats (Sprague Dawley, Janvier, France) as previously described [21]. Neural cells were seeded at a density of 6 × 10⁷ cells/flask and spontaneously aggregate into 3D spheres. The 3D cultures were maintained in a chemically defined medium under constant gyratory agitation (80 rpm) at 37 °C in an atmosphere of 10% CO₂ and 90% humidified air. Media were freshly replenished every third day until day in vitro (DIV) 14, and every second day thereafter until the end of the in vitro culturing period. For experimentation, replicate cultures were prepared by randomizing and aliquoting the free-floating spheres of the original cultures. The day of culture preparation is designated DIV 0. Cultures were exposed at an immature stage of differentiation (from DIV 5) for 24 h (to DIV 6) or 10 days (to DIV 15), or at a matured stage of differentiation, starting at DIV 25 for 24 h (to DIV 26) or 10 days (to DIV 35) to 0.0 (control), 0.5 or 1.0 µM TMT. TMT concentrations were maintained constant with supplementation of TMT with each change of culture media.

2.3. Sample Preparation for Metabolomic Analyses

Metabolomic analyses were performed on two different LC modes (reversed-phase LC (RPLC) and hydrophilic interaction LC (HILIC)) and on positive and negative electrospray ionization polarities. Triplicate samples were obtained from every combination of the studied experimental factors: two maturation states (immature/mature), two exposure lengths (24 h or 10 days), and three TMT concentrations (0.0, 0.5, and 1.0 µM). To stop the enzymatic activity and to perform the extraction of intracellular metabolites, 100 µL of ice-cold MeOH:H₂O 90:10 v/v were added to each sample of neurospheres after removal of the culture medium. Mechanical disruption was achieved by ultrasound treatment on ice (2 × 5 s pulses, 30% amplitude, 2 mm-diameter probes operated on a VibraCell VCX 500 sonicator, Sonics & Materials Inc., Newtown, CT, USA). Samples were then centrifuged (12000× g, 15 min, 4 °C) and supernatants were collected. Equal aliquots of all cell extracts were pooled to generate a post-extraction quality control (QC) sample used for analytical performance evaluation and data treatment. Extracts were diluted 10 times prior to LC-MS analyses using either H₂O:MeCN 95:5 v/v (RPLC) or MeCN:H₂O 95: 5 v/v (HILIC).

2.4. LC-MS Analyses for Metabolomics

Chromatography was performed on a Waters H-Class Acquity UPLC system composed of a quaternary pump, a column manager and an FTN autosampler (Waters Corporation, Milford, MA, USA). For RPLC analyses, samples were separated on a Kinetex C18 column (150 × 2.1 mm, 1.7 μm) and the corresponding SecurityGuard Ultra precolumn and holder (Phenomenex, Torrance, USA). Solvent A was H₂O and solvent B was MeCN, both containing 0.1% formic acid. The column temperature and flow rate were set at 30 °C and 300 µL min⁻¹, respectively. The gradient elution was as follows: 2% to 100% B in 14 min, hold for 3 min, then back to 2% B in 0.1 min and re-equilibration of the column for 7.9 min. HILIC separations were conducted on a Waters Acquity BEH Amide column (150 × 2.1 mm, 1.7 µm) bearing an adequate VanGuard pre-column. Solvent A was H₂O:MeCN (5:95, v/v) and solvent B was H₂O:MeCN (70:30, v/v) containing 10 mM ammonium formate (pH = 6.5 in the aqueous component). The following gradient was applied: 0% B for 2 min, increased to 70% B over 18 min, held for 3 min, and then returned to 0% B in 1 min and to re-equilibrate the column for 7 min (total run time was 31 min). The flow rate was 500 µL min⁻¹, and the column temperature was kept at 40 °C.

In all cases, a sample volume of 5 µL was injected. Samples were randomized for injection, and QC pools were analysed every six samples to monitor the performance of the analytical platform.

The UPLC system was coupled to a maXis 3G Q-TOF high-resolution mass spectrometer from Bruker (Bruker Daltonik GmbH, Bremen, Germany) through an electrospray interface (ESI). The instrument was operated in TOF mode (no fragmentation). The capillary voltage was set at −4.7 kV for ESI+, the drying gas temperature was 225 °C, the drying gas flow rate was set at 5.50 (RPLC) or 8.00 (HILIC) L min⁻¹ and the nebulizing gas pressure was 1.8 (RPLC) or 2.0 bar (HILIC). Transfer time was set at 40 (RPLC) or 60 (HILIC) µs and pre-pulse storage duration at 7.0 (RP) or 5.0 μs (HILIC). For ESI– operation, the capillary voltage was set at 2.8 kV. All the remaining ion source and ion optics parameters remained as in ESI+. Data between 50 and 1000 m/z were acquired in profile mode at a rate of 2 Hz. ESI and MS parameters were optimised using a mix of representative standards fed by a syringe pump and mixed with the LC eluent (mid-gradient conditions) within a tee-junction. Formate adducts in the 90–1247 m/z range were employed for in-run automatic calibration using the quadratic plus high-precision calibration algorithm provided by the instrument’s manufacturer. MS and UPLC control and data acquisition were performed through the HyStar v3.2 SR2 software (Bruker Daltonik) running the Waters Acquity UPLC v.1.5 plug-in.

2.5. Metabolomics Data Pretreatment and Metabolite Identification

Run alignment, peak piking and sample normalization were performed on Progenesis QI v2.3 (Nonlinear Dynamics, Waters, Newcastle upon Tyne, UK) and peaks were identified by matching their retention times, accurate masses and isotopic patterns to those of a library of chemical standards (MSMLS, Sigma-Aldrich, Buchs, Switzerland) analyzed under the same experimental conditions, as described elsewhere [22]. Peaks were evaluated on the basis of the repeatability of their areas on the QC samples and discarded when the coefficient of variation was above 35% within the sequence.

2.6. Sample Preparation for Proteomic Analysis

Duplicate samples of control and 1.0 μM TMT-treated cell pellets (immature cultures, exposed to the toxicant for 10 d) were solubilized in 6 M urea, 1 mM DTE and 0.1 M TEAB. After homogenization steps, solubilized proteins were recovered by centrifugation (supernatants). Reduction of 10 µg of protein per sample was done with TCEP (tris (2-carboxyethyl)phosphine), final concentration of 10 mM, and samples reacted for 30 min at 30 °C. Alkylation was done with iodoacetamide, added to a final concentration of 40 mM, and samples incubated 60 min in the dark at room temperature. Trypsin was added (ratio of 1:25, w/w), and the digestion was performed overnight at 37 °C. Then, a 10-plex tandem mass tags (Thermo Scientific, Rockford, USA) labeling was performed according to manufacturer’s instructions. Tandem mass tags reagent was dissolved in MeCN, and each sample was incubated 60 min at room temperature with a specific tag. For tandem mass tags quenching, 8 µL of hydroxylamine 5% (v/v) were added and incubated with the samples for 15 min. Tandem mass tags-labeled samples were pooled and dried under vacuum. Samples were dissolved in 5% MeCN/0.1% formic acid and desalted with C18 micro-spin columns (Thermo Scientific). Peptides were separated by off-gel electrophoresis, desalted and solubilized in an appropriate amount of 5% MeCN/0.1% formic acid for mass spectrometry analysis.

2.7. MS-Data Acquisition for Proteomics Analyses

For LC-MS/MS, peptides were dissolved in 5% MeCN/0.1% formic acid to a concentration of 0.25 µg/µL. Mass spectrometry experiments were performed on a Q Exactive Plus instrument (Thermo Scientific, San Jose, CA, USA) equipped with an Easy-nanoLC from Thermo. Peptides were trapped on 2 cm × 75 µm ID, 3 µm pre-column and separated on an Easy-spray column, 50 cm × 75 µm ID, PepMap C18, 2 µm (Thermo Scientific). The analytical separation was run for 60 min using a gradient of H₂O (solvent A) and MeCN (solvent B), both containing 0.1% formic acid, at a flow rate of 300 nL min⁻¹. For MS survey scans, the resolution was set to 140,000 and the ion population was set to 3 × 10⁶ with an m/z window from 350 to 2000. Twenty precursor ions were selected for higher energy collisional dissociation (HCD) with a resolution of 35,000, an ion population set to 1 × 10⁵ (isolation window of 0.5 m/z) and a normalized collision energy set to 30%.

2.8. Protein Identification and Quantification

Raw data were loaded on Proteome Discoverer 2.2 software for identification and/or quantification of peptides and proteins. Identification was performed in the UniProt/SwissProt rat database using Mascot (Version 2.5.1, Matrix Sciences, London; 8039 entries). Carbamidomethylation of cysteines, tandem mass tags-sixplex amino terminus, tandem mass tags-sixplex lysine (for tandem mass tags labeled samples) were set as fixed modifications and methionine oxidation as variable modifications. Trypsin was selected as the enzyme, with two potential miss-cleavages. Peptide and fragment ion tolerances were respectively 10 ppm and 0.02 Da. Threshold of the average reporter signal-to-noise ratio (S/N) was set to 2. False discovery ratio (FDR) was set to 1% at peptide-spectrum match (PSM), peptide and protein levels. Only high-confident master proteins with at least two distinct peptide sequences were required for identifications. Obtained fold-changes were tested to be significantly different using two-way ANOVA tests in conjunction with Tukey’s multiple comparison test.

2.9. Data Processing and Analysis

AMOPLS analysis was conducted after unit variance scaling as previously described [14,23] under the MATLAB^® 8 environment (The MathWorks, Natick, MA, USA). The AMOPLS method is able to distinguish specific alterations related to different experimental factors by explicitly taking advantage of the structure of multifactorial designs. It associates ANOVA decomposition of the data matrix and Multiblock Orthogonal Partial Least Squares for an efficient assessment of both main effects and potential interactions related to the factors under study. Dedicated predictive components allow the impact of these effects to be easily displayed for interpretation purpose. Score plots are used to highlight potential groupings of observations according to the factor’s levels, while loading plots help to evaluate the contribution of each metabolite to the observed trends. Moreover, specific Variable Importance in the Projection (VIP²) can be calculated for each effect as a parameter facilitating the interpretation of the sources of variability in the model. A series of 10⁴ random permutations was performed to validate the AMOPLS model and assess the statistical significance of the effects. Metabolites selected according to their VIP² values in the AMOPLS model and proteins selected according to their fold-changes and p-values were passed to MetaCore™ v6.26 (Thomson Reuters) to perform joint pathway enrichment analysis.

3. Results and Discussion

In this work, a composite experimental strategy–a multifactor, multiplatform and multi-omics approach–was adopted to collect the broadest possible amount of information describing the normal physiology of the 3D rat brain cell culture system and the perturbations thereof upon TMT exposure.

3.1. Metabolomic Analyses

3.1.1. Sample Analysis, Data Acquisition and Preprocessing

In previous works we found, in good agreement with most of the existing literature, that protein precipitation with cold hydro-alcoholic mixtures provides a simple and convenient method to disrupt cells, stop enzymatic activity, and solubilize metabolites with a broad variety of polarity values [24,25,26,27]. The chosen biological model contained neurons, as well as glial cell populations (astrocytes, oligodendrocytes and microglia), which spontaneously aggregate into ~1 mm spheres displaying tissue-like organization. Such structures are robust enough as to demand a mechanical aid for complete disaggregation, so an additional sonication step was necessary to achieve complete disruption.

For untargeted metabolomic acquisition, complementary LC and ionization modes were employed. On one side, RPLC provides robust and well-stablished separation patterns for medium- to low-polarity compounds. Core-shell particles were selected based on our previous experience, since they allow for faster separations while keeping backpressure relatively low [28]. As many relevant metabolites present higher polarity values, RPLC was complemented with an amide-based HILIC separation. This approach had been previously applied to study neural metabolome and proved to afford an excellent separation of numerous polar and charged neurotransmitters and their precursors [23].

In both RPLC and HILIC conditions, detection was carried out in ESI+ and ESI– modes with independent analytical runs. For each one of the four LC-MS sequences, chromatograms were imported into Progenesis QI where retention time alignment, peak picking, adduct deconvolution and signal normalization were performed. LC-MS quality was assessed by means of technical quality control samples (QCs) prepared by pooling aliquots of extracts from all the samples in the study, and repeatedly injected during the analytical sequence. Features presenting a coefficient of variation above 35% in the QCs were filtered and not used downstream the metabolomics workflow. This filtering step was found to typically reduce the number of variables in the results matrix to about 10% of the original data. Metabolite annotation was performed using an in-house database containing the properties of more than 600 reference compounds analyzed in the same conditions. By matching retention time, accurate mass and isotopic pattern, highly reliable level 1 annotations of the metabolites were obtained [29]. In the cases where one metabolite was annotated by more than one analytical platform, the final choice was done on the base of the retention factor, peak intensity, signal-to-noise ratio and peak shape. The assignation of each metabolite to its best analytical platform avoids the occurrence of redundant variables that could affect multivariate data analysis. Eventually, a table of 36 samples described by 189 metabolites univocally annotated at level 1 across the four analytical platforms (33, 36, 58 and 62 for RPLC+, RPLC–, HILIC+ and HILIC–, respectively) was obtained.

3.1.2. Multivariate Analysis Using the AMOPLS Algorithm

The first step of the multivariate data analysis (MVA) was to decompose and quantify the specific variability related to the three factors involved in the DOE (maturation state, exposure duration and TMT concentration, as described in Figure 1). Taking into account the typical structure of datasets issued from metabolomics experiments, several approaches were proposed to deconvolute multiple experimental factors with more variables than observations, such as ANOVA-PCA, ANOVA-SCA or ANOVA-PLS [13]. The main limitation of these methods is that the separate analysis of each effect is carried out after ANOVA decomposition. As a consequence, several submodels are generated, thus making the global picture of all the experimental effects challenging to draw. To overcome this limitation and provide both an integrative model and effect-specific components, AMOPLS was recently proposed [14]. In this approach, the ANOVA-based partitioning of the sources of variation (experimental factors) is followed by a joint analysis of the ANOVA submatrices by means of a supervised multiblock algorithm based on the OPLS framework. Several recent publications reported the successful application of the AMOPLS strategy to the analysis of multifactorial omics experiments [23,30,31].

In

Table 1. Contributions of effects and their interactions, Residual Structure Ratio (RSR) and statistical significance (10⁴ permutations) issued from the AMOPLS model built for the metabolomics results.

	Maturation State	Exposure Duration	TMT Conc.	Maturation × Exposure	Maturation × Conc.	Exposure × Conc.	Residuals
Contribution	22.3%	15.3%	6.3%	13.5%	2.1%	3.8%	36.6%
RSR	2.42	1.77	1.37	1.49	1.01	1.07	1.00
p-value	0.01	0.01	0.01	0.01	83.53	29.10	-

, the role of the various effects under study (maturation state, exposure duration, TMT concentration) and their interactions are presented. A first finding is that the contribution from the residual accounts only for about one third (i.e., 36.6%) of the total variability, meaning that the biological model responded to a larger extent to the experimental factors under study (i.e., 63% of the observed variability). This underlines the relevance of well-defined experimental designs and the use of models able to incorporate them into data analysis when dealing with complex questions requiring a multifactorial approach.

The most important contributions came from the maturation state and the exposure duration, accounting for 22.3% and 15.3% of overall metabolic variability, respectively. In addition, the interaction between these two factors corresponds to 13.5%. This can be easily understood if the DOE is considered (Figure 1). Indeed, there is a parallel evolution of the maturation state and the exposure duration, both of them spanning across the culture duration to exert synergetic effects in the “days in vitro” time scale. Interestingly, when a PCA is performed on the data set (Figure 2), it can be seen that samples can be separated on the first principal component according to the age of the cells. This observation, together with the quantitative contribution of each factor and their interaction, made clear that major changes measured in the metabolic phenotype resulted from cell differentiation and maturation.

Even though the contribution of the TMT concentration (6.3%) represents the smallest observed effect among the three factors, it is highly statistically significant with a p-value lower than 0.01% (Table 1). One of the main goals of the study being to determine whether the response of the neural system to the TMT exposure was observed to be dependent on the maturation state or on the exposure duration, we prompted to check the interactions between these factors. Changes associated with the two other interactions, namely TMT Concentration × Maturation state and TMT Concentration × Exposure duration, were not found to be statistically significant. Hence, neither the maturation state nor the duration of the exposure were causing the cells to exhibit a specific metabolic phenotype upon exposure to the essayed TMT concentrations. Thus, we subsequently focused our efforts on the effect of the concentration of TMT.

Using the AMOPLS approach, samples can be neatly separated according to TMT concentration, avoiding the interference of the two other experimental factors (Figure 3). The fourth predictive component of the model (tp4) was mainly associated with TMT concentration (91.9% of the observed variability, Table S1). Since the Concentration main effect was considered at three different levels (0.0, 0.5 and 1.0 µM), two predictive components can actually be used to separate the groups. The second predictive component related to this effect was tp9, although having a notably smaller salience (23.1%). By inspecting the scores plot on Figure 3, it can be observed that samples coming from cells exposed to the same TMT concentration are well grouped, while clusters corresponding to each of the three TMT concentrations are well-separated, as opposed to what was observed when classical PCA analysis was used (Figure 2).

In a previous paper, we showed that in the framework of AMOPLS, VIP² values offer a flexible metrics allowing the quantitative comparison of several effects for a given variable [23]. Figure 4 shows the VIP² values calculated from the loadings of 60 metabolites having the largest contributions to the Concentration main effect in the AMOPLS model. This plot shows how, although the metabolites on the left side have the largest VIP² values for the Concentration effect (for instance N-acetyl aspartic acid, fumaric acid or succinic acid), there are actually other metabolites exerting an even larger global contribution to the ensemble of the model through the remaining effects and interactions (such as O-acetylserine or guanosine-5′-monophosphate). In the case of multifactorial experimental designs, this strategy allows therefore to deconvolute the relevance of each metabolite and experimental factor combination, making possible to focus on specific effects. Using the VIP² and loadings, it is possible to generate a list of identified metabolites that can be used for pathway analysis according to their specific variation towards each effect, as explained below.

To improve the understanding of the untargeted information about differentially expressed metabolites from a toxicological point of view, it is useful to put it in the context of metabolic pathways. By doing so, the biochemical and functional relationships existing between the molecules can be taken into account. This task is usually achieved by using metabolic enrichment or pathway analysis tools. Keeping in mind the possibility to merge the data obtained from different omics approaches for joint enrichment analysis, the MetaCore™ suite was chosen.

The probability of a set of up- or down-regulated metabolites to be erroneously assigned to a certain pathway or biological process is evaluated considering the size of each ensemble. The computed p-values (Figure 5) can then be inspected to rank the processes having the highest probability to be significantly altered under the studied conditions. The exposure of the cells to TMT resulted in dysregulated levels of a number of metabolites that were related to neuron differentiation and neurotransmission processes. Increasing doses of TMT caused a drop in the cellular concentrations of some of the metabolites involved in the tricarboxylic acid cycle, such as fumaric or succinic acids. This is one of the expected effects of exposure to TMT, since this substance is known to be able to exert an inhibitory effect on oxidative phosphorylation [32]. N-acetyl aspartic acid was also found to be dramatically decreased upon exposure to TMT. This metabolite plays a key role in neuronal cells homeostasis participating in osmolarity regulation and myelin synthesis, and supporting the onset of an impaired energetic metabolism [33]. Crucial precursors in neurotransmitter biosynthesis (phenylalanine, tyrosine) or neurotransmitters themselves (such as GABA) were also found in lower concentrations in cells exposed to TMT [34,35]. Modifications in biosynthesis and metabolism of glutamate and of GABA revealed by dysregulated GABA, 4-guanidinobutanoic acid, glutamic acid or 2-oxoglutaric acid, as well as other dysregulated processes involving neurotransmission functionality suggest a predominant deleterious effect of TMT exposure on neuronal cell populations, which is in line with the known neurotoxic effects of the compound [36]. Altered levels of dopamine, glutamic acid, and GABA may result from an adaptation of both neuronal differentiation processes and cholinergic and GABAergic signaling. General biosynthetic pathways involving coupling amino acids to tRNA were affected, as highlighted by the modified concentrations of tyrosine, leucine, glutamic acid, methionine, tryptophan and other amino acids.

3.2. Complementary Proteomic Analyses

In view of the results issued from the metabolomics analyses, complementary proteomics experiments were conducted to gain a deeper insight into the changes taking place due to the Concentration effect. To do so, a tandem mass tag-based proteomic experiment was conducted. The immature state was selected from the Maturation effect, and the level of the Exposure duration was fixed at 10-days with the highest concentration evaluated (1.0 µM). A total of 2422 proteins were identified by at least two unique peptides, with a 1% FDR-threshold at the levels of protein, peptide and peptide-to-spectrum matches (PSM). Systematic comparison of the various tested conditions with the corresponding untreated controls pointed up that the TMT treatment induced an increase of the glial fibrillary acidic protein (GFAP) and the vimentin, which are slightly upregulated (fold-changes of 1.24 and 1.14, respectively). These proteins, considered markers of astrocyte reactivity, underline the onset of the neuroinflammation triggered by TMT treatment. Other individual proteins such as the heme oxigenase 1 (HO-1) and the glutathione peroxidase (GPx1) are upregulated in the TMT-treated cells (fold-changes of 1.28 and 1.23, respectively). The regulation of these two proteins, as well as two metallothioneins, upon TMT treatment suggested a cellular response to oxidative stress. These results from individual protein expression changes in immatures cultures confirm that the exposure is triggering protein changes reflecting the induction of TMT-related biological processes, such as the onset of inflammation, and oxidative stress [37]. Data for the whole ensemble of identified and quantified proteins was subsequently analyzed on the basis of their fold-change by comparison of the 1.0 vs. the 0.0 µM TMT concentrations. Proteins showing a fold-change larger than 1.2 and a p-value below 0.05 were kept for pathway analysis.

3.3. Joint Multi-Omics Analyses

When applied on the same biological samples, each omics aims at studying the underlying phenomena from a different perspective. Numerous strategies were proposed and applied to merge multi-table and multi-omics datasets, each one having their own advantages and limitations [38,39,40,41]. The chosen platform relies on a manually curated database of biological maps, and it is able to simultaneously match various omics datasets against those maps.

Datasets containing significant data from metabolomics (AMOPLS VIP² > 1.1 for the Concentration effect) or proteomics (absolute fold-change > 1.2 and p-value < 0.05 for the two compared situations) were uploaded to the analysis tool. p-values calculation for the highlighted pathways or metabolic maps, along with false discovery rates, are efficient ways to sort and underline biological pathways with the highest biological interest. As a result, the top 10 enriched pathways are shown in Figure 6 and

Table 2. Detailed p-values, False Discovery Rate (FDR) and ratio of significantly altered elements (proteins or metabolites) with regard to the total number of elements contained in each map, for the ten most important metabolic maps found on the joint analysis of both systems biology approaches.

	Proteomics			Metabolomics
Maps	p-value	FDR	Found Elements	p-value	FDR	Found Elements
Neurophysiological process: GABA-A receptor life cycle	3.721 × 10−9	9.208 × 10−7	7/27	1.316 × 10−1	2.466 × 10−1	1/27
Neurophysiological process: Role of CDK5 in presynaptic signaling	4.924 × 10−9	9.208 × 10−7	7/28	9.239 × 10−3	4.070 × 10−2	2/28
Neurophysiological process: HTR2A signaling in the nervous system	1.177 × 10−8	1.761 × 10−6	8/48	2.221 × 10−1	2.784 × 10−1	1/48
Stem cells: Schema: Adult neuron differentiation in the Subventricular and Subgranular Zones	4.166 × 10−2	1.484 × 10−1	2/35	2.227 × 10−8	1.815 × 10−6	6/35
Neurophysiological process: Activity-dependent synaptic AMPA receptor removal	1.221 × 10−7	1.142 × 10−5	8/64	2.847 × 10−1	2.994 × 10−1	1/64
Nicotine signaling in dopaminergic neurons, Pt. 1 - cell body	9.946 × 10−3	6.526 × 10−2	3/48	1.598 × 10−7	7.598 × 10−6	6/48
Gamma-aminobutyrate (GABA) biosynthesis and metabolism	8.397 × 10−2	2.073 × 10−1	2/52	2.608 × 10−7	8.503 × 10−6	6/52
Histidine-glutamate-glutamine metabolism	2.240 × 10−1	3.870 × 10−1	2/96	6.016 × 10−7	1.613 × 10−5	7/96
Neurophysiological process: Constitutive and regulated NMDA receptor trafficking	2.185 × 10−6	1.634 × 10−4	7/65	4.498 × 10−2	1.094 × 10−1	2/65
Neurophysiological process: GABAergic neurotransmission	1.051 × 10−4	3.574 × 10−3	5/51	6.315 × 10−6	1.287 × 10−4	5/51

. Neurophysiological pathways are extensively represented in the resulting enrichment, which allows to highlight three out of the top 10 pathways linked to GABA metabolism or neurotransmission. It is interesting to note that the consequences of GABA perturbations may be different in immature cultures than in the differentiated ones, since from an excitatory and trophic role during early brain development, GABA switches to an inhibitory function in adult brain [42]. AMPA and NMDA-related biological pathways [43] are also emphasized in this enrichment analysis. Such processes are linked to glutamate and show that glutamatergic neurotransmission is altered by TMT treatment, as previously described [44].

In some maps, the ratio of significantly modified proteins was higher than the ratio of the significantly modified metabolites and the reverse was found in other maps, suggesting that the combination of the two omics approaches increases the level of confidence. Although neither the maturation state nor the duration of the exposure was causing the cells to exhibit a specific metabolic phenotype upon exposure to TMT, a previous analysis of proteomics data revealed that a repeated exposure increased the number of significantly modified proteins compared to the acute exposure. However, pathway analysis did not suggest different mechanism of toxicity between immature and differentiated cultures. However, TMT-induced alterations of the maturation process were more prominent in immature cultures and the reverse was observed for TMT-induced glial reactivity [37]. This demonstrates the complementarity of the use of the present workflow and of an unbiased combination of the omics techniques, as summarized in Figure S1, for a more subtle and comprehensive analysis of the mechanism of toxicity.

4. Concluding Remarks

In this work, a complete workflow allowing to explore data issued from multifactorial, multiplatform and multi-omics studies was proposed. By being able to combine complementary biochemical information retrieved by metabolomics, a better insight can be obtained in the study of complex cellular changes taking place upon exposure to molecules of toxicological relevance. The choice of exposing to a known neurotoxicant a 3D neural model undergoing progressive differentiation during time in culture allows to consider the relative effect of maturation stage, of treatment duration and of the concentration. The use of DOE to evaluate the relative contribution of various factors of crucial importance to toxicology studies (cell maturation, exposure duration and concentration) provides a way to comprehensively examine the different responses shown to the same toxic substance. For this purpose, AMOPLS allowed us to separately evaluate the contribution of each effect. Orthogonal LC-MS approaches complement each other to broaden the ability to inspect the sample’s chemical space, thus yielding a richer set of information about metabolic changes. The results revealed that neither the maturation state nor the exposure duration caused a significantly different response in the neural tissues to the different tested concentrations of TMT. The TMT concentration however played a crucial role in the response of the system, impairing multiple cellular processes related with neuronal differentiation, biosynthesis homeostasis and neurotransmission. Finally, metabolomic results were supported by proteomics analyses that could be integrated at the biological level by using a pathway enrichment analysis platform. The presented approach constitutes a good example of how the different tools that systems biology currently provides can be successfully combined in an integrative workflow to address the study of multifactorial toxicology experiments.