Computational model for fetal skeletal defects potentially linked to disruption of retinoic acid signaling

All-trans retinoic acid (ATRA) gradients determine skeletal patterning morphogenesis and can be disrupted by diverse genetic or environmental factors during pregnancy, leading to fetal skeleton defects. Adverse Outcome Pathway (AOP) frameworks for ATRA metabolism, signaling, and homeostasis allow for the development of new approach methods (NAMs) for predictive toxicology with less reliance on animal testing. Here, a data-driven model was constructed to identify chemicals associated with both ATRA pathway bioactivity and prenatal skeletal defects. The phenotype data was culled from ToxRefDB prenatal developmental toxicity studies and produced a list of 363 ToxRefDB chemicals with altered skeletal observations. Defects were classified regionally as cranial, post-cranial axial, appendicular, and other (unspecified) features based on ToxRefDB descriptors. To build a multivariate statistical model, high-throughput screening bioactivity data from >8,070 chemicals in ToxCast/Tox21 across 10 in vitro assays relevant to the retinoid signaling system were evaluated and compared to literature-based candidate reference chemicals in the dataset. There were 48 chemicals identified for effects on both in vivo skeletal defects and in vitro ATRA pathway targets for computational modeling. The list included 28 chemicals with prior evidence of skeletal defects linked to retinoid toxicity and 20 chemicals without prior evidence. The combination of thoracic cage defects and DR5 (direct repeats of 5 nucleotides for RAR/RXR transactivation) disruption was the most frequently occurring phenotypic and target disturbance, respectively. This data model provides valuable AOP elucidation and validates current mechanistic understanding. These findings also shed light on potential avenues for new mechanistic discoveries related to ATRA pathway disruption and associated skeletal dysmorphogenesis due to environmental exposures.

All-trans retinoic acid (ATRA) gradients determine skeletal patterning morphogenesis and can be disrupted by diverse genetic or environmental factors during pregnancy, leading to fetal skeleton defects. Adverse Outcome Pathway (AOP) frameworks for ATRA metabolism, signaling, and homeostasis allow for the development of new approach methods (NAMs) for predictive toxicology with less reliance on animal testing. Here, a datadriven model was constructed to identify chemicals associated with both ATRA pathway bioactivity and prenatal skeletal defects. The phenotype data was culled from ToxRefDB prenatal developmental toxicity studies and produced a list of 363 ToxRefDB chemicals with altered skeletal observations. Defects were classified regionally as cranial, post-cranial axial, appendicular, and other (unspecified) features based on ToxRefDB descriptors. To build a multivariate statistical model, high-throughput screening bioactivity data from >8,070 chemicals in ToxCast/Tox21 across 10 in vitro assays relevant to the retinoid signaling system were evaluated and compared to literaturebased candidate reference chemicals in the dataset. There were 48 chemicals identified for effects on both in vivo skeletal defects and in vitro ATRA pathway targets for computational modeling. The list included 28 chemicals with prior evidence of skeletal defects linked to retinoid toxicity and 20 chemicals without prior evidence. The combination of thoracic cage defects and DR5 (direct repeats of 5 nucleotides for RAR/RXR transactivation) disruption was the most frequently occurring phenotypic and target disturbance, respectively. This data model provides valuable AOP elucidation and validates current mechanistic understanding. These findings also shed light on potential avenues for new mechanistic discoveries related to ATRA pathway disruption and associated skeletal dysmorphogenesis due to environmental exposures.

OPEN ACCESS EDITED BY
Patricia Ruiz, Agency for Toxic Substances and Disease Registry, United States Introduction Retinoid signaling plays an important role in the patterning, differentiation and homeostasis of the developing skeletal system (Janesick et al., 2015;Ghyselinck and Duester, 2019;Mezquita and Mezquita, 2019;Knudsen et al., 2021). All-trans retinoic acid (ATRA), the biologically active form of vitamin A, is an endogenous mediator of skeletal patterning and differentiation as evidenced by various animal models of retinoid depletion or excess during pregnancy (Shannon et al., 2017: Cunningham andDuester, 2015). Sensitive gestational stages correspond to early specification of the body plan during gastrulation and subsequent organogenesis. Specifically, ATRA signaling influences spatial patterning of major body axes (e.g., anterior-posterior, dorsalventral), cranio-facial development, segmentation of the vertebral column, early limb outgrowth and skeletal development (Draut et al., 2019;Knudsen et al., 2021). ATRA deficiency and excess (e.g., via dietary retinol deficiency, exposure to excess ATRA or retinoid compounds, or functional inactivation of key nodes in retinol metabolism and/or signaling in mouse mutant models of ATRA deficiency or excess) have been shown to cause developmental skeletal defects (Cunningham and Duester, 2015;Roberts, 2020;See et al., 2008).
ATRA is the best known endogenous active metabolite of (retinol → RDH → RALDH → ATRA) vitamin A and is considered the cognate ligand for the retinoic acid receptor (RAR). Once inside the cell, ATRA's signal is transduced by specific nuclear hormone receptor complexes consisting of RAR (RARA, RARB, RARG) and retinoid X receptor (RXRA, RXRB, RXRG) heterodimers, which act, in general, as ligand-activated transcription factors at a retinoic acid response element (RARE) (Mark et al., 2006;Niederreither and Dollé, 2008;Blaner et al., 2016). This complex binds DNA and induces changes in gene expression (Shannon et al., 2020) during embryonic development (Cunningham et al., 2013;Shimozono et al., 2013;Chawla et al., 2016;Schubert and Gibert, 2020). Several isoforms of each RAR and RXR exist, with distinct spatial temporal expression patterns. The diverse effects of ATRA/retinoids during normal embryonic development are mediated by the various RAR/RXR heterodimer combinations of each of these isoforms (Chambon, 1994;Kastner et al., 1994;Mark et al., 2009). The classical direct repeat (DR) for RAR/RXR binding has a 5-nucleotide spaced DR (referred to as DR5) (Balmer and Blomhoff, 2005).
In the Organisation for Economic Co-operation and Development's (OECD's) 2012 Detailed Review Paper (DRP) 178 focused on Endocrine Disruptor Screening Program (EDSP), the retinoic acid signaling pathway was ranked second (below PPAR signaling) among seven pathways considered to be susceptible to environmental endocrine disruption and for which relevant endpoints could be measured in a chemical testing battery (OECD 2012(OECD , 2014. Some of the most prominent adverse effects of systemic retinoid exposure, such as fetal skeletal defects, were the target for an OECD workplan to identify chemical hazards to development based on relevant in vitro assays (OECD, 2021).
Conserved cell signaling through ATRA-dependent gene expression has been well documented and shown to have developmental effects on most tissues (Yousefi and Azizzadeh, 2010). Prenatal development is particularly vulnerable to genetic, pharmacological, or chemical disruption of the retinoid pathway during gastrulation-organogenesis (Knudsen et al., 2021;OECD, 2021). This is especially the case for early gestation when the fundamental body plan is established, and subsequent stages when the regional pattern of specific body segments is decoded. At least 12 assays in the ToxCast/Tox21 portfolio map to molecular targets in the retinoid signaling pathway. A preliminary analysis revealed low-/submicromolar bioactivity on one or more target assays for over 100 structurally diverse ToxCast/Tox21 chemicals (e.g., conazoles, organochlorine pesticides, organotins, retinoids, and pharma compounds) suggesting that they can be used to generate models of the retinoid system and provide predictive toxicological information relevant to developmental disruption (Baker et al., 2018).
Skeletal defects are among the most prevalent adverse fetal outcomes associated with prenatal developmental toxicity in EPA's ToxRefDB database (Knudsen et al., 2009). Menegola et al. (2021) presented an adverse outcome pathway (AOP) on disruption of ATRA signaling pathway leading to craniofacial defects. Consequences of environmental disruption of ATRA signaling pathway can lead to stage-and region-specific deficiencies for various skeletal elements (Shenefelt, 1972;Williams and Bohnsack, 2019;Zhang et al., 2015;Qin et al., 2014;Lee et al., 2012). For example, visceral and skeletal anomalies have been observed in various animal models (e.g., rat, mouse, rabbit, chick) and specifically, retinoic acid administration at embryonic day E9.5 led to hypoplasia of the branchial arches, as well as auricular and eye anomalies in mice (Glineur et al., 1999).
With tens of thousands of chemicals currently in commerce or in the environment, New Approach Methods (NAMs) including high-throughput screening/high content screening assays and computational/in silico models aim to identify biological pathways and chemically induced biological activity in human cells. With an emphasis on protecting susceptible populations and lifestages such as embryonic development, NAMs can develop predictive models of in vivo biological response that would ignite a shift from traditional animal endpoint-based testing to human pathway-based risk assessment and hazard identification, (Collins et al., 2008;US EPA, 2021). The Lautenberg Chemical Safety for the 21st Century Act (15 USC 2601 Public Law 114-182) emphasizes the need for reduction and replacement of the use of vertebrate animals in toxicity testing paradigms (2016). To this end, the ToxCast and Tox21 projects are building high-throughput screening and high content screening in vitro datasets, while the Toxicity Reference Database (ToxRefDB; Watford et al., 2019) is a rich database of in vivo data for anchoring predictive models. ToxCast and Tox21 assays cover a wide array of molecular/cellular signaling pathways Thomas et al., 2019), and human stem cell line-based biomarker assays are predictive for developmental toxicity (Zurlinden et al., 2020). ToxRefDB provides prenatal developmental toxicity studies in pregnant animals (e.g., rats, mice, rabbits, and other mammals) (Knudsen et al., 2009;Sipes et al., 2011). Linking the data from all three sources supports statistical models predicting potential for adverse effects (e.g., developmental), as well as identifies potential molecular targets and cellular pathways for incorporation into virtual tissue models simulating cellular dynamics (Knudsen and Kleinstreuer, 2011;Kleinstreuer et al., 2014;Leung et al., 2016;Hutson et al., 2017). Models developed using human-relevant mechanistic information build confidence in silico models' scientific reliability and relevance for regulatory decisions. A detailed predictive computational (in silico) signature model is therefore necessary to understand both normal embryonic skeletal development and how environmental factors may lead to a variety of skeletal developmental defects. Herein, we integrate multiple databases to develop in silico models that identify chemicals associated with ATRA signaling pathway disruption and fetal skeletal defects.
We used the ToxCast/Tox21 in vitro data with ToxRefDB in vivo studies in pregnant vertebrates to build a novel predictive model to identify potential developmental toxicity of chemical compounds associated with ATRA signaling pathway disruptions (RARA/B/G, RXRA/B, DR5, CYP26 Surrogates-onwards labelled CYP surrogate biomarker) and a range of skeletal defects (e.g., axial, appendicular, and cranial defects). Here, we hypothesize that a predictive model of chemical compounds associated with embryonic skeletal developmental toxicity and ATRA signaling pathway disruption will reliably provide scientifically based principles for regulatory decisions regarding chemical use related to pregnant vertebrates. This study uses NAMs data to: 1) compile data available in vitro ToxCast and Tox21 High-Throughput Screening assays on ATRA signaling and metabolism for n = 374 chemicals; 2) cull information from ToxRefDB and select animal studies on prenatal developmental toxicity for n = 370 ToxCast chemicals with skeletal defects identified; 3) identify the relationships of these chemicals and their association with skeletal defects using literature mining (Baker et al., 2018); 4) systematically organize the in vitro and in vivo findings to provide insight into potential molecular initiating events (MIEs) on the ATRA pathway that may lead to testable AOPs. While providing foundational weight-of-evidence, a lack of three-dimensional, dynamic biological systems limits the applicability of AOPs to systemic problems, ushering in the need for more encompassing NAMs.

Phenotypic data compilation
Using Python v.3.3, the Toxicity Reference Database (Watford et al., 2019; ToxRefDB v.2, https://github.com/ USEPA/CompTox-ToxRefDB, accessed November 2020) was mined for chemicals associated with developmental skeletal defects across 2,946 prenatal developmental toxicity studies. Python was used to categorize endpoint targets anatomically by available ToxRefDB skeletal defect annotations of Limb, Bone, Paw/Digits, and Mouth/Jaw. This resulted in the retrieval of 57,198 annotated features linked to skeletal defects ( Figure 1). The distribution of skeletal defects by species included 31,661 in rats, 1,232 in mice, 16,375 in rabbits, and 368 in chinchillas. This was a composite from 363 ToxRefDB chemicals. ToxRefDB Frontiers in Pharmacology frontiersin.org 03 includes bone elements entered as individual targets (e.g., ribs, vertebrae, scapula, sternum, ulna, etc.) for 206 fetal bones. The database also includes the elementary description of the type of defect for each bone (e.g., absent, bent, incomplete ossification, misshapen, etc.) (Knudsen et al., 2009). Based on individual targets, data was parsed into developmental skeletal phenotypes by one of three major anatomical regions-appendicular, axial, or cranial. Anatomical development of skeletal regions has been well described for appendicular and axial skeletal regions; here, axial defects are separated into cranial and post-cranial (trunk) regions. Appendicular defects (representative of 126 bones) were classified into stylopod, zeugopod, or autopod -segments of the fore-and hindlimb skeleton. Effects on the cranial skeleton (28 bones) were classified as neurocranial, orofacial, or viscerocranial regions, and the axial defects (representative of 52 bones) were classified into thoracic cage (ribs, sternum), vertebral, or cauda regions (Knudsen et al., 2021). Lastly, a category named "other" was allocated for 7,562 unspecified developmental skeletal defects. Unspecified developmental skeletal defects are representative of data entries with elementary description of the type of skeletal defect (e.g. absent, misshapen, unossified, etc.), but that do not indicate the individual target bone (Watford et al., 2019). The inclusion of unspecified data may limit the precision of the model but Chemicals associated with 57,198 composite skeletal defects were extracted from in vivo ToxRefDB and ToxCast Databases; defects were categorized into ten skeletal regions (blue boxes). Chemicals that induced in vitro bioactivity in ATRA pathways in 1 or more of 10 relevant ToxCast/ Tox21 assays (red box). In vitro and in vivo data was compared, distinguishing 48 chemicals associated with both skeletal defects and bioactivity for relevant ATRA pathway nodes (purple box and Venn diagram). These chemicals potentially elucidate the role of retinoid signaling in skeletal development for AOPs and predictive toxicology applications.
Frontiers in Pharmacology frontiersin.org 04 provides a more complete dataset. Zurlinden et al. (2020) highlighted 42 ToxCast chemicals that are extensively used by researchers in developmental toxicology non-animal models. A few of these well-annotated developmental toxicants were found in ToxCast but not included in ToxRefDB. We searched PubMed using the Abstract Sifter semi-automated literature mining tool  to derive information from scientific publications related to these 42 chemicals, and identified skeletal defects associated with 7 of these substances. Combined, ToxRefDB and ToxCast chemicals linked to developmental defects in 10 skeletal regions as provided in Figure 1.

Bioactivity data compilation
In the U.S. Environmental Protection Agency's ToxCast high-throughput screening program over 4,000 chemicals have been tested for bioactivity for over 1,400 molecular targets/features accessed through InVitroDB v.3.1 at https://epa.figshare.com/articles/dataset/Previously_ Published_ToxCast_Data/6062551 in November 2020 Richard et al., 2016). The breadth of ToxCast assays include biochemical assays, human cells, mouse embryonic stem cells and zebrafish development platforms (Thomas et al., 2019;Parish et al., 2020;Marty et al., 2022). In the Tox21 project, about 8,500 chemicals were screened, including the over 4,000 ToxCast chemicals. Our project mined a combined 10 assays (Table 1) from ToxCast and Tox21 (Chen et al., 2016) for chemicals that demonstrated bioactivity for relevant ATRA nodes (RARs, DR5, RXRs, and CYP surrogate biomarker). CYP1A1 and CYP2C8 were selected as CYP surrogate biomarker, due to the absence of CYP26 assays in ToxCast and Tox 21 (Baker et al., 2022). Importantly, both CYP surrogate biomarker' assays are biochemical. The binding region of CYP26 is similar to CYP2C8. When CYP2C8 inhibitors were tested against CYP26A1, there was a statistically significant correlation between CYP26A2 and CYP2C8 IC 50s (Foti et al., 2016a;Foti et al., 2016b). The enzyme CYP1A1 is also recorded as metabolizing retinoic acid (Lampen et al., 2000). NovaScreen human unspecified cytochrome P450, family 1, subfamily A, polypeptide 1 assay (NVS_ADME_hCYP1A1) and NovaScreen human unspecified cytochrome P450, family 2, subfamily C, polypeptide 8 assay (NVS_ADME_hCYP2C8) assay data was combined as CYP surrogate biomarker, using the lowest AC 50 if redundant. NovaScreen human unspecified retinoic acid receptor, alpha assay (NVS_NR_hRARa_ Agonist) data was combined with Attagene human HepG2 retinoic acid receptor, alpha assay (ATG_RARa_ TRANS_up) information (no overlap in data occurrences) to compose RARA bioactivity figures utilized (Table 1). While these assays have different dynamic range, both test for increases in RARA bioactivity, and the range of concentrations for which activity were selected were limited by a threshold (<10 µM) chosen for both environmental relevance and specificity to avoid generalized cell stress responses . Along the ATRA signaling pathway, in both databases, relevant assays were measured for interactions between chemicals and receptors and enzymes as molecular targets or assessed downstream effects on reporter gene activity. Chemicals were selected from these assays based upon their potency, half-maximal activity concentrations (AC 50s ), and efficacy in one or more related assays. These chemicals were then compared to the list of ToxRefDB/ ToxCast chemicals associated with skeletal phenotypic defects.

Data analysis and visualization
In addition to the multi-database comparison to identify chemicals associated with skeletal defects as well as ATRA signaling pathway disturbances, ToxCast and ToxRefDB chemicals associated with developmental skeletal defects were prioritized for potency and regional effects through the Toxicological Priority Index (ToxPi ™ ) (ToxPi; Marvel et al., 2018). ToxPi is a software that objectively conducts chemical prioritization via integration across multiple information domains and sources of evidence The "ToxPi visualization" portrays the relative magnitudes of environmental hazards (e.g. exposure of environmental chemicals). Data is transformed into transparent, visual rankings (Reif et al., 2010). A ToxPi for chemicals associated with ATRA pathway disturbances was formulated for further comparison between the two (in vivo vs. in vitro) endpoints ( Figure 3). Separately, for each chemical, a ToxPi was created that was broken into 10 slices, each slice representative of the 10 potential regional skeletal defects ( Figure 2). Separately, Slices are proportionate to the number of each effect caused by a particular compound, which the overall potential ToxPi value equaling one. The width of each slice signifies the number of data points per defect. The larger the radius of the slice from the origin, the greater the chemical's potency, with length representing the overall degree of hazard (National Academy of Sciences, 2014). This approach allows for a multidimensional analysis of the relative potency (in vivo experiments), chemical properties (bioavailability), and perturbation score (pathways). A dimensionless index score, ToxPi, is calculated based on the combination of all the data sources, and then compounds are ranked to prioritize which chemicals should be tested for future toxicological testing, compare toxicity, and identify similarities in predicted compound activity.
For k-means clustering, ToxPi employs a Java port of the R function for means. The k-means clustering module uses agglomerative clustering to organize ToxPi profiles into Frontiers in Pharmacology frontiersin.org 05 clusters (collection of data points) based upon similarity and discover underlying patterns. In ToxPi, k-means clustering uses the centroid to parse the chemicals rather than default sorting by overall priority score (rank). In short, k-means clustering identifies k number of centroids, assigns each data point to the nearest cluster, simultaneously keeping the centroids as small as possible. The major sources of organization were then visualized on a coordinate field by Principal Components Analysis (PCA) (Marvel et al., 2018). This function uses the "Hartigan-Wong" algorithm to execute the clustering. For all k-means clustering procedures, the algorithm was run (nStart) at the default of 15 starting cluster locations, as consistency of results was visible through such a procedure. Effective k-means clustering results in smaller within group sum of squares, which can be achieved by sound selection of the seed of the random number generator, and specified replication (Marvel et al., 2018, Manual v. 2.3). The k-means algorithm grouped similar ToxPi data points together based on distance from centroids, the center of these clusters. These groupings identified underlying patterns. For k-means clustering chemicals involving ATRA signaling pathway activated genes, 4 centroids were designated based upon biological reasoning that the number of clusters needed is one less than the total number of the relevant molecular targets involved in the assays. This allows for the statistical subtraction of one degree of freedom; similarly, 4 clusters were employed for k-means analysis of chemicals associated with ATRA pathway target ( Figure 3). For the skeletal defects ToxPi, data points were assigned to one of 5 clusters (Figure 2). Homogeneous subgroups within the data are identified such that data points in each cluster are as similar as possible according to a similarity measure of skeletal region defects or ATRA signaling pathway node disturbance for the respective k-means cluster.
The Hierarchical Clustering module provides options for organizing ToxPi profiles into clusters based upon similarity, rather than the default sorting by overall priority score (rank). The cluster dendrograms are drawn using an agglomerative hierarchical clustering technique. With ToxPi profiles for individual chemicals, data points are treated as individual clusters. Then with additional iterations, similar clusters merge, until there is one cluster formed of similar data point (Reif et al., 2010;Grimm et al., 2016;Marvel et al., 2018). The Ward.D2 hierarchical clustering tool was utilized to minimize variance and visualized with the ToxPi 2.3 software. The chemicals identified as associated with skeletal defects and ATRA pathway disturbances were selected and displayed in a heatmap using the chemical AC 50 per each of the 8 targets (pheatmap function in R).

Initial workflow
ToxRefDB in vivo and ToxCast/Tox21 in vitro data was compared, identifying 48 chemicals ( Table 2) that were associated with both skeletal defects and bioactivity (<10 µM) for relevant ATRA pathway nodes (Figure 1), purple box and Venn diagram). These chemicals have the potential to elucidate the role of retinoid signaling in skeletal development for hypothesizing and testing AOPs and supporting predictive toxicology applications. The 48 chemicals can also be found in the accompanying Abstract Sifter tool organized by target on the AbstractSifter_Retinoid sheet (Supplementary Material S1), in a simple list on the Notes sheet showing the chemical-specific association with nine skeletal defects phenotypes, and on the Landscape sheet, shown with sample toxicity-related queries and resulting article counts. The Abstract Sifter's Landscape sheet results show that, of the 48 chemicals, retinol and ATRA are associated with the highest number of articles (293 and 247, respectively) describing skeletal defects-specific toxicity (S1). Forty-one of the 48 chemicals on the Landscape sheet have some connections to citations about skeletal defectsspecific phenotypes. Seven chemicals (bromuconazole, buprofezin, fenpyroximate (Z,E), linuron, N-phenyl-1,4-benzenediamine, tebufenpyrad and tributyltetradecylphosphonium chloride) do not have literature associating them with skeletal defects. Frontiers in Pharmacology frontiersin.org 06

FIGURE 2
ToxPi Rank Order (A) and Hierarchical Clustering (B) results for 48 chemicals associated with adverse skeletal outcomes (ASOs). ASOs include autopod, stylopod, zeugopod, cauda, thoracic cage, vertebra, orofacial, neurocranium, viscerocranium, and other. Most potent chemicals associated with the greatest amount of target disruption have the largest ToxPis and highest ToxPi rank (48). Skeletal defects were categorized into three skeletal phenotypes and specified regions and an unspecified region as delineated in the ToxPi key (A). Appendicular regions are shades of (Continued ) Frontiers in Pharmacology frontiersin.org 07 ToxPi, K-means clustering, and hierarchical clustering K-means clustering of skeletal defects provided insight into the most sensitive phenotypic changes following chemical exposure. Analysis of the skeletal defects K-means clustering plot ( Figure 4A) indicates that Principal Component 1 (PC1), the x-axis, identifies variability in the 48 model chemicals, denoting the heterogeneity of the skeletal response. Summarily, the left side of Figure 4A PC1 demonstrates less specificity with 8 phenotypic regions (affected by, e.g., retinol, aspirin) represented, while the right side has more specificity as indicated by the individual thoracic cage slice (affected by, e.g., triflumizole). PC2, y-axis, provides the specificity of the regional phenotype based on chemical. Through PC2, a trend in axial phenotypic changes (thoracic cage, cauda, and then vertebral) occur first, followed by an increase in cranial defects (progressively neurocranium, viscerocranium, and then orofacial regions). Thoracic cage alterations occur without other phenotypic changes preceding its own. Again, looking at the systemic implications of these diagrams, thoracic cage formation is tied to homeobox patterning, hence, alterations in patterning may be leading to an effect on the cauda of the embryo. Among appendicular phenotypes specificity is ordered as stylopod, zeugopod, and then autopod. Interpreting the k-means plot for systemic developmental changes there is a visible progression to appendicular responses. The limb defects increase in frequency following cranial and postcranial axial phenotypic changes, indicating that the limb is more responsive when preceded by other phenotypic regions' alterations. In examining the sensitivity of phenotypic regions, through the fipronil example it is evident that while limb defects have less specificity than cranial responses, appendicular phenotypes are more sensitive than cranial phenotypes when exposed to chemicals inducing any adverse appendicular outcomes. Notably, the retinoids of ATRA and retinol are activating through the biological system throughout the embryo with broad effects in a manner consistent with their status as positive controls.
Hierarchical clustering of skeletal defects highlighted the prevalence and type of phenotypic changes. In the hierarchical clustering of skeletal phenotypes ( Figure 2B), the potency (lower concentrations of chemical compound results in AC 50s more potent than other chemicals at the same relative concentration) of retinol and ATRA are greatest in association with skeletal defects in every zone of interest. These findings reinforce the validity of the chemical compound selection process applied, since ATRA and retinol, two known positive controls on the ATRA pathway, were objectively identified via the workflow in Figure 1. The progression of skeletal defects in the hierarchical clustering is consistent with the parallel skeletal region k-means clustering. All clusters contain one or more axial phenotypes; this suggests that axial phenotypes are the most frequent fetal outcomes and possess the greatest specificity of the skeletal defects associated with the 48 chemicals. Clusters a, b, and c are predominantly thoracic cage defects, with eventual viscerocranial and unspecified defects occuring. The consistency of the unspecified defects fit the patterning associated with axial defects, suggesting they may be axial features. Cluster d starts with the consistent occurrence of defects of the zeugopod phenotype, all axial zones, and neurocranial and viscerocranial regions. This hierarchical cluster demonstrated the specific progression of skeletal defects as the sensitive thoracic cage, progression to appendicular (stylopod and zeugopod) and increasing occurrences of cauda and viscerocranium defects.
K-means clustering of targets on the ATRA signaling pathway identified target disruption related to the 48 chemicals of interest ( Figure 4B). In PC1 moving right to left on the x-axis, the first quintile included those chemicals that activated RAR, RXR, and/or DR5 readouts only. The next quintile adds chemicals that also activated the retinol signaling pathway (RSP) assay (Tox21) node. Furthermore, PC1 and PC2 clearly delineate CYP surrogate biomarker, FIGURE 2 red, axial regions are shades of green, cranial regions are blues, and unspecified skeletal areas are purple. Skeletal Defects Hierarchical Clustering for 48 Chemicals: 'Clusters a-e' represent varying skeletal phenotype alterations associated with up to 48 chemicals of interest. Clusters with chemical ToxPis associated with similar skeletal defect composition. 'Cluster a was primarily thoracic cage defects with occurrences of neurocranium and appendicular defects. A ToxPi rank of 1 is least potent, while a score of 48 is a highly potent associated with significant skeletal defects. Rank order of 'Cluster a' chemicals is 1. Fenpyroximate (Z,E), 2. Chlorothalonil, 3. Forchlorfenuron, 4. Imazalil, 5. Lindane, 6. Myclobutanil,8. Tebufenpyrad,9. Tetraconazole,10. Triflumizole,11. Raloxifene hydrochloride,12. Propiconazole,13. Difenoconazole,14. Flusilazole,15. Fluoxastrobin,and 17. Fipronil (47) and ATRA (48) were strongly associated with multiple appendicular, axial, and cranial phenotypic defects.
Frontiers in Pharmacology frontiersin.org 08 Frontiers in Pharmacology frontiersin.org 09 RARA, and DR5 as the most disrupted targets along the ATRA signaling pathway ( Figure 4B). DR5 is the most sensitive target, with disruption occuring in all k-means quadrants to some degree ( Figure 4B). PC2 provides potency ranking with chemicals having the least potency, but greatest specificity at the top of the y-axis; the greatest potency and broad range of abnormalities occurs moving down the y-axis (e.g., retinol, ATRA, etc.). Figure 3B's hierarchical clustering demonstrates that CYP surrogate biomarker and DR5 are the most specific targets to be activated on the ATRA signaling pathway. The potency of chemicals associated with CYP surrogate biomarker appear to be most potent due to the greater degree of activation of FIGURE 3 and highest ToxPi rank (48). A ToxPi rank of 1 is least potent, while a score of 48 is a highly potent associated ATRA bioactivity. Hierarchical Clustering for 48 Chemicals In Vitro Activity on ATRA Signaling Pathway Targets: 'Clusters a-e' represent varying targets disrupted on the ATRA signaling pathway by up to 48 chemicals of interest. Clusters of chemical ToxPis associated with similar target disruption composition. 'Cluster a' is dominated by DR5 and RARA disruption. Rank order of 'Cluster a' chemicals 1.   K-means clustering of Skeletal Defects (A) and ATRA Bioactivity (B) for 48 chemicals. In image A, skeletal defects examined include autopod, stylopod, zeugopod, cauda, thoracic cage, vertebra, orofacial, neurocranium, viscerocranium, and other. Principal Component 1 (PC1), the x-axis, identifies variability in the 48 model chemicals, indicating heterogeneity in the skeletal response. Moving from left to right, a broader range of phenotypic changes are found, and then the array of skeletal responses diminishes with horizontal progression. The left side PC1 demonstrates less specificity, while the right side has more specificity (i.e., single thoracic cage slice). This is demonstrable by the leftmost cluster group 2 having the (Continued ) Frontiers in Pharmacology frontiersin.org 11 CYP surrogate biomarker when compared to even DR5. The overall severity of CYP surrogate biomarker' disruption, common occurrence of DR5 and RARA in the hierarchical clustering are consistent with the k-means clustering. Of particular interest in both the hierarchical and k-means clusterings there are chemical classification patterns. Organochlorine pesticides (e.g., endosulfan, endrin, lindane, pyridaben, etc) consistently activate DR5. Furthermore, azole fungicides (e.g., triazoles, imidazoles, pyrazoles) have broad association with activation of CYP surrogate biomarker, RARA, and DR5. Organotins biocides preferentially activated RXRs (e.g., triphenyltin, tributyltin hydroxide). Sulfurons (e.g., linuron, forchlorfenuron) consistently activate CYP surrogate biomarker. Chemicals with unspecified or broad categorizations (e.g., SAR 150640, asprin, iprione, etc.) regularly activated DR5, followed by CYP surrogate biomarker, then RARA, and activated the other molecular targets of interest with lower specificity.

Literature & Datamining as basis for MIEs and putative AOPs (pAOPs)
To further support the formation of pAOPs, publicly available databases that describe in vitro assay results of a chemical's activity against one of the ATRA pathway targets of interest were mined (Baker et al., 2022). Databases mined included Protein Data Bank (PDB) (Berman et al., 2003), ChEMBL (Gaulton et al., 2017), and ToxCast/Tox21 Richard et al., 2016). Chemicals with inactivity in the ToxCast and Tox21 assays were excluded from the developed ATRA pathway candidate chemicals (Baker et al., 2022) selected from assays and literature. There were 1188 candidate chemicals (Baker et al., 2022) with activity in one or more ToxCast or Tox21 assays of interest. Furthermore, a multi-rule filter was set to limit chemicals with the following qualifying features: DR5 < 20 μM, CYP <10 μM, RAR <10 μM, RXR <10 μM, RSP (Tox21) < 10 μM, and prevalent RDH10 literature. This resulted in the shrinkage of data for retinoid disruption (potential MIEs) and adverse skeletal outcome (ASO) to 117 chemicals of particular interest for pAOP elucidation.
Supplementary Figure S2 provides k-means clustering ATRA Signaling Pathway for 117 chemicals. Figure 5 is the heatmap dendrogram of 48 chemicals associated with ATRA pathway disruption and skeletal defects in ToxRefDB and ToxCast/Tox21, while S2 is a dendrogram of 117 chemicals with ATRA pathway disruption with potential AOP elucidation. For pAOPs, MIEs were derived from the measured target activity increases (agonist assays) and decreases (antagonist assays) in the Tox21 and ToxCast associated with the chemicals of interest. Adverse skeletal outcomes (ASOs) were based on the phenotypic defects recorded in ToxRefDB. The literature review conducted provided quantitative and qualitative data to inform key event (KE) placement in pAOPs. Table 3 displays 3 pAOPs for ASOs.

Discussion
Forty-eight chemicals were found to represent a subset of the chemical landscape having in vitro (ATRA pathway targets) and in vivo (phenotypic skeletal defects) outcomes defined from ToxCast, Tox21, and ToxRefDB. K-means clustering, hierarchical clustering, and the heatmap results demonstrate the same conclusions about the 48 chemicals associated to both skeletal defects and target disruption on the ATRA signaling pathway. Thoracic cage was the most sensitive and frequently occurring skeletal defect in this model, followed by other axial defects (vertebra and cauda), and cranial defects. Simultaneously, DR5 has the greatest occurrence of target activity; chemicals disrupting DR5 were consistently associated with thoracic cage defects. Thereafter, CYP surrogate biomarker and RARA are the most sensitive targets to disruption on the ATRA signaling pathway. The most sensitive skeletal defects and disrupted targets outcomes are often associated with the same chemicals. The most potent chemicals are similar for both skeletal defects and target disruption (e.g., ATRA, retinol, flusilazole, N,N-dimethylformamide), increasing the likelihood that when disruption of a target on the ATRA signaling pathway results in an adverse outcome (AO) occurring with an associated ASO.
It should be noted that the workflow used in developing this model (Figure 1) was performed in series (supervised → unsupervised learning). Supervised learning culled specific defects (skeletal) and assay targets (ATRA system) from the datasets to localize the overlapping domain. From there, we used unsupervised methodology to classify statistical correlations. most in vivo ToxPi regional phenotype slices at upwards of 8 slices (e.g., retinol, aspirin, etc.), while the farthest right cluster group 3 has solely 1 slice of specified regional phenotypic change (e.g., triflumizole). PC2, y-axis, demonstrates phenotypic separation as deficiencies transition from autopod to zeugopod and upward phenotypic changes are then found in the axial cauda, and thereafter moving downward there is an increase in neurocranial responses. Moving further left and upward along PC2 viscerocranial defects increase. Through PC2, a trend in axial phenotypic changes occur first, followed by an increase in cranial defects. In image B, k-means clustering ATRA Signaling Pathway for 48 chemicals: Disrupted activity in the following targets are of interest: RARA/B/C, RXR A/B, DR5, RSP (Tox21), CYP surrogate biomarker. PC1 and PC2 consistently demonstrate that DR5 is the most sensitive target, with interruption occuring in all k-means quadrants to some degree. Overall CYP surrogate biomarker, RARA, and DR5 are most disrupted along the ATRA signaling pathway in the 10 assays examined.

Frontiers in Pharmacology frontiersin.org
This led to inferences using only input vectors without referring to known, or labelled, outcomes for predictive toxicology (Ciallella et al., 2022); however, it neglects non-skeletal targets as well as chemicals that did not perturb in vitro bioactivity profiles of the ATRA assays available for this analysis, or non-ATRA pathways, which are important for health-protective inferences for DART testing (Rajagopal et al., 2022). The k-means clustering (Figure 4) findings for the 48 chemicals were highly relevant, and consistent with corresponding heatmap findings ( Figure 5). However, for the 117 chemicals, the k-means clustering diagram is biased by the weight of the increased occurrence of Tox21 results, and lack of parallel ToxCast data (S2). This is evident with the straight-lined clustering of Tox21 data; while k-means clustering may not provide sufficient data about the relationships between target disruption on the ATRA signaling pathway, the 117 chemical heatmap (S2) findings are consistent with the 48 chemical heatmap ( Figure 5). In both heatmaps Figure 5 and S2, DR5, CYP surrogate biomarker and RARA are highly perturbed by the chemicals of interest, however there is a notable increase in the RSP activity, likely due to the influx of Tox21 data, that was not complimented by equitable ToxCast data. The Tox21 RSP assay is especially related to RARA, and thus remains consistent with other data. These findings allude to the reality that an increase in data though increased Tox21 data may strengthen the conclusions of our findings. However, an increase in Tox21 assay types and number of chemicals in experimentation will lead to more data, require a complimentary increase in ToxCast data to achieve balanced and correctly scaled ToxPi results. Twenty-eight of the 48 chemicals found both in ToxRefDB and Tox21 and/or ToxCast were candidate reference chemicals (Baker et al., 2022), but had not previously been noted. This finding provides reason for further exploration of the impact of the other 20 chemicals on the ATRA signaling pathway and skeletal defects. All 20 compounds were found to be associated with skeletal defects in literature (S1), but not ATRA signaling pathway disruption. For example, prenatal aspirin exposure in rats induces skeletal anomalies including fused ribs, incomplete ossification of the cervical arch, absent body of thoracic or lumbar vertebra, deformation of lumbar arch, and cartilaginous changes in vertebrae, paw and ribs (Dodo et al., 2009). There are 69 chemicals from the 117 chemicals that are lacking ToxRefDB data. This study focused on in vivo skeletal defects data from ToxRefDB; extending in vivo data to include data from other in vivo databases may show greater consistency between the 48 chemicals and 117 chemicals, the in vivo and in vitro data.
Chemical classification and structure may be associated with target disruption on the ATRA signaling pathway and skeletal defects based on our findings and those from the literature. Our analysis demonstrated that organochlorine pesticides (e.g., endosulfan, endrin, lindane, pyridaben, etc.) consistently activate DR5. Fungicides of the "-azoles" (e.g., triazoles, imidazoles, pyrazoles) family have broad association with activation of CYP surrogate biomarker, RARA, and DR5 (Figure 4). Organotin

FIGURE 5
Heatmap dendrogram of 48 chemicals associated with ATRA pathway disruption and skeletal defects in ToxRefDB and ToxCast/ Tox21 with potential AOP elucidation. Reflected are disruption activity of 8 targets of interest (RARA/B/G, RXRA/B, RSP (Tox21), DR5, and CYP surrogate biomarker). Scale of 0-5, with lightest color indicating less target disruption activity and darker red indicating higher target disruption activity compared to other chemicals for the respective targets. The dominance of activity at DR5, CYP surrogate biomarker, and RARA is consistent with k-means and hierarchical clustering diagrams; however, the heatmap attributes a greater amount of sensitivity of DR5 to chemicals of interest through the consistently darker red shading. CYP26 and RARA follow DR5 as the most disrupted of the ATRA signaling pathway targets.
Frontiers in Pharmacology frontiersin.org 13 biocides preferentially activate RXRs (e.g., triphenyltin hydroxidetributyltin, tributyltetradecylphosphonium chloride). Phenylureas (e.g., linuron, forchlorfenuron) activate CYP surrogate biomarker. Compounds without common chemical classifications (e.g., SAR 150640, asprin, iprione, etc.) regularly activated DR5, followed by CYP surrogate biomarker, then RARA, and occasionally activated the other molecular targets of interest. The similar ToxPi outcomes between chemicals of the same class and comparable structures (analogue) provide evidence that these chemicals may be able to predict bioactivity information and skeletal defect outcomes for other similar compounds (Marvel et al., 2018). This "read-across" method would provide valuable predictive toxicological and hazard index information about chemicals without previous experimentation and at a faster rate than lab work.
We developed a computational model for predictive toxicology of 48 chemicals associated with chemical disruption of the ATRA signaling pathway with concurrent fetal skeletal defects. We found literature relationships between 48 chemicals and skeletal defects in various animal models (S1) consistent with our computational models. For example, mouse and several other animal models showed reductive limb defects resembling those seen in vitamin A deficiency (VAD); truncation and deletion of the long bones (ulna, radius, and humerus in the forelimb) and digital deletions and fusion (Kochhar, 1973(Kochhar, , 1985. Retinoic acid administration to mothers produced severe, multiple skeletal defects and more specific malformations involving the axial skeleton, the fore-and hindlimbs, and cleft palate in embryos (Kistler, 1981). Embryonic exposure to pyrethroids causes craniofacial malformations, body axis curvature in zebrafish embryo model (DeMicco et al., 2010). Chlorothalonil (2,4,5,6tetrachloro-1,3-benzenedicarbonitrile) is a broad-spectrum, nonsystemic, organochlorine fungicide and exposure to this organochlorine fungicide causes skeletal malformations including absence of xiphoid process, malformation of the supraoccipital, absence of caudal vertebrae, incomplete ossification and malformation of xiphoid process in rats (Silva et al., 2020). Furthermore, five organochlorine pesticides, specifically, chlordane, endrin, and endosulfan, activate human RAR-mediated gene transcription via a retinoic acid response element (RARE such as DR5) and these organochlorine pesticides cause skeletal defects in several animal models (Ottolenghi et al., 1974;Lemaire et al., 2005;Kamata et al., 2008;Inoue et al., 2010). Moreover, organotin biocides preferentially activate RXRs and induced skeletal defects in animal models. For example, pregnant mice and rats exposed to triphenyltin hydroxidetributyltin and tributyltetradecylphosphonium chloride showed skeletal malformations including skeleton variations (e.g., poorly ossified skull bones) and malformations (misshapened axis and skull bones) an increased incidence of fusion of ribs, irregular ossification of sternebrae, and bipartite vertebral centers and cleft palate in mice and rat embryos (Sarpa et al., 2007).
Other classified chemicals from the 48 compounds such as phenylureas (e.g., linuron) activate CYP surrogate biomarker and induce skeletal abnormalities including abnormal axial rotation, AOs are phenotypes resulting from stage and positional alterations in the fetal skeleton (Knudsen et al., 2021). MIEs were largely derived from quantitative data, while KEs were developed through the literature, and AOs were gleaned from a combination of literature and datamining.
Frontiers in Pharmacology frontiersin.org and abnormal branchial arch development in rat embryo (McNutt and Harris, 1994), with linuron exposure to zebrafish models inducing spinal curvature and tail specific skeletal defects (Maharaj et al., 2020). Furthermore, CYP26A1(−/−) knockout mouse embryos having limited ATRA degradation had posterior elongation of axial patterning and displayed severe truncation of posterior structures (Rhinn and Dollé, 2012). Several experimental studies have shown that retinoic acid inhibits chondrogenesis and that this inhibition is most likely RARA/B/G-mediated (Eckhardt and Schmitt, 1994;Schroeder et al., 1994). Additionally, studies have shown that an RARAspecific agonist is significantly more potent in inducing skeletal malformations more frequently than RARB-or RARG-specific agonists (Elmazar et al., 1996). RARG is abundantly expressed during chondrocyte maturation, whereas RARAs appears to be expressed in most other tissues of the embryonic limb development. Therefore, under conditions of excessive retinoid where expression of RARB (Mendelsohn et al., 1991;Jiang et al., 1994), RARA, and DR5 may be elevated within the limb mesenchyme region, and can lead to a delay or inhibition of chondrocytes and chondroblast differentiation, resulting in abnormal skeletogenesis. Furthermore, pregnant mice treated with commercially available RARA/B/G-ligands and ATRA in chronological sequences exhibited external, visceral and skeletal malformations by E18 of gestation. In these studies, teratogenic RAR-ligand potency rankings of RARA > RARB > RARG-these are consistent with RAR-ligands' potency ranking found in our study. This study concluded that these retinoids also produced a different spectrum of developmental defects, specifically, the A-ligand and G-ligand induced the mandible, limb malformations, bone ossification deficiencies and defects of the sternebrae and vertebral body-with relative consistency among our compounds associated with these target disruptions and congruent skeletal defects (Elmazar et al., 1996).
The identification of RAREs (e.g., DR5) by RAR-binding assay in animal studies, in vitro cell line transfection assays, and ligand-binding transgenes reporter assay in mouse have been used to identify ATRA targeted genes that may be required for embryonic development (Cunningham and Duester, 2015). There were several studies that report the insufficiency or excess of ATRA due to specific RAREs (e.g., DR5) may lead to developmental defects in embryo (Houle et al., 2003;Kumar et al., 2016). DR5's ability to influence a wide variety of skeletal defects is represented in our study as well as throughout the literature. For example, DR5, functioning as a RARE enhancer, can overactivate the Homeobox protein Hox-a1 (Hoxa1) gene, and disrupt hindbrain development (Langston and Gudas, 1992;Dupé et al., 1997;Houle et al., 2003;Pouilhe et al., 2007). Whereas, RARE that function as a silencer to repress Fibroblast growth factor-8 (Fgf8) expression in caudal region causes axial defects in embryo (Kumar et al., 2016). Hypoactivation of Transglutaminase 2 (Tgm2) by DR5 leads to abnormal interdigital limb development (Nagy et al., 1996;Dupé et al., 1999).
Chemical disruption in the ATRA signaling pathway can be linked to either elevated or diminished retinoid concentrations in target tissues, resulting either in up-or downregulation of RARmediated signaling, respectively. Under these circumstances, the model presented in the current manuscript covers a condition of only elevated ATRA concentration and upregulated RARsignaling resulting from environmental exposure. In contrast, chemicals that disrupt retinol transport or ATRA bioactivation would be expected to invoke developmental phenotypes similar to Vitamin A deficiency (Knudsen et al., 2021). For example, a human biomonitoring study provided evidence associating reduction in circulating retinol levels with the adipose accumulation of persistent organic pollutants (POPs), including some polychlorinated biphenyl (PCB) congeners and dichlofol (Galbán-Velázquez et al., 2021). One mechanistic scenario suggested by those authors was chemical CYP oxidation and induction of an oxidative stress microenvironment that, in turn, lowers retinol levels by interfering with retinol binding protein-4 (RBP4) as a surrogate marker for serum retinol that forms a circulating carrier complex with transthyretin (TTR). The absence of HTS assays for retinol transport and metabolism is a limitation of the current model as it may apply to conditions associated with diminished retinoid concentrations and downregulated retinoid signaling.
Not all relevant biomarkers of retinol delivery were covered in this project. Retinoid delivery to the fetus is frequently mediated by RBP4, however, this is not the only pathway; postprandial retinoids delivered via chylomicrons is also an important pathway by which retinoids are delivered to fetuses (Steinhoff et al., 2021). This additional pathway illuminates the fact that there are multiple pathways of retinol delivery to the fetus, providing the occasion for others to study additional delivery methods through models such as found in this project.
Various CYPs with potential effects on ATRA metabolism could be altered in their expression by activation of AhR, CAR, and PXR pathways that, like RARs also heterodimerize with RXRs (McSorley and Daly, 2000;Esteban et al., 2021). The present study focused on biochemical (cell-free) CYP assays in ToxCast that, unfortunately, do not include CYP26 (this deficiency is currently being updated in the Tox21 portfolio).
Retinoids are inherent component of human and animal diets, with aquatic animals most exposed in early life stages (Kubickova et al., 2021). A broad range of retinoids are used to influence skin structure development cosmetically, as well as for use as dermatological, hematological, and infectious disease treatments (Zasada and Budzisz, 2019;Cosio et al., 2021). Environmental residues from chemicals such as the antifungal azoles may occur by applications on skin or environmentally through wastewater or discarded pharmaceuticals (Daughton and Ruhoy, 2009 Socio-economic factors have been found to influence maternal and fetal exposure to many chemical compounds (Renzo et al., 2015;Sprinkle and Payne-Sturges, 2021). A limitation of this study is high-throughput screening data compiled is at times derived from a limited set of not-diverse cell lines, resulting in a lack of understanding about the implications of our findings on populations that may be most affected by exposure to these chemical compounds. Greater diversity in original tissue and molecular samples would allow for greater direct applicability of study findings to more people, and further illuminate the extent of the impact of environmental factors on developmental toxicity in populations facing environmental injustice.
Heatmaps, k-means clusters, and hierarchical clusters provided clear quantitative associations between chemicals and skeletal defects and target disruption that inform AOP formation. Our pAOPs provide linear, organized, and biologically-relevant perspectives about MIEs determined by quantitative data that lead to literature derived KEs, and concludes with qualitative and quantitative literature and data derived ASOs (Table 3). However, AOPs are deterministic, resulting in less dynamic and vigorous conclusions. AOPs fail to reflect the biological complexity of the realm in which MIE, KE, and AO activities occur. As NAMs become more refined with increased data availability and defining fit-for-purpose criteria, their results will have greater dependability and applicability (Parish et al., 2020). Furthermore, ToxCast and Tox21 data's gene score transformations aims to exclude nonspecific activity due to cell stress and cytotoxicity, however, these factors may lead to skeletal developmental defects . At the same time, cell stress and cytotoxicity exclusion may have resulted in false positives due to global effects such as target bioactivity levels, which were still collected due to data collection methods. We propose the use of Agent Based Modeling (ABM), a dynamic in silico method that will test the effects of these chemicals on skeletal appendicular development. If the data, analytical tools, and AOPs developed here are consistent with ABM findings, they will provide substantial support for the conclusions reached. ABMs allow for high quality sensitivity analysis, further quantification of predictions, and generation of sound hypotheses formulation Knudsen et al., 2020). Computational models also allow for a larger-scale evaluation of high-throughput screening/high content screening data and the ability to derive results that cannot be obtained experimentally (Knudsen et al., 2021).

Author contributions
This work is represented by Tox21 Cross-Partner Project 13 as a collaboration between USEPA and the NIH (NICEATM, NCATS). JP led the data mining and visualization, BA supported data visualization and interpretation, NB led the literature mining, NK participated in the ToxCast-ToxRefDB in vitro-in vivo model inference, MX contributed to Tox21 data inference, and TK was the project lead. All authors contributed to the problem formulation, study design, modeling, referencing, and manuscript draft.

Funding
This project was supported by US EPA's Center For Computational Toxicology and Exposure (CCTE) as part of CSS 5.3. This study was also supported in part by the Interagency Agreement IAA #NTR 12003 from the National Institute of Environmental Health Sciences/Division of the National Toxicology Program to the NCATS, and the Intramural Research Program of the NCATS, National Institutes of Health.

Author disclaimer
The views expressed in this article are those of the authors and do not necessarily reflect the views or policies of the U.S. Environmental Protection Agency. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.

SUPPLEMENTARY FIGURE S2
ToxPi results for 117 chemicals with pAOPs associated with disruption of elements of RA signaling pathway. Genes of interest include RARA/ B/C, RXR A/B, DR5, Tox21 Pathway, CYP surrogate biomarker. Most potent chemicals associated with the greatest amount of gene disruption have the largest ToxPis and highest ToxPi rank (48). K-means Clustering ATRA Signaling Pathway for 117 chemicals.
Disrupted activity in the following targets is of interest: RARA/B/C, RXR A/B, DR5, RSP (Tox21), biomarker. The straight-lined clusters (blue and green clusters) are likely due to activity in RSP (Tox21), but lack of experimentation and complimentary assays with ToxCast. PC2 has greatest potency decreasing from top to bottom. In vitro ToxPis have RSP (Tox21) as most sensitive (larger slices) on the top, with decreasing sensitivity (smaller slices) on the bottom of the PC2. In the PCA plot there is an obvious effect on ToxPi composition, though there is more data needed for a fuller picture. The straight-lined clusters are likely due to activity in Tox21, but lack of experimentation with ToxCast. PCA2 is greatest potency from top to bottom. In vitro ToxPis have RSP (Tox21) as most sensitive (larger slices) on the top and to and decreasing sensitivity (smaller slices) on the bottom of the PCA2. In the PCA plot there is an obvious effect. However, the hierarchical clustering plot is not informative with these 117 chemicals. Heatmap dendrogram of 117 chemicals associated with ATRA pathway disruption and skeletal defects in ToxRefDB and ToxCast/Tox21 with potential AOP elucidation. Reflected are disruption activity of 8 targets of interest (RARA/B/G, RXRA/B, RSP (Tox21), DR5, and CYP surrogate biomarker) in association with 117 chemicals of interest. Scale of 0-5, with lightest color indicating less target disruption activity and darker red indicating higher target disruption activity compared to other chemicals for the respective targets. The heatmaps attributes a greater amount of sensitivity of DR5 to chemicals of interest through the consistently darker red shading. In order, DR5, CYP26, and then RARA are the most disrupted of the ATRA signaling pathway targets in the heatmap.