Machine characterization of rat toxin responses identifies disease states, tolerance mechanisms and organ to whole-body communication

Background Living organisms are constantly exposed to toxic xenobiotics and have therefore evolved protective responses. In mammals, the liver and kidney play central roles in protecting the organism from xenobiotics, and are at high risk of xenobiotic-induced injury. Liver and kidney damage by drugs and industrial toxins have been extensively studied from both classical histopathologic and biochemical perspectives. Methods and Findings We introduce a machine learning approach for the analysis of toxicological response. Unsupervised characterization of physiological and histological changes in a large toxicogenomic dataset revealed nine discrete toxin-induced disease states. Transcriptome analysis showed that some of the machine-identified disease states correspond to known pathology, and to known effects of certain toxin classes, but others were novel. Analysis of dynamics revealed transitions between disease states at constant toxin exposure, mostly in the direction of decreased pathology, which implies induction of tolerance. Tolerance correlated with induction of known xenobiotic defense genes and novel decreased ferroptosis sensitivity biomarkers. These data reinforce emerging evidence that ferroptosis drives organ pathology, and suggest that its downreagulation may promote tolerance and recovery. Lastly, mechanism of body weight decrease, a known primary marker for toxicity, was investigated. Combined analysis of food consumption, body weight, and molecular biomarkers indicated that organ disease states promote cachexia by whole-body signaling through Gdf15 and Igf1, suggesting strategies for therapeutic intervention that may be broadly relevant to human disease. Conclusions Application of machine learning to systematic data collection of physiology, histopathology, transcriptome reveals multiple disease states, tolerance mechanisms and organ to whole-body communication.

states correspond to known pathology, and to known effects of certain toxin classes, but 23 others were novel. Analysis of dynamics revealed transitions between disease states at 24 constant toxin exposure, mostly in the direction of decreased pathology, which implies 25 induction of tolerance. Tolerance correlated with induction of known xenobiotic defense 26 genes and novel decreased ferroptosis sensitivity biomarkers. These data reinforce 27 emerging evidence that ferroptosis drives organ pathology, and suggest that its 28 downreagulation may promote tolerance and recovery. Lastly, mechanism of body 29 weight decrease, a known primary marker for toxicity, was investigated. Combined 30 analysis of food consumption, body weight, and molecular biomarkers indicated that 31 organ disease states promote cachexia by whole-body signaling through Gdf15 and 32 Igf1, suggesting strategies for therapeutic intervention that may be broadly relevant to 33 human disease. 34 values. Then, we found the enzymes' Enzyme Commission (EC) numbers, which was 286 available in org.Rn.eg.db package. All of the enzyme-encoding genes in the 14 287 pathways were either oxidoreductases (EC1), transferases (EC2), or hydrolases (EC3). 288 They require cofactors for the enzymatic functions. Except for EC3, which requires 289 water as cofactor that is abundant in cells, we regrouped EC1 and EC2 enzymes based

Discovery of biomarkers to different cell death phenotypes 299
We previously showed that 2,565 compounds which were relatively cell-line 300 selectively lethal compounds (i.e., compounds whose GI50 values vary across cell lines) 301 in the NCI-60 dataset were classified into 18 mechanistically distinct classes based on 302 the relative growth inhibitory (GI50) profiles [33]. While most of their mechanisms of 303 action were not yet characterized, a few annotated ones were DNA-targeting 304 compounds (DNA), ferroptosis (Fer), and tyrosine kinase inhibitors (TKI). In the paper, 305 we correlated basal microarray expression profiles with drug sensitivity profiles of each 306 mechanism class, where positive and negative correlations can be interpreted as 307 overexpressed in more resistant or sensitive cell lines, respectively. We extended this 308 Severity scores (Kidney) to their names, DS1-9, that are consistent throughout the paper. See also Figures S1-6. 405

406
We next focused on liver and kidney histology data to ask if they support the 407 discrete disease states that emerged from physiology data. Histology is typically 408 employed after physiology to diagnose human disease. It is also part of the routine 409 exercise required in regulatory toxicity assessments. In Open TG-GATEs, these data 410 were recorded as calls from a standard constrained terminology made by expert 411 toxicologic pathologists based on microscopic examinations of H&E stained tissue 412 sections. One compound treatment can induce more than one histopathology 413 phenotype concomitantly, and typically, a tissue experiencing severe injury exhibits 414 multiple histological phenotypes (Fig. S4-5). We computed 'severity scores' based on 415 the number of abnormal histology phenotypes scored. We then color coded the 416 physiology t-SNE map by liver and kidney histopathology severity scores (Fig. 1C). 417 Visual inspection confirmed good correlation between physiology and histology, where 418 small distinct islands exhibited abnormal liver and kidney histology. However, conditions 419 exhibiting similar physiological abnormality did not always agree on the 420 histopathological severity scores. We suspected that this is in part due to lack of reliable 421 severity metrics in histopathology; histopathology calls made by human experts are 422 intrinsically less quantitative than physiology metrics (e.g., low-level changes may be 423 missed) or may suffer from bias and/or inaccuracy due to scoring only a few tissue 424 slices per animal. To resolve such discrepancy between physiology and histopathology, 425 and we decided to primarily rely on physiology rather than histopathology data where 426 there is a disagreement; by imputing the liver and kidney severity scores on the 427 physiology map, we have a more consistency between the two data and we can set 428 unified criteria for identifying pathological conditions relying on both physiology and 429 pathology datasets were set. (Fig. 1D). 430 We next called disease states by identifying discrete clusters in the conjoined 431 map of physiology and histopathology metrics using a density-based clustering 432 algorithm [18]. We iterated the whole procedure from t-SNE to cluster calling 100 times 433 to compensate for the stochastic nature of the t-SNE algorithm [19]. Nine consensus 434 clusters emerged robustly across the 100 runs and were termed "disease state (DS)" 435 clusters (Fig. 1E). While use of the word "disease" implies a negative health impact, it is 436 in principle possible that some "disease states" are beneficial. Of 3,564 conditions, each 437 DS contained 37 to 203 conditions. 2,723 conditions did not belong to any DSs, and we 438 classified these as non-disease state (non-DS) (Table S1). DSs are listed in Table 1 This table summarizes the analysis of physiology (Fig. 2B), histology (Fig. 2C), 443 transcriptome (Fig. 3), toxin class (Fig. S6), and kinetics (Fig. 4). Physiology acronyms 444 are described in Fig. S2A. Refer to the corresponding figures for the details. 445 446

Known and novel disease states identified by unsupervised analysis 447
We next tried to map DSs onto standard pathologies in the toxicology literature, 448 referring primarily to the physiology and histology markers, and in some cases also to 449 subsequent analysis of gene expression changes or toxin mechanisms (described 450 below). Comparison of histopathological severity scores showed which tissues were 451 more affected in each DS ( Fig. 2A). A simplified clustergram highlighted physiology 452 parameters that most strongly defined each state (Fig. 2B). As we hoped, a distinct set 453 of histology phenotypes were associated with DSs (Fig. 2C)

Overrepresentation of drug classes in disease states 496
We next asked whether specific toxin classes reliably induced specific DSs, 497 defining classes as containing multiple compounds with known overlapping biological 498 activity, as summarized in Figure S1C. We found that lipid-lowering drugs mapped to 499 DS1, synthetic hormones mapped to DS2, and non-steroidal anti-inflammatory drugs 500 (NSAIDs) mapped on DS8, each at higher doses, and each more than expected by 501 chance (Table 1, Fig. S6A). NSAIDs are known to cause intestinal bleeding at high 502 doses [24,25], which likely accounts for their mapping to DS8 (Fig. S6B,C). The four 503 lipid-lowering drugs mapped onto DS1 were peroxisomal proliferator-activated receptor 504 alpha (PPARa) agonists (clofibate, fenofibrate, WY-14643) and a cholesterol synthesis 505 inhibitor (simvastatin). PPARa agonists increased peroxisomes, which were recognized 506 as eosinophilic granules in the cells in DS1 [26] (Fig. S6D,E). This class of drugs has 507 been well studied in rats, and its lipid-lowering effect in the short term is often described 508 as beneficial rather than pathological. However, we decided to keep our initial notation 509 of "disease states" even for DS1; long term treatment with these drugs frequently 510 causes liver cancers, possibly due to hyperactivation of metabolism to increase the liver 511 biomass [27]. The NSAIDs-and fibrates-induced disease states of liver were previously 512 characterized by transcriptome-centric analysis of Open TG-GATEs data [28]. 513 Compared to the two drug classes that strongly induced representative physiology and 514 histology phenotypes of the DSs, the effect of synthetic hormones on liver was modest, 515 and lacked conspicuous physiology or histology changes such as hypertrophy (Fig.  516 S6F,G). 517 518

Transcriptome description of disease states 519
To test if each DS was transcriptionally distinct, we performed elastic net 520 classification of liver and kidney transcriptome data (Fig. S1, Fig. S7A). This generated 521 classifiers that attempted to distinguish conditions assigned to each DS from all the rest 522 of conditions using the liver or kidney transcriptome. All classifiers were found to be 523 powerful for separating DSs (all areas under ROC curves were above 0.85, compared 524 to 0.46-0.63 for randomly drawn samples of the same sizes). Thus, each DS has a 525 characteristic transcriptome both in liver and kidney. To examine the functional 526 implications of DS-specific transcriptional states, we assessed whether 914 GO and 527 KEGG pathways changed their activities in DSs, compared to pooled non-DS clusters. 528 We computed "activity scores" by modifying the standard gene set enrichment analysis 529 (GSEA) to capture the significance of systematic changes of pathway activity among 530 conditions assigned to each DS (See Methods for the detail). We interpreted a large 531 positive or negative activity to indicate that a pathway is substantially up or 532 downregulated compared to corresponding vehicle treatments in the DS. Through this 533 analysis, we found that six DSs (DS1,2,5-8) induced systematic changes of pathway 534 activities in the liver transcriptome; in DS9, changes of pathway activities were observed 535 in the common kidney transcriptome but not in the liver transcriptome; and DS3-4 did 536 not capture any common pathway changes either in liver or kidney (Fig. 3A). 537 Hierarchical clustering of pathway activity in both liver and kidney primarily divided nine 538 DSs into two groups (tissue injury or not) (Fig. 3B); consistently, five DSs associated 539 with tissue injury activated pro-inflammatory cytokine signaling, namely cellular 540 response to proinflammatory cytokines (TNF, IFN-g, IL-1) in the corresponding tissues 541 (Fig. S7B).

551
Liver was more transcriptionally responsive to xenobiotic stimuli than kidney 552 across many conditions that induced a DS. This is expected because the liver is 553 primarily responsible for detoxifying xenobiotics, although the choice of toxins in the 554 dataset may also over-emphasize the role of the liver. The kidney only exhibited major 555 transcriptional changes in DS9, which also phenotypically corresponded to kidney 556 damage in the physiology and histopathology (Table 1). 557 To better understand the changes induced by toxin exposure in the six DSs we 558 identified as significantly affecting the liver transcriptome (DS1,2,5-8), we classified the 559 modulated pathways based on their transcriptional activity. Pathways that were 560 transcriptionally activated or suppressed in each DS were mapped onto corresponding 561 nodes on a dendrogram (Fig. 3C, Fig. S7C-E). This allowed us to determine in which 562 contexts the pathways were up-or down-regulated. For example, in all of these six DSs, 563 the liver activated transcription of "cytosolic large ribosomal subunit"; in disease states 564 that reported on tissue injuries (DS5-8) pathways categorized as "cellular response to 565 IL-1 and TNF" was activated in the liver. On the other hand, some pathways such as 566 "p53 signaling" or "collagen biosynthetic process" were activated only in specific states 567 (e.g. hepatocellular injury (DS7)), but not in the other DSs. 568 569

DS2 is a state of drug-induced tolerance 570
A strength of the Open TG-GATEs data is the collection of multiple time points at 571 constant toxin exposure. Kinetic analysis showed that some DSs are mostly or 572 exclusively late-onset, i.e. after 24hrs (DS1,2,8). Others are early-onset or had no 573 particular trends in kinetics (Fig. 4A, Fig. S8A). Of 365 conditions (compound and dose) 574 whose data were collected at all eight time points, 246 (67%) caused at least one DS at 575 some point, and 90 (25%) caused more than one DS over time (Fig. S8B). for each DS were shown. See also Figure S8. 585

586
To inform on causal connections between DSs we focused on conditions causing 587 more than one DS, and classified temporal transitions at constant toxin exposure (Fig.  588   4B). This analysis revealed bi-directional inter-conversion between different DSs 589 corresponding to liver damage, which may be expected given transcriptional and 590 histologic overlaps between these related DSs. A strong, unexpected feature was 591 conversion of multiple liver and kidney DSs to DS2. Since DS2 is not associated with 592 liver or kidney injury by histology, this suggests a time-dependent reduction of organ-593 specific pathology. Fig 4C shows   Toxin-induced tolerance to toxin action, also known as autoprotection [29,30], is 600 poorly understood but is presumably an important component of how xenobiotic 601 resistance evolved, and how modern vertebrates adapt to toxic environments. Analysis 602 of dynamics of DSs and injury biomarkers (Fig. 4B,C) suggests that DS2 is a state of 603 induced tolerance. To determine molecular mechanisms that might drive tolerance and 604 result in DS2, we first re-examined genes that are selectively regulated in this DS. 605 Xenobiotic catabolism genes were strongly and selectively induced in DS2 compared to 606 all other DS and non-DS conditions (Fig. 3D). While this finding might be expected, it 607 emphasizes the function of xenobiotic defense genes, and supports our characterization 608 of DS2 as a state of tolerance. 609 610

Tolerance is associated with induced resistance to ferroptosis 611
Xenobiotic metabolism is a multi-step reaction: in phase I, cytochrome P450 612 monooxygenases (Cyp450) conjugate xenobiotic compounds with oxygen using 613 NADPH; in phase II, the products of this reaction are conjugated with hydrophilic groups 614 such as sugars (e.g. glucuronic acid) and glutathione (GSH) to facilitate excretion. 615 When we regrouped xenobiotic catabolism genes based on co-factors, we found that 616 genes encoding NADPH-and GSH-utilizing enzymes were among the most highly 617 expressed genes in DS2 (Fig. 5A, Fig. S9A). Activation of redox metabolism functions is 618 also involved in ferroptosis, a form of cell death that has been implicated in drug-619 induced liver injury [31]. Thus, we also suspected that the liver acquires resistance to 620 ferroptosis in DS2. 621  expression signatures that serve as markers for sensitivity (sen) or resistance (res) for 639 each of these treatments. We converted these six signatures to rat orthologs and 640 evaluated them across all nine DS compared to non-disease samples. Fer-res was 641 strongly and uniquely upregulated by DS2 (Fig. 5C). Conversely, DSs associated with 642 liver injuries downregulated Fer-res and upregulated Fer-sen, consistent with the role of 643 ferroptosis in liver injury. DS2 does not increase DNA-res or TKI-res, suggesting that 644 toxin-induced tolerance is associated with acquiring resistance to ferroptosis, but not to 645 other cell death mechanisms. 646 There was no substantial overlap between genes involved in Fer-res and the 14 647 pathways activated exclusively in DS2 (Fig. S9B). Nevertheless, the transcriptional 648 activity of Fer-res among 3,528 liver transcriptome is most highly correlated with the 649 scores of the 10 GO and KEGG pathways exclusively upregulated in DS2 while the 650 activity of Fer-sen is somewhat anticorrelated with them (Fig. S9C,D). Our analysis 651 supports the hypothesis that DS2 not only overexpresses xenobiotic metabolism genes 652 but also acquires resistance to ferroptosis. 653 654

Whole body response to toxins 655
Toxin responses have a whole body impacts in addition to organ-specific effects, 656 but these have been much less studied in the toxicology literature, or the human 657 liver/kidney disease literature. Loss of body weight is a well known indicator of chronic 658 toxicity in rats, though its etiology is poorly understood. It likely corresponds to cachexia 659 phenotypes in man, which are well-recognized contributors to disease burden [34]. 660 Body weight decrease was the most significant physiological descriptor for DS2, but it 661 was not unique to DS2: all DSs associated with tissue injury (DS5-9) exhibited body 662 weight decrease to a similar extent to DS2 (Fig. 6A,B), while tolerance-associated 663 pathways were activated in DS2, but not in DS5-9, highlighting the differences between 664 the states (Fig. 6C). Toxin-induced weight loss could result from direct responses of 665 peripheral tissues such as muscle and adipose tissues to toxin, or (more likely) from 666 their indirect responses to signals from toxin-exposed tissues that alter global 667 metabolism or/and feeding behavior. The liver and kidney communicate with the whole 668 body through concentrations of multiple metabolites in blood, and also through secreted 669 signaling proteins. We therefore analyzed metabolite and secretome changes in 670 conditions that caused body weight loss, irrespective of whether tolerance was induced. 671