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.

   Keywords: metabolomics, proteomics, pathway analysis, multifactorial
   experiments, AMOPLS, multiplatform omics, toxicology, trimethyltin

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 [[42]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 [[43]2,[44]3,[45]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
   [[46]5,[47]6]. Proteins need a longer period to be either synthetized
   or degraded [[48]7,[49]8], while metabolomics can give a picture of
   both long-term and instantaneous reactions of the system induced by the
   stressor under study [[50]9,[51]10,[52]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 [[53]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
   [[54]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 [[55]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 ([56]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 [[57]15,[58]16,[59]17,[60]18]. Taking previous in-house
   experience with this substance in the 3D rat model [[61]19,[62]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 [[63]21]. Neural cells were seeded at a density of 6 × 10^7
   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[2] 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[2]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[2]O:MeCN 95:5 v/v (RPLC) or MeCN:H[2]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[2]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^−1, 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[2]O:MeCN (5:95, v/v) and solvent B was H[2]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^−1, 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^−1 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 [[64]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[2]O (solvent A) and MeCN (solvent B), both containing
   0.1% formic acid, at a flow rate of 300 nL min^−1. For MS survey scans,
   the resolution was set to 140,000 and the ion population was set to 3 ×
   10^6 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^5 (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 [[65]14,[66]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^2) can be
   calculated for each effect as a parameter facilitating the
   interpretation of the sources of variability in the model. A series of
   10^4 random permutations was performed to validate the AMOPLS model and
   assess the statistical significance of the effects. Metabolites
   selected according to their VIP^2 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 [[67]24,[68]25,[69]26,[70]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 [[71]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
   [[72]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 [[73]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 [74]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 [[75]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
   [[76]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 [[77]23,[78]30,[79]31].

Figure 1.

   [80]Figure 1
   [81]Open in a new tab

   Experimental sequence of the metabolomics study presented in this work.
   Numbers correspond to the Days In Vitro (DIV).

   In [82]Table 1, 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.

Table 1.

   Contributions of effects and their interactions, Residual Structure
   Ratio (RSR) and statistical significance (10^4 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 -
   [83]Open in a new tab

   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 ([84]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 ([85]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.

Figure 2.

   [86]Figure 2
   [87]Open in a new tab

   Scores plot obtained when the metabolomics experimental data is
   submitted to conventional Principal Component Analysis (PCA). Data has
   been centered and UV-scaled.

   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% ([88]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 ([89]Figure 3). The fourth predictive component of
   the model (tp4) was mainly associated with TMT concentration (91.9% of
   the observed variability, [90]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 [91]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 ([92]Figure 2).

Figure 3.

   [93]Figure 3
   [94]Open in a new tab

   Scores plot obtained from the AMOPLS model when the predictive
   components bearing the variability coming from the Concentration effect
   are selected (tp4 and tp9, according to [95]Table S1). Note how the
   samples are neatly separated into the three levels for the desired
   effect, as opposed to [96]Figure 2.

   In a previous paper, we showed that in the framework of AMOPLS, VIP^2
   values offer a flexible metrics allowing the quantitative comparison of
   several effects for a given variable [[97]23]. [98]Figure 4 shows the
   VIP^2 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^2 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^2 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.

Figure 4.

   [99]Figure 4
   [100]Open in a new tab

   VIP^2 values calculated for all of the statistically significant
   effects and interactions in the model. Only the 60 metabolites showing
   the largest VIP^2 values for the TMT Concentration effect are
   presented.

   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
   ([101]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 [[102]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 [[103]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 [[104]34,[105]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 [[106]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.

Figure 5.

   [107]Figure 5
   [108]Open in a new tab

   Main biological processes altered as a consequence of TMT exposure,
   independently of the maturation state or the exposure duration.
   p-values were calculated from the metabolomics data using the MetaCore™
   pathway analysis tool.

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 [[109]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 [[110]38,[111]39,[112]40,[113]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^2 >
   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
   [114]Figure 6 and [115]Table 2. 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 [[116]42]. AMPA and NMDA-related biological pathways
   [[117]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
   [[118]44].

Figure 6.

   [119]Figure 6
   [120]Open in a new tab

   List of the top-10 metabolic maps as ranked according to the joint
   contribution of proteomics and metabolomics results. Bars represent the
   p-value for a particular map and for each systems biology approach.
   These values correspond to the probability for mistakenly accepting
   that a map is altered upon evaluation of the altered metabolites or
   proteins assigned to each map.

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
   [121]Open in a new tab

   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 [[122]37]. This
   demonstrates the complementarity of the use of the present workflow and
   of an unbiased combination of the omics techniques, as summarized in
   [123]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.

Acknowledgments