Zebrafish husbandry

All experiments involving zebrafish were approved by the Institutional Animal Care and Use Committee at the U.S. EPA National Health and Environmental Effects Research Laboratory and performed in accordance with appropriate guidelines and regulations. A mixed wild type adult zebrafish line (Danio rerio) was used for all studies48. Adult zebrafish were housed in 6 L tanks (at a density of ~ 8 fish/l) and fed Gemma Micro 300 (Skretting) once daily and shell free E-Z Egg (Brine Shrimp Direct) twice daily during the week. On weekends, zebrafish were fed both food sources once daily. Adults were maintained on a 14 h:10 h light:dark cycle at 28.5 °C. For spawning, 60–100 adults were placed in 10 L angled static breeding tanks overnight. The following morning, adults were transferred to new bottom tanks containing treated reverse osmosis water (60 mg/l sodium bicarbonate and 0.4 g/l Crystal Sea Bioassay Formula Marine Mix) and embryos were collected after 45 min.

Chemicals

17β-estradiol (E2) was purchased from Sigma Aldrich (#E8875, CAS: 50-28-2; Lot: SLBP6339V). Stock solution aliquots (40 mM) were prepared by dissolving neat chemical into anhydrous dimethyl sulfoxide (DMSO) and stored at −80 °C. For each 10-day experiment, fresh working solutions were prepared by thawing 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-13C 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.

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.

Axenic derivation and conventionalization with fish facility microbiota

Axenic (AX; microbe-free) zebrafish were generated as previously described19,39,40 (Fig. 1). Briefly, embryos were resuspended in FS-10% HBSS containing antibiotics (amphotericin B (0.25 μg/ml), kanamycin (5 μg/ml) and ampicillin (100 μg/ml)) for four hours at 26 °C. At 5 hours post fertilization (hpf), embryos were treated in a 15 ml conical tube with 0.5% poly(vinylpyrrolidone)-iodine (PVP-I; CASRN 25655-41-8) for 2 min and 0.05% bleach for 20 min and sorted into sterile T25 tissue culture flasks (15–30 embryos per flask, depending on experiment) in 25 ml of FS-10% HBSS. As a control for the derivation process, a subset of axenic embryos was conventionalized with aquaculture facility microbiota at 1 dpf to generate the axenic colonized on day 1 (AC1) zebrafish cohort. Embryos were conventionalized by 80% media change which included 10 ml fish facility water (FRW, fish room water; 5 µm syringe-filtered) and 10 ml of FS-10% HBSS (Fig. S8). To maintain consistency between cohorts, axenic and conventionally colonized (colonized on day 0, not treated with antibiotics, PVP-I, or bleach) flasks also underwent an 80% media change at 1 dpf. Similar to axenic colonized on day 1 zebrafish, conventionally colonized embryos received 10 ml FRW and 10 ml of FS-10% HBSS. Axenic embryos received 10 ml of 0.2 µm filter-sterilized fish facility water (FS-FRW) and 10 ml of FS-10% HBSS (Fig. S8).

For experiments involving all three zebrafish cohorts, media sterility was tested as previously described19. 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 previously48.

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 described19. 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™ DNA/RNA MiniPrep Kit (Zymo Research #D7003). Total DNA yield was quantified and samples were stored at −80 °C. DNA (250 ng) from each sample was added to triplicate PCR reactions along with barcoded primers72 specific for the V4 region of the 16S rRNA gene and amplified with the Roche FastStart High Fidelity PCR System (Sigma-Aldrich, #4738292001). PCR reactions were run at 95 °C for 2 min, followed by 25 cycles of 95 °C for 30 sec, 55 °C for 30 sec and 72 °C for 1 min, with a 10 min final extension at 72 °C. Triplicate reactions were pooled and products were purified and normalized with the SequalPrep Normalization Plate Kit (ThermoFisher, #A1051001). Samples were pooled by volume and DNA was sequenced using the Illumina MiSeq Reagent Kit v2 (500 cycles, #MS-102-2003) and Illumina MiSeq instrument. Positive and negative PCR control reactions were run with every 30 samples19.

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 v773. Reads were analyzed using the QIIME 1.9.0 software package74,75. The 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 predictions76 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 described19. 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 annotations77 and relative abundances of predicted functional pathways were formatted as previously described78. 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 described19,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-13C 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-13C 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-13C 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. Total run time was 30 min. An Agilent 6530B Accurate-Mass Quadrupole Time-of-Flight Mass Spectrometer (Palo Alto, CA) with a Dual AJS ESI source was operated under positive and negative electrospray ionization in full scan mode (100–1000 m/z). The following spray settings were employed: capillary voltage (±3500 V), nozzle voltage (500 V), fragmentor (135 au), skimmer (65 au), octupole RF peak (750 au), gas temperature (300 °C), gas flow (1 L/min), nebulizer pressure (40 psig), sheath gas temperature (350 °C), and sheath gas flow (11 L/min). Data were acquired under 2 GHz extended dynamic range mode.

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 database80. 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.