Microbiota alter metabolism and mediate neurodevelopmental toxicity of 17β-estradiol

Estrogenic chemicals are widespread environmental contaminants associated with diverse health and ecological effects. During early vertebrate development, estrogen receptor signaling is critical for many different physiologic responses, including nervous system function. Recently, host-associated microbiota have been shown to influence neurodevelopment. Here, we hypothesized that microbiota may biotransform exogenous 17-βestradiol (E2) and modify E2 effects on swimming behavior. Colonized zebrafish were continuously exposed to non-teratogenic E2 concentrations from 1 to 10 days post-fertilization (dpf). Changes in microbial composition and predicted metagenomic function were evaluated. Locomotor activity was assessed in colonized and axenic (microbe-free) zebrafish exposed to E2 using a standard light/dark behavioral assay. Zebrafish tissue was collected for chemistry analyses. While E2 exposure did not alter microbial composition or putative function, colonized E2-exposed larvae showed reduced locomotor activity in the light, in contrast to axenic E2-exposed larvae, which exhibited normal behavior. Measured E2 concentrations were significantly higher in axenic relative to colonized zebrafish. Integrated peak area for putative sulfonated and glucuronidated E2 metabolites showed a similar trend. These data demonstrate that E2 locomotor effects in the light phase are dependent on the presence of microbiota and suggest that microbiota influence chemical E2 toxicokinetics. More broadly, this work supports the concept that microbial colonization status may influence chemical toxicity.

Essentially all known multicellular organisms are naturally colonized by a diverse array of microbial communities comprised of bacteria, viruses, archaea, fungi and protozoa. Current estimates suggest that the adult human gastrointestinal tract alone harbors more than 100 trillion bacteria from over 1,000 different species 1 . The metagenome of these gut microbiota encodes 100-150 times more genes than the human genome, with rapid plasticity and unique biological functions 2 . Microbiota can interact with their hosts directly or through the production of various microbial products, such as extracellular enzymes and cell wall components 3 to influence immune and nervous systems, metabolism, and behavior [4][5][6] . Microbiota can also interact with drugs and environmental chemicals in ways that can shift their health effects or toxicity profiles [7][8][9] . Such findings indicate that host microbiota can be an important determinant of both homeostasis and disease risk.
Neurodevelopmental disorders like autism-spectrum disorder and attention-deficit/hyperactivity disorder have increased in prevalence over the last decade 10 . Emerging evidence indicates that microbiota may mediate or influence neurologic development. During early life, different stages of brain development coincide with the establishment and development of the gut microbiome 11 . Recent studies have associated some neurodevelopmental disorders with altered microbiota profiles 12 , and experiments that use axenic (microbe-free) animals or antibiotics to deplete host-associated microbiota indicate that microbiota influence key aspects of early neurodevelopment. Affected processes include formation of the blood-brain barrier 13 , myelination 14 , neurogenesis 15 , neurotransmitter levels 16 , and behavior [17][18][19] . Collectively, these interactions are driven by the microbiota-gut-brain axis, which involves bidirectional communication between intestinal microbiota and the central nervous system. E2 exposure does not impact zebrafish-associated microbial communities. To identify whether altered behavior following E2 exposure co-occurred with changes in microbial community structure or predicted microbial community function, conventionally colonized zebrafish were exposed to five non-teratogenic E2 concentrations, and 16S rRNA gene sequencing was used to identify operational taxonomic units (OTUs) present in each sample. Non-metric multidimensional scaling (NMDS) analysis revealed that global compositional profiles of microbiota did not cluster by E2 treatment (ANOSIM OR = 0.064) (Fig. 3A). Bray-Curtis similarity indices for all between-treatment group comparisons (e.g. DMSO vehicle control vs. 0.34 µM E2) were similar and ranged on average from 55.8-68.1% (Fig. S1). Taxonomic analysis revealed that Proteobacteria was the most abundant phylum at all concentrations tested (data not shown). In addition, 24 bacterial families captured >95.6% of all OTUs identified represented, and the relative abundances of family-level taxa were similar across all DMSO control and E2-exposed larvae (Fig. 3B). E2 exposure also did not impact any alpha diversity metrics including total number of species, Margalef 's species richness index, species evenness, Shannon's diversity index, or Simpson's diversity (Fig. S2). Relative abundances of PICRUSt-generated Level 2 KEGG predicted functions were similar across all treatment groups (Fig. 3C). To assess whether E2-exposure impacted function at a more granular level, linear discriminant effect size (LEfSe) analysis was subsequently used to identify differentially abundant Level 3 KEGG functions (Table S1). All 11 identified functions were not significantly altered in E2-exposed larvae relative to DMSO vehicle controls (Fig. S3).
E2-induced locomotor effects in the light phase are dependent on microbiota. To examine whether hypoactivity observed in conventionally colonized larvae following E2 exposure (Fig. 2B,C) was related to microbial colonization of the host, three cohorts of zebrafish (conventionally colonized, axenic colonized on day 1, and axenic) 19 (Fig. 1) were continuously exposed to 0.4 or 1.2 µM E2. As expected, DMSO control axenic zebrafish were hyperactive in the dark phase relative to conventionally colonized or axenic colonized on day 1 zebrafish (p < 0.0001, Figs 4A-I and Fig. S4). At 10 dpf, conventionally colonized larvae exhibited significant hypoactivity in the light phase at 0.4 and 1.2 µM E2 (p < 0.00001 and p = 0.0407, respectively, Fig. 4B), similar to www.nature.com/scientificreports www.nature.com/scientificreports/ the initial behavior experiment in conventionally colonized larvae (Fig. 2BC). In the dark phase, significant hypoactivity was observed at 0.4 µM E2, but not 1.2 µM E2 (p = 0.0043 and p = 0.5862, Fig. 4C). In axenic colonized on day 1 larvae, significant hypoactivity was observed in both the light and dark phases at 0.4 and 1.2 µM E2 (light: Figure 3. E2 exposure does not impact community structure or predicted function of zebrafish-associated microbiota. 16S rRNA gene sequencing was performed on 10 dpf larvae developmentally exposed to E2, and PICRUSt was used to predict metagenomic function. Non-metric dimensional scaling (NMDS) ordination plots are shown for microbial community structure. The Bray-Curtis similarity index was used to assign percent similarity values. Analysis of similarities (ANOSIM) ordered R (OR) statistic is shown. Percent relative abundances of (B) family-level taxa or (C) Level 2 KEGG predicted functions in 10 dpf larvae are shown. n = 4 biological replicates with 10 larvae per replicate. p = 0.0033 and p < 0.0001, respectively, dark: p = 0.0259 and p = 0.0004, respectively) ( Fig. 4E,F). In contrast to colonized larvae, axenic larvae exhibited no significant change in locomotor activity during the light phase compared to DMSO vehicle controls (Fig. 4H). In the dark, axenic larvae behaved like colonized larvae exhibiting locomotor hypoactivity at 0.4 and 1.2 µM E2 (p = 0.0041 and p = 0.0026, respectively) ( Fig. 4I).
Axenic zebrafish contain more parent E2 than colonized zebrafish. To compare measured E2 concentrations in exposure media and whole-body tissue homogenates (that included intestinal tracts) obtained from conventionally colonized, axenic colonized on day 1, and axenic larvae continuously exposed to 0.4 or 1.2 µM E2, samples were collected and analyzed using targeted chemistry techniques. Immediately following dosing on day 1, measured E2 media concentrations were similar across colonization statuses and ranged from 0.242-0.453 µM and 0.910-1.485 µM for the 0.4 and 1.2 µM E2 exposure groups, respectively (Dose, p < 0.0001, Fig. 5A and Table S2). In day 10 media samples, significant main effects of dose (p < 0.0001) and colonization status (p < 0.0001) were observed (Fig. 5B). Average measured media concentrations ranged from 0.155-0.509 µM for the 0.4 µM E2 group, and 0.975-1.76 µM for the 1.2 µM E2 group ( Fig. 5B and Table S2). Individual doses were not compared across colonization status because no significant interaction between dose and colonization status was observed. Analysis of day 10 zebrafish tissue revealed distinct colonization-dependent differences in measured E2 concentration in zebrafish exposed to 1.2 µM E2. In addition to significant main effects of dose (p < 0.0001) and colonization status (p < 0.0001), a significant interaction between these two terms was observed (p < 0.0001), indicating that the pattern of measured internal E2 doses across all concentrations tested differed by colonization status (Fig. 5C). Subsequent pairwise comparisons revealed that following exposure to a nominal concentration of 1.2 µM E2, measured E2 tissue doses in axenic larvae were 2.5-3.7 times higher than in conventionally colonized or axenic colonized on day 1 larvae (Fig. 5C and Table S2). At a nominal concentration of 0.4 µM E2, tissue doses were comparable across colonization status ( Fig. 5C and Table S2).
Axenic zebrafish contain higher abundances of some E2 metabolites. As a hypothesis-generation strategy, a non-targeted approach was initially selected to identify any chemical features (e.g. E2 metabolites or microbial metabolites) that may be altered in conventionally colonized, axenic colonized on day 1, or axenic zebrafish following exposure to 0.4 or 1.2 µM E2. Following initial feature filtering, 94 features were retained www.nature.com/scientificreports www.nature.com/scientificreports/ in the dataset (Fig. S5). Comparison of features that met fold-change threshold criteria (≥2-fold) within colonization status (e.g. axenic 0.1% DMSO vs. axenic 1.2 µM E2) yielded numerous significant features in all three cohorts ( Fig. S6A-C). Between-status comparisons at 1.2 µM E2 (e.g. conventionally colonized 1.2 µM to axenic 1.2 µM E2) revealed that 33 features were concordant (i.e. increased or decreased) in both conventionally colonized and axenic colonized on day 1 larvae compared to axenic larvae. Of these 33 features, 26 were significantly decreased and seven were significantly increased ( Fig. S7A-C and Table S3). However, using this stringent filtering approach, only one feature was identified as a predicted E2 metabolite (E2 sulfate) (Fig. 6B), and other differentially represented unknown metabolites were not identified. Therefore, to specifically identify additional E2 metabolites that were differentially represented following exposure to the parent compound, the raw LC-MS data were screened for known metabolites in the E2 metabolism pathway (Fig. 6A). Five chemicals consistent with known metabolites were tentatively identified: E2 sulfate, E2 glucuronide (A and B, two isomers), estriol (E3) sulfate, estrone (E1) sulfate, and E1 glucuronide (Fig. 6A, Table 1). Integrated peak areas for all metabolites were significantly correlated (p < 0.0001) with measured E2 tissue concentrations. The strongest correlations were observed for E2 sulfate (r = 0.84) and the two isomers of E2 glucuronide (r = 0.87 and 0.90) ( Table 1). Boxplots for these E2 metabolites ( Fig. 6B-D) show very similar trends to that observed for parent E2 concentrations (Fig. 5C). Specifically, peak areas of these metabolites were noticeably elevated in the axenic 1.2 μM E2 group. Correlations were less pronounced for E3 sulfate (r = 0.73), E1 sulfate (r = 0.69), and E1 glucuronide (r = 0.69) (Table 1)  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Recent work suggests that estrogens and host-associated microbiota, particularly those within the intestinal tract, to influence a variety of health outcomes including brain function, behavior, obesity, and cancer 28 . In this study, we used colonized and axenic zebrafish as an experimental system to characterize the effects of E2, a potent estrogen receptor agonist, on microbiota composition and predicted metagenomic function. In addition, we asked whether E2 exposure impacted locomotor activity in a colonization status-dependent manner. We report that E2 did not impact microbial community structure or predicted function. However, E2 exposure did result in significant hypoactivity in the light phase in colonized larvae only. Measured parent E2 concentrations were ~3 times higher in axenic relative to colonized zebrafish. Predicted E2 metabolites E2 sulfate and E2 glucuronide were also higher in axenic relative to colonized zebrafish. In contrast, E3 sulfate, E1 sulfate and E1 glucuronide were similar in colonized and axenic larvae. To our knowledge, this is the first study to demonstrate that microbiota interact with an exogenous chemical to influence light phase locomotor effects. Taken together, these data support the concept that microbiota and environmental chemicals can interact to influence host physiology.
In mammals, the estrogen signaling pathway has been implicated as an important regulator of various cognitive and behavioral responses. For example, low estrogen levels have been associated with a decline in memory-related functions in humans and animal models 41,42 . Similar to the hypoactivity data presented here, exposure to E2 in ovariectomized mice resulted in locomotor hypoactivity in a light/dark transition test 43 . The Figure 6. Axenic zebrafish contain higher concentrations of some putative E2 metabolites. LC-MS was performed on whole zebrafish tissue. (A) E2 metabolism reference schematic is shown. Black boxes indicate putative E2 metabolites identified in this study. Boxplots with integrated peak area for (B) estradiol sulfate, (C) estradiol glucuronide A, (D) estradiol glucuronide B, (E) estriol sulfate, (F) estrone sulfate and (G) estrone glucuronide in conventionally colonized, axenic colonized on day 1, and axenic larvae are shown. A concentration of 0 µM refers to 0.1% DMSO vehicle controls. Spearman rank correlation coefficients were calculated to assess the relationships between integrated peak areas for each detected E2 metabolite and measured E2 tissue concentrations from targeted analysis (Table 1). n = 2-4 (10 larvae per biological replicate).
www.nature.com/scientificreports www.nature.com/scientificreports/ same study also reported increased time spent in the outer zones in an open field test 43 . These two widely used mammalian tests measure anxiety-like behavior, suggesting that, in animals lacking endogenous estrogens, E2 exposure increased anxiety-like behavior and heightened fear response, both of which are acquired adaptive survival behaviors 44 . In a female rat study, E2 treatment during the neonatal period resulted in improved learning and enhanced memory performance in adulthood 45 . These effects were hypothesized to be related to E2 effects on γ-aminobutyric acid (GABA) receptor expression and function in the adult brain 45 . Together, these studies suggest that exogenous E2 exposure has widespread effects on nervous system function in mammalian models.
The circuits that underlie the patterned development of stereotypical locomotor behaviors are highly regulated, and it is now understood that microbiota play a critical role in the development of many normal locomotor phenotypes and regulate neurophysiological behaviors through alteration of neural, endocrine, and immune pathways 46 . We have previously shown that axenic animals lacking microbiota exhibit a consistent increase in dark phase locomotor activity 19 . Using a standard light/dark assay to investigate whether exposure to 0.4 or 1.2 µM E2 caused a colonization status-dependent behavioral response, we showed that only colonized larvae exhibited light phase locomotor hypoactivity. While both colonized and axenic larvae exhibited locomotor hypoactivity in the dark phase, these results are consistent with other studies that demonstrate phase-specific behavioral effects. For example, in other larval zebrafish studies that used similar light/dark assays, developmental exposure to the flame retardants 2,2′,4,4′-tetrabromodiphenyl ether (BDE-47), tricresyl phosphate (TMPP), and phenol, isopropylated, phosphate (3:1) (IPP) caused dark phase-specific hyper-or hypoactivity 47,48 . Developmental 2,2′,4,4′,5-pentabromodiphenyl ether (BDE-99) exposure also resulted in differential phase-specific effects, with hyperactivity observed in the light phase and hypoactivity observed in the dark phase 47 . In addition, Kokel et al. showed that embryonic zebrafish exhibit shared and specific behavioral "barcodes" based on their locomotor phenotype in response to a photic stimulus over four different phases following exposure to a large-scale set of drugs 49 . This approach allowed for the classification of putative biological mechanisms of action based on barcode similarity to pharmacological agents with known modes of action. Collectively, these data show that zebrafish light/dark swimming behavior is broadly affected by xenobiotic exposures and that environmental chemicals and drugs produce a diverse set of both phase-shared or phase-specific locomotor effects. This evidence further suggests that changes in light phase-specific locomotor behavior arise from multiple biological pathway-level perturbations. While more work is needed to understand the biology that controls zebrafish swimming behavior in response to light or dark stimuli, the data presented here are in line with previous studies that show differential phase effects of xenobiotic exposures on locomotor activity [47][48][49] . In addition, our results are the first to show that changes in the stereotypical locomotor response to a light stimulus following exposure to an exogenous environmental chemical are influenced by microbiota.
Microbiota can interact with environmental chemicals via toxicodynamic (alterations in microbial composition or function) or toxicokinetic (microbial xenobiotic biotransformations) mechanisms. In this study, due to logistical challenges associated with our unique concentration-response design, we did not specifically analyze the larval gut microbiome. However, metagenomic analysis revealed that composition of the DMSO vehicle control zebrafish microbiome in this study was similar to that reported from larval gut microbiomes of 10 dpf fish in Stephens et al. 37 . Specifically, at the phylum and class levels, microbiota obtained from both whole zebrafish larvae or larval intestinal tract dissections at 10 dpf was dominated by Proteobacteria and Gammaproteobacteria, and contains families that align with core genera identified in adult laboratory and wild zebrafish in Roeselers et al. (i.e. Aeromonadaceae, Shewanellaceae, Enterobacteriaceae, and Comamonadaceae) 50 . These similarities are likely because the intestinal tract contains the most bacteria relative to other organs, like in humans 36 . Some differences in composition between microbiota profiles obtained from whole larvae or dissected intestinal tracts were also noted. Namely, the Firmicutes phyla comprised on average 4.2% of all bacterial taxa at 10 dpf in Stephens et al. In our study, although we did detect Firmicutes (20/232 OTUs), the relative abundance of this phylum was <1% of all sample totals. At the family level, higher relative abundances of Chromatiaceae were also observed in our study compared with Stephens et al. However, higher levels of Chromatiaceae were not observed in all DMSO controls across consecutive experiments done previously in our lab 51 , suggesting that microbiota variability exists within and across aquaculture systems and can be due to a variety of factors including rearing temperature, zebrafish strain, and diet 35 www.nature.com/scientificreports www.nature.com/scientificreports/ structure and putative function of zebrafish-associated microbiota were not altered. Microbiota are accustomed to interacting with endogenous E2 and play a significant role in the uptake and metabolism of estrogens 52 . Therefore, it is plausible that exposure to a key endogenous hormone like E2 does not modify the composition of host-associated microbiota in larval zebrafish in the current study. Although one previous study in adult zebrafish using a static-renewal system demonstrated disruption of microbiota following exposure to a nominal concentration of 1 µg/l (3.6 nM) E2 53 , these contrasting results are not surprising given that microbial composition changes with life stage 37 .Taken together, the results presented here suggest that toxicokinetic mechanisms (e.g. microbial-mediated E2 metabolism) may be driving the E2-induced light phase alterations in locomotor activity.
In the current study, measured internal E2 dose varied depending on colonization status. Paradoxically, 10 dpf axenic zebrafish contained significantly more (2.5-3.7 times) parent (unconjugated) E2 compared to conventionally colonized or axenic colonized on day 1 larvae. Non-targeted analysis revealed that peak areas of putative metabolites E2 sulfate and E2 glucuronide were also elevated in axenic larvae, while peak areas of E3 sulfate, E1 sulfate, and E1 glucuronide were similar in colonized and axenic zebrafish. In general, estrogen metabolism occurs primarily in the hepatobiliary system and the intestinal "estrobolome" (the bacterial genes that are capable of metabolizing estrogens) 30 , resulting in potential subsequent detoxification or activation. One explanation for the observed locomotor phenotype following chemical exposure is that E2 is deconjugated and then converted to a more active metabolite, either by microbial enzymes or a more metabolically competent host. Previously published work shows that microbiota bioactivate pharmacological agents like diclofenac, indomethacin, ketoprofen, and sulfasalazine 54,55 . As estrogenic activity is often cell, tissue, model, or endpoint-specific, a more bioactive estrogen profile could be driven by increased deconjugation and bioavailability 56 or metabolites with enhanced local ER binding or activation 57,58 . For example, sulfonated estrogens, particularly E2 sulfate, can circulate in the blood and reach estrogen target tissues where they may be deconjugated and converted to receptor-active free estrogens 59 . In addition, a study examining the competitive binding activity of dehydroepiandrosterone (DHEA), an abundant steroid hormone that serves as a precursor to estrogen and androgen, and it's sulfonated metabolite (DHEA-S), showed that DHEA-S, but not DHEA, was competitive with E2 for ERɑ and ERβ binding and stimulated MCF-7 cell proliferation 60 . Parent E2 and E2 glucuronide have also been shown to bind to and activate ERα in the liver 61 . Furthermore, under different physiological conditions, E2 metabolites can differentially activate ERɑ and ERβ signaling pathways, elicit downstream gene activation, and trigger intracellular signaling cascades 62 , which may influence neurodevelopment.
Important limitations of the non-targeted approach used in this study included the inability to detect the parent compound, unknown differences in extraction efficiency, column affinity, breakdown products, ionization efficiency 63 , and a lack of analytical standards for E2 metabolites. While these limitations precluded direct comparisons between levels of E2 metabolites within an exposure group and across colonization statuses, these data support the targeted results and suggest that colonization status influences E2 metabolism. In future follow-up studies, quantification using a targeted approach would allow for more detailed comparisons of E2 metabolite profiles in colonized and axenic zebrafish. It is also possible that E2 metabolites potentially contributing to the observed behavioral effect in colonized larvae in the light phase may not have been identified. For instance, we did not detect hydroxylated E2 metabolites that have been reported in mammals 64 . While more work is needed to enhance detection and determination of E2 metabolites in zebrafish, the application of non-targeted analysis led to the simultaneous identification of five unique sulfonated or glucuronidated E2 metabolites with unknown toxicity profiles that may contribute to colonization-specific alterations in locomotor activity. In addition, chemical toxicokinetics is not simply defined by biotransformation (e.g. bioactivation). As this study did not fully assess absorption, distribution, or excretion of E2 in colonized or axenic zebrafish, other potential mechanisms may also be related to the novel behavioral phenotype observed in the light phase in colonized larvae described here.
Host-associated microbiota are recognized as important modulators of the complex interactions between factors such as diet, sex, race, and life-stage and many human diseases 65 . An increasing body of evidence suggests that microbiota also mediate health effects of xenobiotic exposure including gastrointestinal toxicity and immune responses [66][67][68][69] . Our results indicate that microbiota are required for E2-induced light phase locomotor hypoactivity, further supporting the idea that microbiota and chemicals may interact to influence biological processes. This work also raises an interesting point related to hazard identification, suggesting that microbiota may be an important factor for characterizing chemical interactions with a host organism. Continued investigation of how estrogens and other environmental chemicals interact with the microbiota-gut-brain axis may uncover the biological mechanisms by which microbiota communicate with the nervous system and influence brain development and behavior. www.nature.com/scientificreports www.nature.com/scientificreports/ stock solution aliquots and performing serial dilutions in DMSO. Working solutions were stored in 4 ml amber glass vials at room temperature. All exposure groups contained a final concentration of 0.1% DMSO (v/v). Vehicle controls received 0.1% DMSO only. For targeted chemistry, E2-3,4-13 C 2 (internal standard, Product#: CLM-803-S; Lot#: SDDI-011) was obtained from Cambridge Isotope Laboratories, Inc. (Tewksbury, MA). Stock solutions were prepared in methanol and acetonitrile and stored at −20 °C. Intermediate standards were prepared fresh daily from stocks. All reagents and solvents were reagent grade or high-performance liquid chromatography (HPLC) grade.

Methods
exposures. E2 exposures were performed in T25 tissue culture flasks incubated at 26 °C on a 14 h:10 h light:dark cycle. Zebrafish were continuously exposed in a semi-static system from 1 to 10 days post-fertilization (dpf). For all experiments, chemical exposures began on day 1 (Fig. 1). All flasks were housed statically through 6 dpf, followed by daily renewal of chemical exposure solutions in concert with an 80% media change (0.2 µm filter-sterilized 10% Hanks Balanced Salt Solution (FS-10% HBSS) from 6-9 dpf. On days 6-9, each flask also received 75 kGy gamma-irradiated Gemma Micro 75 as a food source at a final concentration of 0.04% (v/v) (Fig. 1). Dead embryos or larvae were removed from each flask during media changes. For the developmental toxicity experiment, conventionally colonized (CC, colonized with aquaculture facility microbiota on day 0) zebrafish were exposed to 0.05-30 µM E2. These concentrations were selected based on previous zebrafish studies and zebrafish assays within the U.S. EPA ToxCast Dashboard (https://actor.epa.gov/dashboard/) 25,70,71 . For 16S rRNA gene sequencing, conventionally colonized zebrafish were exposed to non-teratogenic concentrations ranging from 0.34-3.5 µM E2. For 3-cohort behavior and mass spectrometry, conventionally colonized, axenic colonized on day 1, and axenic zebrafish were exposed to 0.4 or 1.2 µM E2.
For experiments involving all three zebrafish cohorts, media sterility was tested as previously described 19 . Briefly, at 1 and 10 dpf, two tryptic soy agar (TSA) plates (Sigma, #22091) were inoculated with 10 µl of media from each flask. At 10 dpf, the sterility of media from axenic flasks was further tested by inoculating 100 µL of flask media into tubes of Nutrient Broth (Sigma, #70122), Brain Heart Infusion Broth (Sigma, #53286), or Sabouraud Dextrose Broth (Sigma, #S3306). Plates and tubes were incubated at 26 °C under aerobic and anaerobic conditions for at least seven days. Contaminated flasks were excluded from the study.
Behavior testing. At 10 dpf, morphologically normal larvae were removed from flasks using a sterile transfer pipet and placed into 48-well plates containing 500 µl of FS-10% HBSS per well. Plates were sealed with Microseal A film (BioRad, #MSA5001), wrapped in Parafilm, and placed in the dark in a temperature controlled behavior testing room at 26 °C, for at least 2 hr prior to testing. For testing, microtiter plates were placed on a Noldus tracking apparatus and locomotor activity was recorded for 40 mins. The light program consisted of a 20-min acclimation period in the dark (0 lux) followed by a 10 min light period (5 lux) and a 10 min dark period (0 lux). All tests were carried out between 11:15 am-4:00 pm. Videos were analyzed using Ethovision software Version 12 (Noldus Information Technology) as described previously 48 .
DNA sequencing of 16S rRNA gene. Whole-body zebrafish homogenates were used to evaluate changes in microbial community structure and predicted function. DNA extraction and sequencing of the 16S rRNA gene was completed as previously described 19 . Briefly, 10 dpf larvae (n = 4 biological replicates, 10 larvae per replicate) were collected, anesthetized, and homogenized. DNA was isolated from each sample using a ZR-Duet ™ Analysis of 16S rRNA gene sequences. For microbial community structure and predicted metagenomic function analyses, paired-end sequences were trimmed at a length of 250 base pairs and quality filtered at <0.5% expected error using USEARCH v7 73 . Reads were analyzed using the QIIME 1.9.0 software package 74,75 . The www.nature.com/scientificreports www.nature.com/scientificreports/ total number of reads for each sample can be found in Table S4. A cutoff of 500 reads per sample was used for downstream analyses. One sample did not meet this criterion (3.5 µM E2_R4). A closed-reference OTU table was generated within QIIME-1.9.0 and used to assess microbial composition and generate PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) functional predictions 76 with the Greengenes 13_8 reference database (https://catalog.data.gov). Alpha and beta diversity analyses were performed using PRIMER 7 software (Primer-E v7.0.11) for total number of species, Margalef 's richness index, Pielou's evenness index, Simpson's index, Shannon's diversity index, Bray-Curtis similarities, and non-metric multidimensional scaling (NMDS) ordination plots as previously described 19 . Stacked bar plots were generated using OTUs or predicted functions that contributed at least 1% to any of the sample totals. The KEGG orthology classification scheme was used for functional annotations 77 and relative abundances of predicted functional pathways were formatted as previously described 78 . Kruskal-Wallis and pairwise Wilcoxon tests were conducted (ɑ = 0.05) using E2 as the main categorical variable. Linear Discriminant Analysis (LDA) scores >2.0 were used to assess functional enrichment. Microbial community structure was also analyzed with mothur-generated OTU tables and the Silva reference database as previously described 19,79 and provided for future reference as a taxonomic record (Table S5). sample generation and preparation for mass spectrometry. Conventionally colonized, axenic colonized on day 1 and axenic zebrafish larvae were exposed to 0.1% DMSO or 0.4 or 1.2 µM E2 as described above. Blank (fish-free) flasks containing gamma-irradiated Gemma Micro 75 only were included as a control. On day 1, 1 ml of media was collected from all flasks immediately after dosing to measure initial E2 concentrations. At 10 dpf, larvae were transferred from each flask to a vial containing a final volume of 500 µl FS-10% HBSS (n = 2-5 biological replicates, 10 larvae per replicate). Larvae were anesthetized on ice, flash frozen in liquid nitrogen, and stored at −80 °C. Media samples (1 ml) from each flask were also collected at 10 dpf and stored at −80 °C. Larval tissue samples were homogenized in 200 µl of 0.1% formic acid in acetonitrile containing the E2-3,4-13 C 2 internal standard using a FastPrep 24 homogenizer (MP Biomedicals, Santa Ana, CA) and ~150 mg, 1.0 mm diameter zirconia/silica beads. The homogenate was centrifuged at 14,000 rpm for 15 minutes at 4 °C. For analysis, 150 µl of the supernatant was diluted in a liquid chromatography (LC) vial with 150 µl of HPLC grade water. Exposure media was prepared for analysis by diluting 200 µl with 50 µl methanol containing the E2-3,4-13 C 2 internal standard. Calibration and verification standards were prepared in solvent: 0.1% formic acid in 50% acetonitrile and water for zebrafish tissue samples and 20% methanol in water for exposure media samples. Blanks and quality control samples were prepared by adding appropriate amounts of tissue homogenate or exposure media. targeted mass spectrometry and analysis. Targeted chemistry analysis was conducted by liquid chromatography-tandem mass spectrometry (LC-MS/MS) on an Accela UPLC and TSQ Quantum Ultra triple quadrupole mass spectrometer (Thermo Fisher Scientific, Waltham, MA) equipped with atmospheric pressure photo ionization (APPI) source with krypton lamp operated in positive ionization mode. LC separation was achieved using Kinetex EVO, 50 × 2.1 mm, C18, 2.6 µm, 100 å LC column (Phenomenex, Torrance, CA) with gradient elution at a flow rate of 360 µl/min using A (5% methanol and water) and B (5% water and methanol). Initial LC conditions, 80% A and 20% B, were held for 15 seconds followed by a linear ramp to 100% B at 8 minutes, and held at 100% B for 2 minutes. The column was equilibrated for 2 minutes at initial conditions. Total run time was 12 minutes. Toluene was added to the LC flow post column at 40 µl/min using a PHD Ultra syringe pump (Harvard Apparatus, Holliston, Massachusetts). Source conditions were optimized for the precursor [M-H2O + H] + ions of E2. Collision energy and collision pressure were optimized for product ions of E2 and E2-(3,4-13 C 2 ).
Integration, calibration, and quantitation were performed using Xcalibur 3.0.63 (Thermo Fisher Scientific). The observed measured range in zebrafish tissue was 486 fmole/larva -15 picomole/larva with a 56 fmole/larva method detection limit (MDL). Recovery of 80% was observed for samples spiked at 1.9 pmol/larva and matrix effects were observed to be less than 10%. The observed range in exposure media was 45.9 nM -2.31 µM with a 4 nM MDL. 86% recovery was observed for samples spiked at 459 nM and matrix effects were observed to be less than 10%. Batch results were accepted based on the following criteria: 1) The six standard calibration curve had a correlation coefficient >0.99 and accuracy tolerance ≤20%, 2) >75% of all individual quality control standards had sample accuracy tolerance ≤30% and %RSD (precision) ≤20%, and 3) blank response was < MDL.

Non-targeted mass spectrometry and analysis.
To identify endogenous and xenobiotic metabolites associated with E2 exposure and microbial status, non-targeted mass spectrometric analysis was performed within a four-week period on the same samples used for targeted mass spectrometry. Sample extract (10 µl) was injected in triplicate. Analyte separation was accomplished using Waters Acquity UPLC ® BEH C 18 (2.1 × 50 mm, 1.7 µm) connected to an Agilent 1290 Infinity II Liquid Chromatography system (Palo Alto, CA) equipped with a degasser, binary pump, and autosampler. The mobile phase flow rate for gradient elution was set at 200 µl/min using A (aqueous 0.1% formic acid) and B (acetonitrile with 0.1% formic acid). The gradient started with 90% A and 10% B for 2 min, followed by a linear ramp to 100% B after 15 min. This condition was kept for 5 min before returning to starting mobile phase conditions at 21 min. The column was re-equilibrated for the next injection for 9 min. www.nature.com/scientificreports www.nature.com/scientificreports/ Raw data were initially processed using Agilent Profinder vendor software (v.8.00) for molecular feature extraction and integration, with additional filtering using in-house scripts (Supplemental Methods). Chemical features (unique empirical mass-to-charge ratio (m/z) and retention time (min)) were tentatively assigned based on hits against the EPA Chemistry Dashboard (https://comptox.epa.gov/dashboard/) and an Agilent Personal Compound Database and Library (PCDL) of the METLIN database 80 . Monoisotopic ion masses for components of the E2 metabolism pathway (see Fig. 6A) were then separately extracted for each data file using a 10 ppm mass window. The extracted ion chromatograms were smoothed with a 9-point moving average and integrated with Agilent's Agile 2 algorithm. For samples which lacked an obvious peak, a time window corresponding to the average peak window was manually integrated to capture a noise value. If a sample showed no chromatographic noise at the elution time, the value was recorded as zero. statistical analyses. The AC 50 value for E2 (i.e. concentration that elicits a 50% inhibitory response) was obtained from the zebrafish developmental toxicity assay using the log(agonist) vs. normalized response-Variable slope equation (GraphPad Prism 7) 71 . A one-way analysis of variance (ANOVA) with Tukey pairwise comparisons was used for analysis of the zebrafish developmental toxicity assay. An analysis of similarities (ANOSIM) test (concentration = ordered factor) was used to assess changes in microbial composition. For comparisons of Bray-Curtis similarity scores, a permutational one-way analysis of variance PERMANOVA with pairwise Monte Carlo comparisons was used (p < 0.05). Alpha diversity metrics were analyzed using a Kruskal-Wallis and Dunn's pairwise multiple comparisons test (p < 0.05). All behavior data were analyzed using SAS v9.4 software using a mixed effects repeated measures model. Each individual fish was used as a subject for the repeated measures, with the movement measures across phase (10 min light or 10 min dark period) and time (five time points within each phase, reflecting each two min period). Plate, flask and experiment effects were tested as random factors, found to be not significant (ɑ = 0.05), and subsequently removed from the model. Main effects of each fixed factor (i.e. concentration, phase, or time), and any interaction between or among the factors, were tested within each status. Backwards stepwise elimination was used to identify the most parsimonious model. If a significant 3-way interaction (concentration*phase*time, p < 0.05) or 2-way interaction (concentration*phase, p < 0.05) was observed, differences between concentration groups were tested in each phase using t-tests with a Tukey-Kramer adjustment for multiple comparisons.
For targeted chemistry data analysis, multiple linear regression models were used to identify significant predictors of E2 tissue concentration (pmole/larva) or E2 exposure media concentration (µM). Backwards stepwise elimination was used to identify the most parsimonious model using square root-adjusted values to satisfy modeling assumptions related to normality and homoscedasticity. The effects of concentration, status, or interaction on zebrafish tissue or exposure media concentrations or abundance were assessed (p < 0.05). For zebrafish tissue, any negative values (five non-detects, only found in DMSO controls) were assigned the lowest non-negative value (0.0003 pmol/larva) in the dataset. If a significant interaction between dose and status was observed, pairwise comparisons across groups were evaluated using differences of least squares means and Bonferroni-adjusted p-values (p < 0.05). For non-targeted chemistry analysis, Spearman rank correlation coefficients were calculated to assess the relationships between integrated peak areas for each detected E2 metabolite and measured E2 tissue concentrations from targeted analysis.

Data Availability
The datasets generated during the current study are available by searching for the manuscript title at https:// catalog.data.gov.