Antagonistic Interactions in Mitochondria ROS Signaling Responses to Manganese

Antagonistic interaction refers to opposing beneficial and adverse signaling by a single agent. Understanding opposing signaling is important because pathologic outcomes can result from adverse causative agents or the failure of beneficial mechanisms. To test for opposing responses at a systems level, we used a transcriptome–metabolome-wide association study (TMWAS) with the rationale that metabolite changes provide a phenotypic readout of gene expression, and gene expression provides a phenotypic readout of signaling metabolites. We incorporated measures of mitochondrial oxidative stress (mtOx) and oxygen consumption rate (mtOCR) with TMWAS of cells with varied manganese (Mn) concentration and found that adverse neuroinflammatory signaling and fatty acid metabolism were connected to mtOx, while beneficial ion transport and neurotransmitter metabolism were connected to mtOCR. Each community contained opposing transcriptome–metabolome interactions, which were linked to biologic functions. The results show that antagonistic interaction is a generalized cell systems response to mitochondrial ROS signaling.


Introduction
Technological advances to process large biological data sets are projecting a dynamic and critical change into the global analysis of environmental stressors and disease [1]. These advances transform mechanistic research by enabling the study of effects on multiple targets as determinants of outcome rather than limiting hypotheses to singular determinants. In human, murine and cell studies, this multi-targeted approach reveals candidate regulatory determinants with intervention potential. Integrated omics studies show that a synergistic stress response hub was identified for combined fungicide (maneb) and herbicide (paraquat) exposure linked to Parkinson's disease. In this published study, PPAR-γ and NRF2 genes were associated with protective responses, while pro-apoptotic genes were associated with damage response. These adverse and beneficial outcomes were interpreted based on known gene functionality [2] and not by measured functional parameters within the cell. Other network-based studies showed immune networks of genes and metabolites that correlate with later adaptive response to herpes zoster vaccination [3]; that network structures connect inflammation to lung metabolomics following viral infection [4,5]; and that redox proteomics in cadmium toxicity links fatty acid metabolites to respective mitochondrial fatty acid enzymes and transport systems [6,7]. In principle, such approaches can disentangle complex relationships among beneficial and harmful signaling responses to environmental exposures, an emerging need for exposome research [8][9][10].
Mitochondrial connections with other cell functions are especially important because mitochondria are susceptible to environmental agents and are essential for genetic and metabolic signaling to maintain and integrate energy, fatty acid and amino acid metabolism [11,12]. Although often conceptualized as occurring through discrete pathways and studied by isolating and targeting specific mitochondrial sites, stressors involving reactive oxygen species (ROS) impact multiple targets within and outside of mitochondria. This results in multiple signals, some with opposing beneficial and adverse effects. Such an antagonistic interaction is well-recognized in genetics [13,14] and aging [15][16][17][18], but has been little studied in exposome research.
Manganese (Mn) is an essential trace element with characteristics suitable to study antagonistic interactions from environmental exposures. Mn is required at low doses as a co-factor for the antioxidant mitochondrial protein superoxide dismutase-2 (SOD2), as well as by many other enzymes. Mn at high doses causes neurotoxicity (Manganism, neurodevelopmental defects, intellectual cognitive deficits, parkinsonian symptoms) due to occupational and environmental exposures [19][20][21][22][23]. Within the cell, Mn undergoes redox reactions and stimulates mitochondrial oxidant production [24][25][26][27]. Mn is similar in size and charge to other metals, such as iron and calcium, leading to alteration in their functions [28][29][30] and inducing further oxidative stress. Mn also interacts with electron transport chain complexes and impedes mitochondrial oxidative phosphorylation [24,27,31,32], increasing mitochondrial oxidative stress. Previous studies from our group have shown linear responses of the external trigger Mn directly on genes [33] and metabolites [34] separately, without considering both linear and nonlinear functional parameters as intermediate signaling and their role in gene-metabolite interactions.
To test for antagonistic interactions, we used a data-driven approach with xMWAS [12,35] to integrate metabolomics and transcriptomics with functional measures of mitochondrial activities in human neuroblastoma SH-SY5Y cells treated with increasing Mn concentrations for 5 h. Raw data from three previously published distinct studies [27,33,34] were used to identify mitochondrial-cellular signaling (MCS) subnetwork structures in response to Mn ( Figure S1). The resulting ranked correlations of metabolites and transcripts with functional parameters such as cellular Mn, cellular thiols, mitochondrial ROS (mtROS), mitochondrial H 2 O 2 (mtH 2 O 2 ), SOD2 activity, mitochondrial basal oxygen consumption rate (OCR) and proton leak OCR (mtOCR) were used with community detection tools and targeted analyses to reveal subnetwork structures that mapped to known adaptive and adverse toxicologic responses. The results establish that the causative agent Mn affects functional responses, resulting in adverse and beneficial mitochondria-cell signaling mechanisms that occur simultaneously and provide a generalizable approach to study antagonistic interactions at a systems level.

Cell Culture, Mn Dose Treatment, Big Data Acquisition
Raw data (mitochondrial cellular parameters, metabolomics, RNA seq) from three previously published distinct studies [27,33,34] were used with xMWAS to identify MCS subnetwork structures in response to a Mn dose in human neuroblastoma cells (SH-SY5Y). The study design is provided in a schematic ( Figure S1). Cellular Mn concentrations were measured by graphite furnace atomic absorption spectroscopy (GFAA) (Thermo Scientific, Waltham, MA, USA, iCE3000 Series), total cellular thiols were measured using Ellman's reagent, mitochondrial function assay was measured using WST-1(Roche), global mitochondrial oxidant production was measured using Mitotracker Red CM-H2XRos (Thermo scientific), mitochondrial H 2 O 2 production was measured using mitochondrial peroxy yellow 1 (MitoPY1) (Tocris Biosciences, Minneapolis, MN, USA) and mitochondrial respiration was analyzed using a Seahorse Bioscience XF96 extracellular flux analyzer (North Billerica, MA, USA). Median value for each functional parameter at each Mn dose was obtained, and the resultant functional output matrix was used in downstream omics integration analysis ( Figure S1).
Metabolomics analyses with three technical replicates were performed with ultrahigh-resolution mass spectrometry with hydrophilic interaction liquid chromatography  [36] and xMSanalyzer [37] to yield datasets with mass to charge (m/z), retention time (RT, s) and intensity for each mass spectral feature for all samples. The data were further filtered so 80% of all the samples had an intensity, quantile normalized and log2 transformed. Median value for each Mn dose was obtained, and the resultant metabolome data matrix was used in downstream omics integration analysis ( Figure S1).
RNA extraction, RNA-Seq analysis and further data processing to obtain transcriptomics read counts has been previously described [33]. Raw data are made accessible to the public at Geo-NCBI with accession GSE129336. Briefly, total RNA was extracted using miRNeasy Mini Kit (Qiagen) and processed for RNA-Seq analysis following quality control checks for each sample. The platform used for sequencing was the Illumina HiSeq 2500. The data were normalized by variance stabilization transformation as implemented in the DESeq2 R/Bioconductor package [38]. Median value for each Mn dose was obtained, and the resultant transcriptome data matrix was used in downstream omics integration analysis ( Figure S1).

Integrated Multi-Omics Network Analysis
The MCR dataset, HRM dataset and transcriptome dataset for Mn dose response were integrated using an xMWAS v0.54 tool based on the partial least squares (PLS) regression method for pairwise data integration [12]. xMWAS is a web-based data integration software used for differential network analysis of two or more omics platforms with biochemical responses available at https://github.com/kuppal2/xMWAS accessed on June 2020. The 3 datasets used as input for xMWAS included data for 5 Mn doses along with control and consisted of-1. Mitochondrial-cellular parameters with 7 measures (x) including cellular responses such as intracellular Mn concentration (Mn), total cellular thiols (thiols) and mitochondrial-specific responses such as ROS production (mtROS), H 2 O 2 production (mtH 2 O 2 ), SOD2 activity, basal oxygen consumption rate (bOCR) and proton leak dependent oxygen consumption rate (pOCR) measures; 2. Metabolome data matrix (y) with 6296 metabolic features; 3. Transcriptome data matrix (z) with 19,965 transcripts ( Figure S1). Correlation threshold criteria for determining significant associations were set at (|r| > 0.8) and p < 0.05. Community detection was also performed using the multilevel community detection method [39], wherein a community consists of tightly connected biomolecules, herein genes, metabolites and the MCR parameters. The resultant multi-data integrative network provided by xMWAS allowed for identification and visualization of associations between mitochondrial functional outputs, genes and metabolites.

Pathway Analysis for Selected Metabolic Features and Genes
Metabolic features, selected based on the correlation threshold criteria, were subjected to permutation testing (p < 0.05) in pathway enrichment analysis using mummichog v1.0.9 [40]. Only metabolic pathways including greater than 3 overlap sizes were selected for further analysis. Genes selected based on the correlation threshold criteria were examined for enrichment of known biological processes in human species using DAVID v6.8, the Database for Annotation, Visualization and Integrated Discovery. The functional classification tool used generated a gene-to-gene similarity matrix based on shared functional annotation by using over 75,000 terms from 14 functional annotation sources [41].

Metabolite Annotation and Identification
Metabolites selected based on pathway analysis were annotated using an in-house library of confirmed metabolites and xMSannotator v1.3.2 [42] with Human Metabolome Database v3.5. Level 1 metabolite identification by the criteria of Schymanski et al. (Schymanski et al., 2014) was performed by co-elution relative to authentic standards and ion dissocia-tion mass spectrometry. Additional annotations (levels 3, 4 and 5 metabolite identification criteria by [43]) were made using medium to high confidence scores from xMSannotator and using the KEGG (Kyoto Encyclopedia of Genes and Genomes) [44] and HMDB (Human Metabolome Database) [45] databases at m/z match of 10 ppm tolerance. Confidence scores for annotation by xMSannotator were derived from a multistage clustering algorithm based on correlation, co-elution, network modularity, and adducts/isotopes patterns. Additional information on metabolites discussed in the figures is plotted in Table S9.

Subcellular Fractionation and Western Blotting
To examine nuclear translocation of p65 NF-κB by Mn and Mn + MitoQ, cells were exposed to Mn doses (0, 10 and 100 µM for 5 h) with or without prior treatment with 1 µM MitoQ or vehicle control for 30 min. Post-exposure, cells were isolated for subcellular fractionation using a nuclear and cytoplasmic extraction kit (Thermo Scientific, Rockford, IL, USA). Isolated fractions were confirmed by Western blotting using antibodies against p65 NF-κB. Controls were determined by probing nuclei for lamin A/C. All antibodies were obtained from Cell Signaling Technology, Boston, MA, USA. An IRDye 800 conjugated affinity purified anti-rabbit (Rockland Immunochemicals, Gilbersville, PA, USA) or Alexa-Fluor-680conjugated anti-mouse (Invitrogen, Waltham, MA, USA) was used as secondary antibody. Bands were visualized using an Odyssey scanner (Li-Cor Biosciences, Lincoln, NE, USA) and Odyssey 2.1 software (Li-Cor, Lincoln, NE, USA) and quantified using ImageJ v1.48.
All statistical analyses, boxplots and heatmaps were generated in R v3.2.3. Schematic figures were created using online software at BioRender.com. Abbreviations used throughout the text are provided in Supplementary Text S1.

Data-driven integration of mitochondrial activity with transcriptome and metabolome.
Partial least squares (PLS) regression in xMWAS provided ranked correlations for the interactions of seven mitochondrial and cellular parameters (MCP) ( Figure S1). From these correlations, community detection allowed visualization of three major communities linked to cellular Mn (Community one, Figure 1A), mitochondrial oxidative stress [mtOx (mtROS, mtH 2 O 2 , SOD2 activity)] and total cellular thiols (Community two; Figure 1B), and mitochondrial oxygen consumption rate [mtOCR (basal OCR and proton leak OCR)] (Community three; Figure 1C), which were preserved at increased stringency (|r| > 0.9; Figure 1, center). The composition of the communities at |r| > 0.8 included 480 metabolites and 947 genes (Figure S1), with genes and metabolites differing for Community one (Tables S1 and S2), Community two (Tables S3-S6) and Community three (Tables S7 and S8). Community two had two sub-clusters centered on mtOx (Community 2a; mtROS, mtH 2 O 2 and SOD2 activity) and cellular thiols (Community 2b), which were largely negatively associated with each other.
From a top-down view, Community one included transcripts and metabolites most closely associated with cellular Mn, perhaps reflecting non-mitochondrial mechanisms of toxicity and homeostasis. The structures of Community two contained elements linked to mitochondrial oxidative stress, the main mechanism of cell toxicity. To validate this interpretation, we performed a targeted experiment with mitoQ, an inhibitor of mitochondrial oxidative stress, and found that NF-κB translocation to the nucleus, a well-established response to oxidative signaling, was blocked (see below) ( Figure S2). Community two had opposing substructures for mitochondrial ROS (Community 2a) and cellular thiol content (Community 2b), consistent with known causal mechanisms in which thiols are oxidized in response to oxidative stress. Community three included elements more strongly associated with mitochondrial respiration than with Mn or oxidative stress. For every gene and metabolite, the direction of positive and negative correlation with the functional measure is provided as [+] and [−], respectively, throughout the text. Directionality is also provided for the rest of genes and metabolites represented by correlation coefficient, presented in the supplementary tables as +|r| or −|r|. in response to oxidative stress. Community three included elements more strongly associated with mitochondrial respiration than with Mn or oxidative stress. For every gene and metabolite, the direction of positive and negative correlation with the functional measure is provided as [+] and [−], respectively, throughout the text. Directionality is also provided for the rest of genes and metabolites represented by correlation coefficient, presented in the supplementary tables as +|r| or −|r|.

Figure 1. Data-driven integrative analysis defined 3 communities with transcriptome and metabolome interaction with mitochondrial-cellular parameters [MCP] in response to Mn dose.
Three major communities generated through PLS regression method for pairwise data integration using xMWAS v0.54 at |r| > 0.9. (A) Community 1 (yellow)-environmental stressor subnetwork comprising intracellular Mn positively associated with 79 genes and 51 metabolites and negatively associated with 144 genes and 116 metabolites. (B) Community 2 (green)-oxidative stress subnetwork further bifurcated into two sub-communities at |r| > 0.8-communities 2a and 2b. Community 2a (dark green)-mitochondrial oxidative stress (mtOx) subnetwork community comprising SOD2, mtROS and mtH2O2 positively associated with 163 genes and 88 metabolites and negatively associated with 268 genes and 111 metabolites. Community 2b (light green)-cellular oxidative stress subnetwork comprising total cellular thiols positively associated with 270 genes and 56 metabolites and negatively associated with 115 genes and 107 metabolites. (C) Community 3 (blue)-oxidative phosphorylation subnetwork (mtOCR) comprising basal OCR and proton leak dependent OCR, positively associated with 120 genes and 109 metabolites and negatively associated with 172 genes and 56 metabolites. Seven MCPs are represented as rectangles, genes as triangles and metabolites as circles. Red edges represent positive, while blue represents negative associations with their respective MCPs. |r| > 0.8, p < 0.05.

Community one: Cell Mn was associated with biological processes for biosynthesis, morphogenesis, amino acid metabolism, differentiation and apoptosis.
We examined each community for gene-metabolite associations that could be traced to signaling mechanisms. Previous research showed adaptive responses to non-toxic Mn exposure in amino acid metabolism and protein secretion pathways, while toxic Mn exposures disrupted fatty acid and energy metabolism [33,34]. Integrated TMWAS showed

Figure 1. Data-driven integrative analysis defined 3 communities with transcriptome and metabolome interaction with mitochondrial-cellular parameters [MCP] in response to Mn dose.
Three major communities generated through PLS regression method for pairwise data integration using xMWAS v0.54 at |r| > 0.9. (A) Community 1 (yellow)-environmental stressor subnetwork comprising intracellular Mn positively associated with 79 genes and 51 metabolites and negatively associated with 144 genes and 116 metabolites. (B) Community 2 (green)-oxidative stress subnetwork further bifurcated into two sub-communities at |r| > 0.8-communities 2a and 2b. Community 2a (dark green)-mitochondrial oxidative stress (mtOx) subnetwork community comprising SOD2, mtROS and mtH 2 O 2 positively associated with 163 genes and 88 metabolites and negatively associated with 268 genes and 111 metabolites. Community 2b (light green)-cellular oxidative stress subnetwork comprising total cellular thiols positively associated with 270 genes and 56 metabolites and negatively associated with 115 genes and 107 metabolites. (C) Community 3 (blue)-oxidative phosphorylation subnetwork (mtOCR) comprising basal OCR and proton leak dependent OCR, positively associated with 120 genes and 109 metabolites and negatively associated with 172 genes and 56 metabolites. Seven MCPs are represented as rectangles, genes as triangles and metabolites as circles. Red edges represent positive, while blue represents negative associations with their respective MCPs. |r| > 0.8, p < 0.05.

Community one: Cell Mn was associated with biological processes for biosynthesis, morphogenesis, amino acid metabolism, differentiation and apoptosis.
We examined each community for gene-metabolite associations that could be traced to signaling mechanisms. Previous research showed adaptive responses to non-toxic Mn exposure in amino acid metabolism and protein secretion pathways, while toxic Mn exposures disrupted fatty acid and energy metabolism [33,34].  [48,49] varied negatively with imidazole acetate [−] and N-acetylglutamate [−] ( Figure 2C). The gene-metabolite pair is based on highest correlation values between the two and their enrichment in the pathway analysis ( Figure 2C). The results show an apparent antagonistic interaction associated with imidazole acetate and N-acetylglutamate and illustrate that Community one consists of opposing transcriptome-metabolite interactions ( Figure 2D). At the transcriptome level, Community one showed that cellular Mn was associated with changes in 223 genes ( Figure 1A), including 176 mRNA, 31 long non-coding At the transcriptome level, Community one showed that cellular Mn was associated with changes in 223 genes ( Figure 1A), including 176 mRNA, 31 long non-coding (lncRNA), 15 miRNA (microRNA) and 1 sncRNA (Y RNA) (Table S1). Biologic processes associated with mRNAs included cell proliferation, anatomical morphogenesis, growth regulation, cell migration and apoptosis (Figure 2A). Potentially adverse changes included the decrease in PAICS and CDH22 as well as FYN and FRK members of Src family of tyrosine kinases functioning in cell growth and proliferation [50,51] (Table S1).
Intracellular Mn was associated with changes in 167 metabolites ( Figure 1A and Table S2), which showed broad impact by pathway enrichment analysis on amino acid pathways (His; Lys; Gly, Ser, Ala and Thr; Trp; Ala and Asp; Met and Cys) and the nicotinamide pathway ( Figure 2B). Imidazole acetate is involved in His metabolism, and N-acetylglutamate is an activator of carbamoyl phosphate synthase I in the urea cycle ( Figure 2B) [52,53]; both were decreased by Mn. Hexose-bisphosphate, a precursor for glycolysis and the pentose phosphate pathway, and riboflavin, involved in energy pathways, were also decreased by Mn (Table S9).
Community two: Mitochondrial oxidative stress (mtOx) was associated with inflammation and fatty acid metabolism.
Mitochondrial oxidative stress represents a key mechanism in neurotoxicity caused by excess Mn. In the xMWAS, all three mitochondrial oxidative stress measures (mtROS, mtH 2 O 2 and SOD2) were tightly correlated in a sub-community, Community 2a. Cell thiol loss, a downstream consequence of increased ROS, partially separated as Community 2b with negative associations from Community 2a. Examination of the TMWAS for Community 2a showed that the proinflammatory transcript IKBKG (inhibitor of kappa light polypeptide gene enhancer in B-cells, kinase gamma) [+] encoding a protein that activates NF-κB [54] ( Figure 3A) was positively associated with the anti-inflammatory glycerophosphoinositol [+] ( Figure 3B,C). As a proof of concept experiment to show that data-driven associations without using metabolic perturbations can be validated in a system, additional experiments were conducted. IKBKG (NEMO) encodes a protein that activates NF-kB, which leads to NF-kB translocation to the nucleus and induces inflammatory signal [55]. Our current network analysis demonstrates Mn-dependent mt-Ox measures are positively associated with IKBKG gene, suggesting that Mn-dependent mitochondrial ROS-production could lead to NF-kB translocation. To demonstrate this, Mn exposed cells were treated with and without MitoQ, a mitochondrial targeted antioxidant, to decrease mitochondrial oxidative stress. Cells exposed to Mn indeed showed NF-kB translocation, which was further blocked by Mito-Q. These results demonstrate and validate our network analysis interpretation that Mn-dependent mitochondrial ROS-production indeed leads to NF-kB translocation and potential inflammatory signaling ( Figure S2). In addition to IKBKG, transcript data for other NF-kB components are presented ( Figure S3). Thus, network analysis can reveal hidden relationships just by functional measures integrated with omics data.
In the opposite direction, the essential fatty acid linolenate [−] associated negatively with IKBKG [+] ( Figure 3C, Table S4). Opposite effects were also seen when comparing glycerophosphoinositol [+] with the phospholipase gene LIPG [−] ( Figure 3C, Table S4). CPT1A, encoding a key protein to support mitochondrial fatty acid oxidation, was increased with mtOx. Thus, pro-inflammatory processes were accompanied by opposite effects on fatty acid and lipid metabolism.
In addition to the general pattern of increased inflammatory gene transcripts with mtOx, including IKBKG ( Figure 3A (Table S3), other transcripts in Community 2a (Table S3) were negatively associated with mitochondrial measures of oxidative stress ( Figure 1C). These included decreased WNT1, a transcript for a protein that inhibits inflammation and is neuroprotective [56] (Figure 3A). Community 2b was centered on cell thiols, and elements within this subcommunity were frequently associated negatively with elements in Community 2a. TMWAS for Community 2b also included opposing responses, such as for CLIP1 (CAP-Gly domain containing linker protein 1) [+], encoding a protein involved in microtubule dynamics and  Table S6). GPX8 (glutathione peroxidase 8) [−], encoding a protein for the elimination of H 2 O 2 from the ER, and DPYSL3 (dihydropyrimidinase-like 3) [+], encoding a protein functioning in cytoskeletal organization and axonal guidance, varied in opposite directions in association with α-aminoadipate ( Figure 4C). FER (tyrosine kinase) [+], a gene supporting synapse organization and generation of excitatory postsynaptic currents, was also negatively associated with thiols (Table S5).

PEER REVIEW 9 of 18
central toxicologic mechanism of Mn, therefore show patterns of antagonistic interaction within signals associated with increased ROS and within signals associated with loss of thiols.  At the metabolomic level, pathway enrichment showed that the metabolites associated with cell thiols ( Figure 1B, Table S6) were present in amino acid (e.g., alanine, glycine), niacin (nicotinate, nicotinamide) and purine and pyrimidine pathways ( Figure 4B, Table S6). The overall results for Community 2 ( Figure 4D), which represents the central toxicologic mechanism of Mn, therefore show patterns of antagonistic interaction within signals associated with increased ROS and within signals associated with loss of thiols.
Community three: Mitochondrial respiration was associated with ion transport and neurotransmitter metabolism.
Maintenance of mitochondrial function and cellular energetics is critical to protect against Mn toxicity. Basal OCR and proton leak OCR each had a biphasic response to Mn with stimulation at adaptive concentrations and decrease at toxic concentrations [27]. In xMWAS, genes and metabolites associated with basal and proton leak OCR clustered together in a tight mtOCR Community three ( Figure 1C). TMWAS for Community three showed that when mtOCR was highest, the abundance of CACNA1E [−], a subunit of a high voltage activated calcium channel that regulates neurotransmitter signaling ( Figure 5A), and the neurotransmitter, GABA [−] ( Figure 5B, Table S8), were lowest. CACNA1E modulates neuronal firing patterns and is associated with encephalopathies [57][58][59]. In addition to neurotransmitter signaling, an imbalance in bioenergetics also leads to consumption of alternative fuel sources such as GABA to produce TCA cycle intermediates through the GABA shunt [60]. The decrease in CACNA1E abundance while the neurotransmitter GABA was also decreased with increased mitochondrial oxidative phosphorylation, suggests an antagonistic response in the regulation of neuronal firing, and neurotransmitter-related metabolic and intramitochondrial function.
Negative associations were observed for the fatty acid, stearate and transcripts for iron sulfur assembly (IBA57) [−] and NADH reductase (CYB5R3, cytochrome b5 reductase 3) [−] ( Figure 5C), showing that opposing responses integrate neurotransmitter systems and other homeostatic mechanisms ( Figure 5D). Transcripts with protective functions that were positively associated with mtOCR included UPP1 (uridine phosphorylase 1) [+], an ATPdependent transcriptional regulator involved in pyrimidine salvage (Table S7). Transcripts that were negatively associated with the mtOCR included KCNC4 [−], a potassium channel, and SIRT6 [−], an NAD-dependent protein deacetylase involved in DNA repair (Table S7). Pathway enrichment analysis showed associations with multiple adaptive response pathways, such as cell cycle, intracellular signal transduction, transcription and ion transport ( Figure 5A). Pathway enrichment analysis for metabolomics demonstrated glutamate as the topmost pathway, followed by fatty acid metabolism, vitamin A, pyrimidine, urea cycle and butanoate metabolism ( Figure 5B and Table S9).
Overall, Community three associations linked to mitochondrial bioenergetic activity represent positive and protective response to increased Mn (Community one), as well as to oxidative stress (Community 2), but the results show that this occurs functionally through response structures with antagonistic signaling. Thus, the results from this community also support the conclusion that antagonistic interactions occur at sub-network levels downstream to the independent variable, Mn.
x FOR PEER REVIEW 10 of 18 addition to neurotransmitter signaling, an imbalance in bioenergetics also leads to consumption of alternative fuel sources such as GABA to produce TCA cycle intermediates through the GABA shunt [60]. The decrease in CACNA1E abundance while the neurotransmitter GABA was also decreased with increased mitochondrial oxidative phosphorylation, suggests an antagonistic response in the regulation of neuronal firing, and neurotransmitter-related metabolic and intramitochondrial function. Negative associations were observed for the fatty acid, stearate and transcripts for iron sulfur assembly (IBA57) [−] and NADH reductase (CYB5R3, cytochrome b5 reductase 3) [−] ( Figure 5C), showing that opposing responses integrate neurotransmitter systems and other homeostatic mechanisms ( Figure 5D). Transcripts with protective functions that were positively associated with mtOCR included UPP1 (uridine phosphorylase 1) [+], an ATP-dependent transcriptional regulator involved in pyrimidine salvage (Table S7). Transcripts that were negatively associated with the mtOCR included KCNC4 [−], a potassium channel, and SIRT6 [−], an NAD-dependent protein deacetylase involved in DNA repair (Table S7). Pathway enrichment analysis showed associations with multiple adaptive response pathways, such as cell cycle, intracellular signal transduction, transcription and ion transport ( Figure 5A). Pathway enrichment analysis for metabolomics demonstrated glutamate as the topmost pathway, followed by fatty acid metabolism, vitamin A, pyrimidine, urea cycle and butanoate metabolism ( Figure 5B and Table S9).
Overall, Community three associations linked to mitochondrial bioenergetic activity represent positive and protective response to increased Mn (Community one), as well as to oxidative stress (Community 2), but the results show that this occurs functionally through response structures with antagonistic signaling. Thus, the results from this community also support the conclusion that antagonistic interactions occur at sub-network levels downstream to the independent variable, Mn.   oxidative phosphorylation subnetwork associated with signaling markers for ion transport, neurotransmitter metabolism and iron sulfur cluster assembly. Pathway enrichment analysis of genes and metabolites associated with basal OCR and proton leak dependent OCR (mtOCR) subnetwork were studied. (A) Top pathways enriched by mtOCR-dependent genes at p < 0.05, represented as −log 10 P. Expression of individual genes associated with biphasic mtOCR from these enriched pathways presented as boxplots include genes involved in ion transport channels for calcium, CACNA1E; potassium, KCNC4; iron sulfur cluster assembly gene, IBA57 and NADHdependent reductase of redox cyclers, CYB5R3. (B) Top pathways representing mtOCR-dependent metabolites at p < 0.05, represented as −log 10 P. Raw intensity of individual metabolites from these enriched pathways presented as boxplots include neurotransmitter, GABA and long chain fatty acid, stearic acid. (C) Antagonistic interactions of genes and metabolites in mitochondrial oxidative phosphorylation cluster. (D) Schematic illustrating biphasic mtOCR-associated key drivers of ion transport, neurotransmitter metabolism and iron sulfur cluster assembly. The mtOCR intensity on the x-axis is indicative of Mn exposure low-med-high from left to right. Detailed information on genes and metabolites from Community 3 provided in Tables S7-S9. |r| > 0.8, p < 0.05.

Antagonistic patterns of non-protein coding transcripts associated with neurodegenerative diseases.
The capability to delineate opposing signaling processes raised the possibility that the integrative omics analysis using xMWAS described here could be used to evaluate antagonistic interaction in therapeutics. To test the feasibility of this concept, we considered Mn dosing as a model for a therapeutic treatment. We examined non-coding RNAs (ncRNA) associated with each community (highlighted in Tables S1, S3, S5 and S7) for ones involved in neurological dysfunction. ncRNAs, such as microRNA (miRNA) and lncRNA, are cellular regulators in various disease states including neurological dysfunction, cancer and oxidative stress [61,62]. Community-one-associated miRNA involved in regulation of autophagy-related protein PINK1 (PTEN-induced kinase 1, linked to Parkinson's disease), MIR27A, increased with intracellular Mn ( Figure 6A).
Community-2a-associated ncRNAs such as stroke-related MIR149 and Alzheimer'srelated lncRNA, KIAA0125, increased with the increase in mtOx, and neurobehavioralassociated lncRNA, PAX8-AS1, declined with the increase in mtOx ( Figure 6B). Community 2b was associated with Alzheimer's-related lncRNA, EPHA1-AS1, which declined with the decrease in total cellular thiols ( Figure 6C). Community-three-associated ncRNAs included Huntington's disease-related MIR9-1, which was negatively associated with mtOCR, and anti-inflammatory-associated MIRLET7A1 and cell migration lncRNA DEPDC1-AS1, which were positively associated with mtOCR ( Figure 6D). The results show that application of xMWAS to interrogate antagonistic interactions has the potential to allow deconvolution of very complex biologic processes linked to ncRNAs and provide a strategy to identify disease-specific therapeutic targets, as well as monitor therapeutic response to intervention. ones involved in neurological dysfunction. ncRNAs, such as microRNA (miRNA) and lncRNA, are cellular regulators in various disease states including neurological dysfunction, cancer and oxidative stress [61,62]. Community-one-associated miRNA involved in regulation of autophagy-related protein PINK1 (PTEN-induced kinase 1, linked to Parkinson's disease), MIR27A, increased with intracellular Mn ( Figure 6A).

Discussion
Complex molecular interactions drive biological processes in response to environmental exposures. Utilizing an advanced multi-omics approach with data-driven bioinformatics analysis, the present study demonstrates that metabolomics and transcriptomics changes can be combined with cellular-and mitochondrial-specific responses to provide a systems view of underlying mechanisms relevant to neurotoxicity. We show that at multiple levels of network structure, the genes and metabolites involved with specific functional traits exhibit characteristics which we interpret as antagonistic interaction, a characteristic previously found for mitophagy [63], viral infections [64], antioxidant and redox systems and aging [18,65]. The present bioinformatics approach shows that opposing gene and metabolite responses occur directly related to the causal agent, Mn. Some of these antagonistic interactions reflect mechanisms of mitochondrial oxidative stress, the consequent oxidation of cell thiols and the downstream protective responses of mitochondrial bioenergetics.
Antagonistic interactions occur because Mn is a trigger that generates both positive and negative responses. Within the Mn network structure, genes such as NFIB, which is critical for cell survival, differentiation and migration in cortical and forebrain development [66], also include genes such as BCL6 functioning in apoptosis. Complex interactions of antagonistic responses in a high voltage activated calcium channel controlling neuronal firing patterns that are associated with migraines and epileptic encephalopathies [57][58][59] associated with neurotransmitter, GABA, within the mitochondrial oxidative phosphorylation structure (mtOCR). These associations do not identify which is cause or effect; however, previous studies show that increasing calcium concentration directly inhibits mitochondrial oxidative phosphorylation in a dose-and time-dependent manner [67]. Other examples illustrate opposing associations that are not easily predicted. These include changes in an iron sulfur cluster assembly protein, IBA57, required for respiratory chain complexes I and II assembly [68], and cytochrome b5 reductase 3, required for cellular lipid homeostasis, [69,70], both of which negatively associated with lipid metabolites.
Multiple transcript-metabolite associations with oxidative stress (mtOx) have previously been implicated in mitochondrial dysfunction and neuroinflammation as mechanisms of neurological diseases, including Alzheimer's disease [71,72] and progressive multiple sclerosis [73]. Some of the reactions that appear to represent antagonistic interaction may actually represent generalized disruption of cell function. For instance, decreased LIPG and CYP27A1 could reflect a generalized dysregulation of lipid metabolism in response to mtOx. The causal relationship between phenotypes, however, will have to be determined. CYP27A1 encodes a mitochondrial oxidase functioning in cholesterol metabolism. In addition, a decrease in linolenic acid and elevation in CPT1A could indicate a generalized disruption of lipid metabolism. Such an interpretation is consistent with CYP27A1-knockout mice and patients with CYP27A1 loss of function mutations, having hypertriglyceridemia in the brain characterized by behavior and motor dysfunction [74][75][76].
Similarly, opposing signals associated with the oxidation of cellular thiols could reflect generalized disruption of the redox proteome. Thiols are the main determinants of total antioxidant capacity of a cell [77], and previous studies show metals decrease cellular thiols and cause actin cytoskeletal remodeling [78]. The current results show an increase in several cytoskeletal genes, but do not show whether these are specific or non-specific responses to thiol oxidation. These include genes for collapsin response mediator protein 4 (CRMP4), a highly expressed protein in the adult nervous system regulating axonal guidance [79], FER, a cytoplasmic protein tyrosine kinase regulating microtubules [80], and CLIP1, a microtubule protein in vesicular trafficking and associated with intellectual disability [81]. Our findings build upon existent knowledge where similar pleiotropic signaling by reactive oxygen species has been previously described [82,83] While systems biology research has emphasized the need for network approaches to understand complex systems, extensions to define opposing responses within subnetworks are limited. The communities identified in the present study were obtained from systematic variation in a single independent variable, Mn, and future studies will need to address how multiple interacting exposures impact response communities. In the present analysis, only Community one was directly related to Mn, and the other communities were associated with the dependent variables mtOx and mtOCR. Simultaneous exposures to other oxidants, antioxidants or respiratory fuels or inhibitors could alter these response structures. This exposure also includes the environment in which cells grow. Neuronal cells in the brain sense oxygen, leading to changes in its activity [84]. Response structures to Mn and dependent variables may be different under different oxygen-dependent conditions, such as hypoxia. Although not discussed, other associations between communities were present (see Supplemental Tables S1-S9). Additionally, associations below the feature selection threshold |r| > 0.8, were not visible in the data shown, which may constitute proteins encoded by mitochondrial DNA. The central conclusion that antagonistic interaction is a common feature of the complex systems response to excess Mn is supported by a wealth of evidence (Tables S1-S8). Targeted location of these system networks also becomes important as proteins are transferred from cytosol to mitochondria, especially mitochondrial proteins encoded by the nucleus or repair enzymes produced outside the mitochondria. Additional understanding will be derived from future identification of functions of miRNAs and lncRNAs, from experimental evidence in other cell types, integration with other omics platforms such as epigenetics and by establishing the interactions of pathways with Mn in primary versus immortalized cell lines, which will be needed for a global understanding of cellular and mitochondrial responses [85].
A strength of the approach is that the integration of multiple omics platforms can be used to improve mechanistic knowledge for any natural or anthropogenic environmental chemical exposure, drug or therapeutic intervention. The approach can be applied to any cell or animal model, as well as human studies. Recognition that opposing reactions are inherent within response structures aids in defining positive and negative responses, such as those that occur in bi-directional signaling between mitochondria and cell nuclei [86]. Additional research will be needed to understand which of these opposing reactions, whether functionally dependent or independent, are truly beneficial versus harmful responses, as distinguished from those that might represent rectifying pleiotropy, i.e., opposing responses that have evolved to allow homeostatic recovery after initial response to challenge. Such responses are well-known following inflammation, where resolvins are produced as silencers of proinflammatory signaling [87].
In conclusion, a data-driven bioinformatics analysis of mitochondria-cell phenotype in neuroblastoma cells dosed with Mn showed three major network communities, each of which had extensive evidence for opposing response structures, which we interpret as antagonistic interaction. While these can be segregated into individual responses, they likely work in concert to cause Mn-related neurological disorders. The approach using xMWAS to identify these opposing reactions is general and can be broadly applied to model systems and human health research. The implications for biomedical research are that all manifestations of health exposures are likely to include opposing detrimental and beneficial responses, which will need to be addressed at a systems level for effective intervention.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Raw data from three previously published distinct studies [27,33,34] were used.

Conflicts of Interest:
The authors declare no conflict of interest.