Chemical structure drives developmental toxicity of alkyl-substituted naphthalenes in zebrafish

TL;DR

This study shows that the chemical structure of alkyl-substituted naphthalenes affects their toxicity in zebrafish embryos, using a new method to test semi-volatile compounds.

Contribution

A novel vial-based method was developed to assess developmental toxicity of semi-volatile naphthalenes in zebrafish.

Findings

The vial method revealed variable developmental toxicity not seen in standard 96-well plate exposures.

LogKow and body burden were poor predictors of toxicity, suggesting a non-narcotic mode of action.

Transcriptomic analysis indicated disruption of glucocorticoid signaling pathways.

Abstract

Naphthalene and its alkyl-substituted derivatives are among the most abundant polycyclic aromatic hydrocarbons (PAHs) in environmental and human exposure studies, yet their developmental toxicity and mode of action remain poorly understood due to challenges in testing semi-volatile compounds. This study developed a vial based, high throughput method to effectively assess the activity of naphthalenes and a set of 24 alkyl-substituted naphthalenes. Early life stage zebrafish were exposed to a concentration series of each chemical (0–50 μM) in rotating sealed glass vials to minimize volatilization. Benchmark concentration (BMC50) values were calculated for morphological endpoints and lowest effect levels were determined for behavioral effects. The data were assessed for evidence of a narcotic mode of action using body burden measurements for select chemicals and logKow modeling. Targeted transcriptomics at a single concentration and timepoint as well as in silico molecular docking were conducted to generate mode of action hypotheses. The vial method enabled detection of highly variable developmental toxicity not previously observed using standard 96-well plate exposures. LogKow and body burden were poor predictors of toxicity, suggesting a non-narcotic mode of action. Transcriptomic analysis revealed evidence for the disruption of glucocorticoid signaling pathways. Molecular docking identified potential protein targets (e.g., CYP1A2, NT5E, FOLR1) that may mediate observed effects. This study demonstrates the importance of appropriate exposure methods for semi-volatile compounds, reveals structure-dependent toxicity among alkyl-substituted naphthalenes, and provides a foundation for further mechanistic studies and improved risk assessment of alkyl-substituted PAHs.

Figures

B](#F4)). Using body burden concentration instead of logK_ow_ as a proxy for membrane partitioning, we further tested the validity of the narcosis mode of action hypothesis. Body burden concentration at 24 hpf did not correlate with the estimated percent Any effect at 20 μM (from concentration-response curves described in Section 2.3.2; Spearman’s rank correlation p-value = 0.23 and ρ value of 0.7. The body burden of 2,4,5-trimethylnaphthalene was less than half that of 1,5- and 2,7-dimethylnaphthalene but had a much higher percent effect at 20 μM than either of the latter two chemicals ([Fig.

B](#F5)). *Fkbp5* was the most frequent DEG, associated with 18 of the 22 chemicals and with an average log_2_FC of −3.57. Additionally, *fkbp5* expression was negatively correlated with increasing potency determined by Spearman’s rank correlation analysis (Table S4). *Fkbp5* expression was strongly correlated with BMC_50_ with a rho of 0.63 and an adjusted p-value of 0.01. *Fkbp5* is part of the FK506 binding protein family and an immunophilin with peptidylprolyl-isomerase activity. It is known to play critical rolls in protein folding, immune regulation and is central to glucocorticoid recep

Full text

Introduction

Polycyclic aromatic hydrocarbons (PAHs) are a ubiquitous class of environmental contaminants characterized by two or more fused aromatic rings. PAHs are in high abundance in petroleum products and coal and are also produced through the incomplete combustion of organic materials like tobacco and wood (Casillas et al., 2025). While PAHs are naturally occurring chemicals, their presence in the environment is largely attributable to anthropogenic activities (Howsam & Jones, 1998). Naphthalene is the smallest of the PAH chemical class, comprised of two fused aromatic rings. Naphthalene has numerous industrial applications. Isolated from petroleum and coal tar, it is used in the production of polyvinyl chloride (PVC) plastics, dyes, resins, and insect repellent (Preuss et al., 2003). As a member of the US Environmental Protection Agency’s 16 priority PAH list, naphthalene has been included in many environmental monitoring studies (Hudson-Hanley et al., 2021; Keith, 2015; Souza et al., 2023). One study, using data collected from the National Health and Nutrition Examination Survey, measured the presence of oxygenated metabolites of four parent PAHs, including naphthalene in urine (Hudson-Hanley et al., 2021). Metabolites of naphthalene were 20–30 % more prevalent than the other PAHs measured and increased over the fourteen-year study (Hudson-Hanley et al., 2021). Environmental monitoring has routinely found that naphthalene is the most abundant PAH in air and water samples (Moran et al., 2023; Yu et al., 2021). While alkyl-substituted PAHs are much less commonly monitored, when measured, they are routinely more abundant than parent or otherwise substituted PAHs (Golzadeh et al., 2021; Moradi et al., 2022; Qiao et al., 2020; Wnorowski et al., 2022). Alkyl-substituted naphthalenes are typically the most abundant amongst alkyl PAHs (Peng et al., 2023; Wnorowski et al., 2022). For example, comprehensive review of PAH derivatives found that 2,3- and 2,6-dimethylnaphthalene alone made up 10 % of all alkyl PAHs measured (Peng et al., 2023).

There is a large body of knowledge concerning the toxicity of naphthalene, while research on alkyl-substituted naphthalenes, beyond 1- and 2-methylnaphthalene, is limited. Naphthalene was designated by the National Toxicology Program (NTP) as reasonably anticipated to be a human carcinogen in 2004 thus, carcinogenicity is typically the main concern for human health. Naphthalene can also cause respiratory and hepatic injury in a club cell specific manner (Chen et al., 2018; Van Winkle et al., 1997). The mechanisms underlying these effects are linked to the in-situ production of toxic metabolites (eg. Naphthalene 1,2-epoxide) by non-AHR induced cytochrome P450s (CYPs) (Casillas et al., 2025; Genter et al., 2006; Warren et al., 1982). Monomethylated naphthalenes are singly alkylated with similar production sources and uses and are associated with similar health effects as unsubstituted naphthalene (Casillas et al., 2025). The mechanisms underlying these effects are also linked to toxic metabolite related oxidative stress. Though ring oxidation leading to epoxide formation does occur with these chemicals, oxidation at the methyl group is the primary metabolic pathway (Casillas et al., 2025; Wang et al., 2020). Thus, any differences in toxicity between naphthalene and monomethylated naphthalenes are often attributed to the reduction in the formation of reactive epoxides (Höke & Zellerhoff, 1998). However, evidence suggests that 2-methylnaphthalene may also act through non-CYP dependent mechanisms, as CYP inhibition is not as effective in reducing toxicity as it is for naphthalene (Lin et al., 2009).

There is limited research on the toxicity of alkyl-substituted naphthalene compounds beyond monomethylated derivatives. In mice, 2-isopropylnaphthalene and 2,6-diidopropylnaphthalene are minimally active compared to naphthalene and its monomethylated derivatives (Honda et al., 1990). This difference in activity was attributed to the further reduction in ring oxidation due to the increased alkyl substitution of these compounds, and minimal formation of reactive epoxides (Höke & Zellerhoff, 1998; Honda et al., 1987; Kojima et al., 1982). Metabolite formation of a number of other alkylated naphthalenes beyond those previously mentioned, was conduct by Wang et al. (2020), who found that alkylated derivatives of naphthalene favored side chain oxidation over ring oxidation with negligible ring oxidation on naphthalenes with six or more alkyl substitutions (Wang et al., 2020). The authors suggested that the reduction in ring oxidation of more alkylated naphthalenes would equate to a reduction in bioactivation and toxicity (Wang et al., 2020). This hypothesis is yet to be confirmed in vivo. Another frequently proposed hypothesis for the activity of PAHs and alkyl-substituted PAHs more generally is the narcosis model. In it, chemicals act through baseline toxicity induced by basic partitioning into membranes, leading to membrane disruption (Klüver et al., 2016, 2019; Lambert et al., 2022). There is little experimental work that has tested these hypotheses.

The only prior large scale screenings of alkyl-substituted naphthalenes for vertebrate biological activity were our previous publications (Geier et al., 2018a; Morshead et al., 2025; Shankar et al., 2019). In these studies, 24 alkyl-substituted naphthalenes were screened for biological activity in early life stage zebrafish in a high throughput, 96-well plate format (Geier et al., 2018a; Morshead et al., 2025; Shankar et al., 2019). Embryos were statically exposed in 100 μL of exposure medium and sealed at 6 h post fertilization (hpf) until morphological and behavioral data were collected at 24 and 120 hpf (Geier et al., 2018a; Morshead et al., 2025; Shankar et al., 2019). The study showed no morphological effects and only minimal behavioral effects (Geier et al., 2018a; Morshead et al., 2025). Chemical fate analysis and modeling in the 96-well plate system determined that volatilization was the primary explanation for the minimal uptake of naphthalene by embryos (Geier et al., 2018b; Rude et al., 2025). The challenge of volatilization is not unique to naphthalenes in our previous exposure system, and has been documented with other semi-volatile and volatile compounds in diverse aqueous exposure scenarios (Jochum et al., 2024; Proença et al., 2021; Stadnicka-Michalak et al., 2021). We hypothesized that a more tightly sealed approach to developmental exposures would reduce volatilization and present a truer picture of alkyl naphthalene bioactivity (Geier et al., 2018a; Morshead et al., 2025).

Here we tested our hypothesis in rotating 2 mL glass vials with lined screw caps filled with 2 mL of exposure medium. The technique limited headspace and sealed much more effectively than a 96 well plate could be sealed. We screened naphthalene and the 24 alkyl-substituted derivatives tested in our previous 96-well plate studies for developmental toxicity in zebrafish. The data were used to assess whether increased ring oxidation equated to reduced toxicity, as previously proposed (Wang et al., 2020). We assessed whether alkyl-substituted naphthalenes acted through narcosis by assessing how well logK_ow_ correlated with the observed potency, as previous narcosis models have done (Di Toro et al., 2000; Klüver et al., 2016, 2019). To support our modeling work, body burden concentrations were measured in a subset of chemicals and compared to bioactivity. Targeted transcriptomics and reverse molecular docking data generated hypotheses for the mechanisms underlying the variability in potency of bioactivity. Transcriptomic responses were assessed for gene expression correlation with morphological effects. Reverse molecular docking to a library of ~30,000 predicted human protein binding sites enabled hypotheses about molecular initiating events. Proteins of interest were further assessed utilizing crystal structures and zebrafish protein structure predictions.

Methods

All data analysis, unless otherwise stated, was performed in R studio (R version 4.4.2) (R Core Team, 2023).

Chemicals

2.1.

Analytical grade standards were obtained from the Chemical Mixtures Core of the Oregon State University Superfund Center from various vendors (Table S1). All standards were analytically verified and dissolved in dimethyl sulfoxide (DMSO) if obtained neat or solvent exchanged to DMSO (10 mM concentration). A detailed chemical list with CASRN and stock concentration can be found in Supplemental Information 1 (SI-T1).

Zebrafish husbandry

2.2.

Adult tropical 5D line zebrafish (Danio rerio) were maintained under Institutional Animal Care and Use Committee protocol 2021–0227 at the Oregon State University, Sinnhuber Aquatic Research Laboratory (SARL, Corvallis OR). Fish were kept in densities of ~ 500 fish/35-gallon tank at 28 °C in a 14 h:10-hour light:dark cycle, tank water was recirculated, filtered water was supplemented with Instant Ocean salts. Adult fish were fed twice a day with Sparos Zebrafeed 300 μm and spawned by placing a specialized gridded platform and funnel in the tank the night before, spawning was initiated upon first light in the morning. The resulting embryos were collected, screened for malformed or unfertilized individuals, and sorted to the same developmental stage into large glass petri dishes containing embryo medium (EM) (Kimmel et al., 1995).

Toxicity testing

2.3.

When embryos were 6 hpf they were bleached in groups of 1000 or less to limit bacterial growth in sealed vials. Embryos were placed in a strainer for transfer between solutions in 500 mL containers. First embryos were placed in system fish water for 1 min, then transferred to a 2.25 % sodium hypochlorite solution for five min, then transferred to a 70 ng/mL sodium thiosulfate solution for 1 min and rinsed in EM for five min and transferred to a clean glass petri dish. Pools of eight embryos were placed into 2 mL clear glass screw cap exposure vials with polyvinyl-faced pulp liner (Avantar Radnor, PA, United States) using a glass Pasteur pipette. Excess EM was removed, and 2 mL of fresh EM was immediately added. A total of 10 uL of 200x working stock (in DMSO) was added to each vial, then immediately capped and inverted three times to mix. No precipitation was observed for any chemicals at the concentrations tested. Vials were loaded into a custom 3D printed vial holder for slow end-over-end rotation in a dark 28 °C incubator overnight. Rotation of the vials served to reduce chemical partitioning into the headspace and ensure uniform respiration of embryos (Bedell et al., 2025). At 24 hpf, embryos were removed from exposure vials and singulated into 96-well plates loaded with 100 μL chemical-free EM. Each exposure vial group was maintained in plate columns. After loading into 96-well plates, embryos were assessed for 24-hour mortality. For the remainder of the exposure, plates were kept stationary at 28 °C in the dark, until 120 hpf when morphology screening and the larval photo motor response (LPR) assay were conducted, section 2.3.2 and section 2.3.3, respectively. 96-well plates were labeled with a unique plate identifier and all data collected for the endpoints explained in the following sections were connected through a custom laboratory information management system (Zebrafish Acquisition and Analysis Program (ZAAP)) (Truong et al., 2016). A visual summary of this exposure method is shown in Fig. 1.

Experimental design

2.3.1.

All chemicals were tested at five concentrations and a control (0, 3.125, 6.25, 12.5, 25, 50 μM) and normalized to 0.5 % DMSO, a concentration that did not cause morphological effects but fostered solubility (Hoyberghs et al., 2021). Each treatment group was replicated over four vials, resulting in 32 embryos per treatment, split over two replicate 96-well plates per chemical. Some chemicals lacked intermediate response data (between background and 100 % effect) which hindered our ability to accurately model the sigmoidal concentration-response curve and estimate benchmark concentration 50 % (BMC_50_) values. For these chemicals, follow-up exposures were conducted at four intermediate concentrations between background and 100 % effect to improve sigmoidal concentration–response modeling. These concentrations were 12.8, 16, 20, and 25 μM for chemicals lacking data between 12.5 and 25 μM and 30, 37.5, 45, and 50 μM for chemicals lacking data between 25 and 50 μM. Vial groups of plated larvae were tracked to assess the impact of the vial on the results. All reported concentrations were nominal water concentrations and are included in SI-T4.

Morphology

2.3.2.

At 24 hpf embryos were assessed for mortality (MO24). At 120 hpf embryos were assessed for abnormality in nine morphological endpoints: mortality (MORT), edemas (EDEM), bent axis (AXIS), touch response (TCHR), craniofacial (CRAN), muscular/cardiovascular (MUSC), lower trunk (LRTK), brain (BRAN), skin (SKIN), and notochord malformations (NC). Morphological assessments were made by a team of three to four individuals. The team was blinded to the identity, plate layout, and concentrations of the chemicals, thereby minimizing bias in morphological screening assessment. To ensure the assessments were consistent, concentration–response reproducibility of a parathion ethyl quality control plate was performed along with each experimental exposure set up. Details of all measured endpoints can be found at Developmental toxicity zebrafish phenotype atlas. The binning of morphological endpoints for assessment was done as previously described in (Wilson et al., 2022). Each endpoint was scored as a binary presence/absence and scores were compiled into an “Any effect” endpoint, the presence of abnormality in any individual endpoint.

Vial effect assessment.

2.3.2.1.

Morphological results were assessed for the impact of vial exposure groups on morphological outcomes. This assessment used generalized linear mixed effects models and variance source partitioning for each chemical using lme4::glmer() function in R (Bates et al., 2015). Scaled concentration (log-plus-one scaled) was used as the fixed effect to predict the binomial presence or absence of the Any effect endpoint and the vial group of a given observation was assigned as the random effect. R^2^ values for the model fixed effect (concentration) and random effect (vial group) for each model were obtained using the partR2::partR2() function, remaining error was calculated as 1-fixed R^2^-random R^2^. The ratio of fixed R^2^:random R^2^ was calculated, a ratio > 1 indicated that the impact of the vial grouping explained less of the variability observed in the data than the concentration, thus it was appropriate to treat embryos as individual observations. For chemicals with a fixed R^2^:random R^2^ ratio < 1 the ratio of random R^2^:remaining error was used to determine the magnitude of variability explained by the vial groupings. Low random R^2^:remaining error ratios (<1) indicated that the amount of variability explained by the vial groupings was negligible in explaining the overall observed variability in the data and thus it was appropriate to treat embryos as individual observations for the following analysis. The result of this analysis is available in SI-T2. All chemical data sets met the criteria to treat embryos as individual observations for all downstream analysis.

Morphology summary statistics.

2.3.2.2.

The Any effect percent was used to calculate the BMC_50_, 50 % change relative to the background for each endpoint and chemical. All values were reported as the nominal exposure concentration. BMC_50_ values were calculated as previously described (Truong et al., 2022). Briefly, a unrestricted three-parameter log-logistic model for dichotomous data with “extra risk” was fit following guidelines described in the EPA BMDS v3.2 manual, using the maximum likelihood estimation to find the three parameters (US Environmental Protection Agency (EPA), 2020). These curves were used to calculate the reported BMC_50_ values which were required to be within the exposure concentration range, 95 % confidence intervals (BMC upper (BMCu) and BMC lower (BMCl)) were computed using the likelihood-ratio method described in the EPA BMDS v3.2 manual (US Environmental Protection Agency (EPA), 2020). BMC values were reported if the BMCu:BMCl ratio was < 40. BMC_50_, BMCl and BMCu and binomial response data are available in SI-T3 and SI-T4, respectively.

LogKow modeling.

2.3.2.3.

Average predicted logK_ow_ values for each chemical were obtained from the Comptox dashboard (version 2.5.2) when available (Williams et al., 2017). The logK_ow_ value for 2,4,5-trimethylnaphthalene was unavailable on the Comptox dashboard and was predicted using the KOWWIN module in EPI suite^™^Estimation Programs Interface (US EPA, 2024). The BMC_50_ values for chemicals with calculated values for the Any effect endpoint and predicted logK_ow_ values were used to build a simple linear model in R. Model performance was assessed using p-value, R^2^, and average relative error. The average relative error was calculated for each chemical in the model using the equation, |residual|/observed BMD50 and averaged. Model parameters are available in Table S2 and predicted logK_ow_ values are included in SI-T5.

Larval photomotor response (LPR)

2.3.3.

The LPR assay was conducted at 120 hpf. This was done in 96-well plates using a light–dark cycle and video tracking software via the Viewpoint ZebraBox behavioral system software version 5.29.0.70 (Viewpoint Life Sciences, Lyon, France). Larval movement was tracked throughout minutes 6–24. Animals dead or morphologically abnormal were excluded from the analysis, and if more than 30 % of the wells in a treatment group were excluded, the treatment was removed.

LPR summary statistics.

2.3.3.1.

Animal removal due to mortality and morphology effects precluded estimation of concentration–response curves, instead the lowest effect level (LEL) was used. Movement values were binned into 6s intervals, log normalized and used to calculate area under the curve (AUC). AUC values from the first cycle after calibration were used to calculate the LEL for each period (light and dark). AUC values were compared to the day-matched control using a Student’s t-test (Holm-corrected p value ≤ 0.05). An LEL value was only reported where the next highest concentration also showed a significant effect in the same direction, except for the highest concentration. This minimized the risk of false positives. LEL values and individual animal movement data are available in SI-T6 and SI-T7 respectively.

Body burden measurements

2.4.

Body burden was measured for exposure to five chemicals at two timepoints: 24 and 48 hpf. These chemicals (1,5-dimethylnaphthalene, 1,7-dimethylnaphthalene, 2,7-dimethylnaphthalene, 2,4,5-dimethylnaphthalene, and naphthalene) were chosen to span a range of chemical toxicity based on morphological BMC_50_ values and to draw comparisons between chemicals with similar structures and variable toxicity. The exposure concentration was chosen as it was the maximum possible concentration that ensured enough viable animals for collection at 48 hpf and allowed direct comparison of chemical uptake at the same exposure concentration.

Tissue collection

2.4.1.

Chemicals were exposed at 20 μM following the exposure procedure detailed in section 1.3, 24 vials per chemical. At 24 hpf, embryos were pooled in groups of 32 embryos (combining embryos from 4 vials) into 1.5 mL polypropylene RINO-safe lock tubes (n = 3 per chemical treatment, n = 3 per control). Tubes were briefly placed on ice to reduce embryo movement (<1 min); the remaining medium was removed; the tube was then flash frozen in liquid nitrogen and stored at −20 °C until extraction. Embryos in the remaining vials were plated as previously described in section 1.3 until collection at 48 hpf following the same procedure as 24 hpf collections (n = 3 per chemical treatment, n = 1 control).

Tissue extraction

2.4.2.

Tissue extraction was performed as previously reported (Rude et al., 2025). Each sample was spiked with 25 μL internal standard in ethyl acetate and 100 μL of 1 mm glass beads (Next Advance, Troy, NY United States) were added to the RINO-safe lock tubes containing thawed tissue samples and left on ice for 10 min 250 mg of Na_2_SO_4_ and 200 μL of ethyl acetate were added to each tube and vortexed. Tubes were homogenized at speed 8 for 6 min in a Bullet Blender (Next Advance, Troy, NY, United States). 300 μL of ethyl acetate was added to each tube and the homogenization step was repeated until no pieces of non-homogenized fish tissue remained. Samples were then centrifuged for 10 min at 4 °C, 1600g. Sample supernatant (0.5 mL) was collected into amber vials and stored at 4 °C until chemical analysis.

Chemical analysis

2.4.3.

Samples were analyzed for quantitation and identification of 161 PAHs. The method was performed with an Agilent 7890 gas chromatograph (GC), with a 7000C triple quadrupole mass spectrometer (MS/MS) using an Agilent J&W PAH select (30 m × 250 μm × 0.15 μm GC column). Calibration for each analyte utilized at least a 4-point (range: 4–7, median 6 points, n = 161) curve with R^2^ ≥ 0.99. Specific instrument conditions were modified from Anderson et al. (2015) and are detailed in Table S3. GC–MS/MS data was analyzed using MassHunter Quantitative Analysis v.B.06.00 SP1 build 6.0.388.1 software (Agilent Corp. Wilmington,DE) (Anderson et al., 2015). Dynamic muti-reaction monitoring mode was used during analyses, and all PAHs were identified using analytical standards obtained from multiple sources (Table S1). Collison energies were optimized for each individual unsubstituted PAH, alkylation series (C1, C2, etc.), and when possible, specific alkylation group (i.e. methyl, ethyl, etc.).

A standard containing all 161 target analytes was used for continuing calibration verification (CCV) and was analyzed before and after each batch of samples with additional checks for system performance approximately every 10 samples. The data quality objective of at least 80 % of the 161 PAHs (i.e. 129 compounds) within 30 % of the expected concentration was met for all CCVs. In addition, specific PAHs used for exposures were all within 30 % of the known value. Solvent blanks were analyzed after every fifth sample and after each CCV, and all target PAHs were below detection limits. PAHs not used for exposures were below detection limits in all samples. To demonstrate precision, one extract was analyzed in duplicate; duplicate analysis displayed an average relative percent difference (RPD) of 0.7 % across detected PAHs. Finally, as a check against potential matrix influence regarding accuracy, one sample extract was over-spiked with a solution containing all 161 PAHs, and the resulting average percent recovery was 104 % with a range of 78–123 %. Resulting values are available in SI-T8.

Statistical analysis

2.4.4.

Body burden concentrations based on average embryo wet weight at each timepoint were determined by the equation below. Total measured chemical was calculated from ng/mL concentrations times sample volume (0.5 mL).

$$[eqn]$$

Average body weights at 24 hpf and 48 hpf were assumed to be 0.31 mg and 0.26 mg, respectively as determined previously (Rericha et al., 2024). Statistical analysis to compare time points within chemicals was performed using a two-sided Student’s t-test (p-value ≤ 0.05). An ANOVA with a follow-up Dunn’s test was used to compare body burden between chemicals at 24 hpf (p-value ≤ 0.05). Body burden concentrations at 24 hpf were compared to the estimated percent effect of each measured chemical at the exposure concentration used for body burden analysis (20 μM). These estimates were made using the fitted concentration-response curves previously described in section 2.3.2.2.

The concentration uptake ratio of naphthalene at 24 and 48 hpf (nmol/embryo)/(nmol of nominal chemical exposure), was compared to the 96-well exposure method using previously published results (Geier et al., 2018b; Rude et al., 2025). Concentration uptake ratios between vial exposures and 96-well plate exposures were compared using a two-sided Student’s t-test, p-values ≤ 0.05 were considered significant. Concentration data as well as metadata highlighting differences in experimental design are included in SI-T9.

Targeted transcriptomics

2.5.

Targeted transcriptomics were performed using the Biospyder^™^ TempO-Seq^™^ 96-well plate platform. This platform used detector oligos for a gene set of 3112 probes selected through data driven approaches and expert curation (Balik-Meisner et al., 2019). The targeted transcriptomics approach was selected to maximize the number of chemicals that could be included in the analysis while managing cost and data analysis requirements. Targeted transcriptomics were limited to 96 samples, three chemicals without morphological effects were excluded from the analysis to ensure n = 4 per chemical and n = 6 controls, excluded chemicals are indicated in SI-T3.

Exposures and RNA collection

2.5.1.

A single exposure concentration of 20 μM was chosen as the maximum possible concentration to elicit a biological response that also ensured enough viable animals for RNA collection for every chemical. The same concentration per chemical was chosen as opposed to a phenotypically anchored concentration due to variability in 48 hpf survival across chemicals at the same 120 hpf percent effect and to allow for a direct comparison of chemical potency and transcriptomic response. 48 hpf was chosen to ensure the measurement of transcriptomic responses prior to the emergence of observable phenotypic effects. This timepoint has been used extensively and enhanced the ability to compare our results across datasets. Exposure was conducted over two days, chemicals were split into day groups by numerical order of randomly generated internal labels, 11 chemicals and control each day, n = 48 embryos across six vials per treatment. Exposure was conducted as previously described in section 2.3. At 48 hpf embryos were collected for RNA extraction. Embryos from two exposure vials were kept for morphological effect assessments which were as expected based on screening results. From the remaining wells, six phenotypically normal embryos were collected and pooled in their vial groups into 1.5 mL Eppendorf tubes, placed on ice (~30 sec), excess medium was removed, placed back on ice (~30 sec), and 200 μL RNA shield was added. Five samples of n = 6 pooled embryos were collected from each treatment.

RNA extraction and quantification

2.5.2.

RNA was extracted from exposure groups on day one and day two. 100 μL of lysis buffer was added to each tube containing pooled embryo samples and 200 μL of RNA shield, the full contents of the tube were then transferred to the homogenization plate (2 mL deep 96-well plate pre-filled with 400 μL 1 mm zirconium beads and 500 μL of lysis buffer). The plate was homogenized for 1.5 min using the Mini-BeadBeater 96 (BioSpec Products, Bartlesville, OK, United States) and centrifuged for five min at 3000 rcf. 550 μL of the supernatant was distributed over two plates, each pre-filled with 275 μL of ethanol. RNA extraction was performed using the KingFisher Apex System (ThermoFisher Scientific, Waltham, MA, United States) and Quick-RNA Magbead kit (Zymo Research Corporation, Irvine, CA, Untied States). Briefly, this automated process uses Mag beads to collect RNA from the sample, after a series of washes and treatment with DNase the RNA is eluted in 50 μL of ultra-pure water. RNA quality was assessed using the 4150 TapeStation System (Agilent, Santa Clara, CA, United States). All samples used had RNA integrity number (RIN) scores >9.8. RNA concentration was determined using Quant-iT RNA Reagent and assay kit and Varioskan ALF multi-mode plate reader (ThermoFisher Scientific, Waltham, MA, United States) in duplicate, samples with high variability (coefficient of variance (CV) > 10) between reads were repeated in triplicate. The average of the quantification reads was used for RNA concentration normalization. Of the five replicates, samples with the highest concentration and lowest CV were selected (n = 4 per chemical treatment n = 3 per day match control). All 94 samples were normalized using ultra-pure water to the lowest concentration sample − 14.9 ng/μL.

Library preparation and sequencing

2.5.3.

1 μL per sample of normalized RNA or ultra-pure water for two blank wells, was used as input for library preparation following the manufacturers protocol for the 96-sample TempO-Seq S1500 + Zebrafish Surrogate Assay panel (BioSpyder Technologies, Inc., Carlsbad, CA, United States). Briefly, samples were added to the assay plate with an equal volume of 2x enhanced lysis buffer and processed through annealing with detector oligos, digestion, and ligation. Proper amplification was verified using real-time PCR green fluorescence. Library purification was performed using the NucleoSpin Gel and PCR Cleanup Kit (Macherey-Nagel Inc., Allantown, PA, United States) following adjustments detailed in the manufacturer’s instructions. PCR amplification products were pooled into a single sequencing library of 5 μL each into a 1.5 mL Eppendorf tube. The pooled and purified libraries were sequenced using one lane of the NovaseqX 10B flow cell (BioSpyder Technologies, Inc., Carlsbad, CA, United States).

Data analysis

2.5.4.

Pre-processing.

2.5.4.1.

Demultiplexing and alignment was performed using the manufacturer’s platform TempO-SeqR (BioSpyder Technologies, Inc., Carlsbad, CA, United States). The output table of gene counts was used for downstream analysis. Total reads per sample ranged from 7.4 million to 16.8 million with a mean of 10.6 million reads. Sample naming metadata and raw read counts are available in SI-T10 and SI-T11, respectively. All samples passed quality control metrics where sequence depth > 300 K, >40 % reads aligned to probe, >30 % of reads aligned uniquely to the probes, >300 K aligned reads, >30 % of the probes had at least five reads, per base sequence quality was high, and per base N content was low. Counts of genes with multiple probes were summed as instructed by the manufacturer resulting in 3069 unique genes (Everett et al., 2022). Read counts were normalized by upper quartile normalization using the calcNormFactors function from the Bioconductor package in EdgeR with the 75th percentile of reads per samples as the scaling factor. This normalization was chosen following recommendations for TempO-Seq data (Bushel et al., 2020). Principal component analysis (PCA) of normalized read counts identified a control sample from day one as an outlier which was removed from future analysis. PCA was used to determine that a batch effect was present based on exposure and RNA extraction day (nested by experimental design). The batch effect was removed using the ComBat_seq function from the SVA R package as previously applied to data of this type (Lee et al., 2024). Briefly, this function applied a negative binomial regression model to estimate and remove batch effects from count data (Zhang et al., 2020). This function worked with unnormalized read counts and produced batch corrected but not normalized counts. Upper quartile normalization previously described was applied following batch correction. Normalized and batch corrected read count data are available in SI-T12. A PCA was used to check the results of the batch correction (Fig. S1). Targeted transcriptomic data are available through NCBI’s Gene Expression Omnibus accession number GSE296410 (Edgar et al., 2002).

Differential expression analysis.

2.5.4.2.

Differential expression analysis between controls was performed to determine if a global control would be appropriate. Differential expression analysis was done using the limma package by fitting linear models with empirical Bayes smoothing. P-value was adjusted using the Benjamini-Hochberg method and the adjusted p-value was used for all analysis (Ritchie et al., 2015). There were no differentially expressed genes (DEGs) between the day-matched controls (p-value ≤ 0.05) therefore, a global control was used for analysis. DEGs for the chemical exposures were calculated using the method described above compared to global controls (p-value ≤ 0.05 and log_2_ fold change (log_2_ FC) ≥ |1|) (SI-T13).

Correlation analysis.

2.5.4.3.

Correlation analysis was performed to examine the relationship between the log_2_FC of the top five most frequently differentially expressed genes (DEGs in at least 25 % chemical exposures) and morphological BMC_50_ values. This was done using a Spearman’s rank correlation analysis log_2_FC values for the five genes meeting these criteria and the BMC_50_ values (chemicals without a BMC_50_ value were assigned 50 μM). Spearman’s rank correlation was used as opposed to Pearson’s to include chemicals without a calculated BMC_50_ value. Correlations were considered significant if Holm-adjusted p-value ≤ 0.05 and rho (ρ) ≥ 0.5, results available in Table S4.

Reverse molecular docking

2.6.

Naphthalene and alkyl-substituted naphthalene structures were docked to predicted binding sites of AlphaFold2 predicted structures (Sim et al., 2023; Varadi et al., 2022). Description of the generation of the binding pocket list used in this screen was previously reported (Kong et al., 2024). Briefly, binding pocket predictions were acquired from the HPRoteome-Bsite database (Sim et al., 2023) which contains binding pockets for 16,554 proteins covering 80 % of the UniProt entries in the AlphaFold human proteome database (Varadi et al., 2022). While many proteins had one predicted domain, some had multiple. For each predicted protein domain, both the highest-ranking sequence-based and 3D-structure-based predictions were used in the screen, resulting in 33,446 total binding pockets. Reverse molecular docking screening and validation using zebrafish ortholog sequences were conducted with Autodock Vina v. 1.2.3. Naphthalene and alkyl-naphthalene ligand structures were generated by conversion of smiles sequences obtained from ChemDraw version 21.0.0 (Revvity Signals Software, Inc.) to pdbqt files using Open Babel (O’Boyle et al., 2011). For each protein, the value from the binding pocket and prediction method with the lowest binding energy value was selected for further analysis (all binding energy values available in SI-T14). For each chemical, the binding energy values for each protein were used to create a rank with 1 being the protein with the lowest binding energy value for each respective chemical. These ranks were compared as opposed to predicted binding energy values due to variability in mean binding energy of each chemical.

Follow up docking of NT5E was conducted using a zebrafish derived AlphaFold structural prediction and two experimentally derived structures of human NT5E. The two 3D crystal structures were obtained through the Protein Databank (PDB, www.rcsb.org); PDBIDs were 4H2I and 7P9N. Preparation of the structures for docking was done using Chimera. Ligands and solvents were removed during docking prep with default settings including metal ion charges (Zn^2+^, Ca^2+^, and Cl^1+^) (Pettersen et al., 2004). Docking was conducted using the AutoDock Vina plug-in. The docking boxes were chosen to encompass the area containing the experimental ligand (7P9N centroid: 9.575, 11.878, 36.011 size: 30 × 30 ×30 Å; 4H2I centroid: −18.924,21.495, −30.287, size: 18 ×18 × 16 Å).

The zebrafish-derived protein prediction (zf-Nt5e) was obtained from the AlphaFold protein structure database (alphafold.ebi.ac.uk/download). The binding site location was predicted using ProteinsPlus (proteins.plus) and adenosine monophosphate (AMP) as the representative ligand (Centroid: 3.699, −1.229, 2.578, size: 20 × 20 × 20 Å) (Schöning-Stierand et al., 2022). Docking prep and docking were conducted following the same procedures as the two crystal structures. Follow-up docking scores can be found in SI-T15.

Results and discussion

Morphology and Behavior

3.1.

The vial exposure method revealed morphological effects from alkyl-substituted naphthalene exposure that were not detected from 96-well static exposure screening in previous studies (Geier et al., 2018a; Morshead et al., 2025; Shankar et al., 2019). Vials were more effective for this set of semi-volatile compounds, and we would expect the same for other semi-volatile compounds. Morphological effects in controls were well below the 20 % quality control threshold used in our 96-well exposure method. Seven of the chemicals tested did not have an Any effect BMC_50_ value and those that did had values ranging from 46.6 μM to 13.8 μM (Fig. 1).

The heatmap in Fig. 2 shows the BMC_50_ values for each chemical across all endpoints which had a BMC_50_ value for at least one chemical. MO24 and MORT clustered with Any effect while the non-mortality endpoints grouped together indicating that the Any effect endpoint was driven by mortality for many of the chemicals. Chemical clustering based on the phenotypic profiles revealed five distinct groupings. Cluster one contained only 1,7-dimethylnaphthalene, which had the lowest BMC_50_ value for Any effect of the dimethyl-substituted naphthalenes tested. This chemical had a unique morphological profile in that it was the only chemical with a BMC_50_ value for the MUSC endpoint and was not associated with mortality. Cluster two contained nine chemicals which had the lowest response, with the majority having no BMC_50_ value for any endpoint. The two chemicals in this cluster with BMC_50_ values for the Any effect endpoint were close to the maximum concentration tested (50 μM). Cluster three contained two chemicals, 1,5-dimethylnaphthanlene and 2,7-di-tert-butylnaphthalene which both had no MO24 BMC_50_ and BMC_50_ values for variable endpoints at 120 hpf. Cluster four chemicals (1,4,6,7-tetramethylnphthalene, 1,6,7-trimethylnaphthalene, and 2-butylnaphthalene) had the most diverse morphological effects with BMC_50_ values for all endpoints except MUSC. Cluster five contained 10 chemicals, all of which had BMC_50_ values for Any effect, MO24 and MORT and varying presence of activity in other endpoints (CRAN, AXIS, EDEM, and LTRK). 2-ethylnaphthalene specifically had less potent effects over the same endpoints as those in cluster 4.

Unsubstituted naphthalene as well as the monomethylated 1-methylnaphthalene and 2-methylnaphthalene all had no morphological BMC_50_ values for any morphological endpoints. While the majority of more highly alkyl-substituted naphthalenes caused morphological effects. Notably, 2,6-diisopropylnaphthalene did not; it was one of few alkyl-substituted naphthalenes previously tested in mice, and shown to cause minimal damage to lung tissue (Honda et al., 1990). The lack of activity for this chemical may indicate that reduction in ring oxidation in alkyl substituted naphthalenes reduces their toxicity (Höke and Zellerhoff, 1998; Lin et al., 2009; Wang et al., 2020). While our data for 2,6-diisopropylnaphthalene agrees with previously reported toxicity, our more extensive screen identified many alkyl substituted naphthalenes with potent toxicity. This suggests that lower ring oxidation among higher alkyl substituted species may not correlate with reduced toxicity or that their toxicity is not mediated through ring oxidation metabolites.

Chemicals with behavioral effects (bolded chemical names in Fig. 2) did not follow morphological effect groupings. Most LPR effects were hypoactivity in the dark period with increasing LPR effects at higher concentrations, apart from 1,2-dimethylnaphthalene which had the opposite trend (Fig. 3). One chemical, 1,4,5-trimethylnaphthalene, caused hyperactivity during the light period, which was only observed at the highest concentration (Fig. 3). Larval movement tracks over time for chemicals with LPR LEL values are shown in Fig. S2. Nine of the chemicals were previously screened in the 96-well plate format where all but 2,6-dimethylnaphthalene caused LPR effects (Geier et al., 2018a). Here seven of those previously screened were not active in the LPR assay, but 1-methynaphthalene and 1,2-methylnaphthalene were. This unexpected result was likely due to the more stringent criteria for behavioral activity used here which required that multiple consecutive concentrations have significant LPR effects for estimation of an LEL.

LogKow modeling

3.2.

Predicted logK_ow_ values for all 25 chemicals and their BMC_50_ values were used to build a simple linear model for the BMC_50_ concentrations (Fig. S3). These models are commonly used to determine if adverse effects are due to narcosis (Di Toro et al., 2000; Klüver et al., 2016, 2019). Narcosis generally lowers biological activity via narcotic disruption of cell membrane functions; the potency of a narcotic increases as membrane partitioning increases (estimated using logK_ow_). It is important to note that chemicals without an estimated BMC_50_ value for Any effect were not used to build either of the models presented. The R-squared value for model A, containing all chemicals with a BMC_50_ value, was 0.21, indicating that logK_ow_ explained 21 % of the variability observed in the BMC_50_ values. Model A did not meet the significance criteria with p = 0.053 (Fig. S3A). The predicted logK_ow_ of 2,7-di-tert-butylnaphthalene was identified as an outlier by interquartile range. Removal of 2,7-di-tert-butylnaphthalene as linear model B (Fig. S3B) improved the model to an R-squared and p-value of 0.61 and 0.00021, respectively. The model B predicted logK_ow_ explained 61 % of the variability observed in the data. Additionally, the average relative error (the average percent error between the predicted and observed BMC_50_ values) was 17.3 %, driven by many chemicals having similar predicted logK_ow_ around 4.3 and highly variable Any effect BMC_50_ values (Fig. S3). The domain of applicability of model B is limited as it can only be used for alkyl-substituted naphthalenes with logK_ow_ between 4.23 and 5.13 for which there is an Any effect BMD_50_ value below 50 μM. These limitations make the model impractical as a risk assessment tool. Several chemicals with no Any effect BMC_50_ and therefore not used in modeling, had high logK_ow_ > 5 which further emphasized the poor correlation between chemical potency and logK_ow_. We concluded that logK_ow_ did not predict alkyl-substituted naphthalene potency and that narcosis likely does not drive their biological activity.

Body burden

3.3.

Body burden measurement of select chemicals was also used to test the narcosis hypothesis and understand alkyl-substituted naphthalene fate. Measurements were made at 24 and 48 hpf for five chemicals selected for variable potency and structural comparisons. Exposures were done at 20 μM, the maximum with enough viable embryos for collection at 48 hpf.

Body burden was significantly reduced for all chemicals from 24 to 48 hpf (Fig. 4A). In previous 96-well plate exposures, the body burden of naphthalene reached its maximum around 24 hpf. Though we did not have earlier time points to confirm, it was likely that body burden peaked around 24 hpf in vial exposures as well. In vial exposures the body burden loss by 48 h was 46 %-98 %, while in 96-well plate exposures most of the reduction in body burden occurred between 48 and 72 hpf, i.e. once metabolism became active (Rude et al., 2025). An earlier and more rapid reduction with the vial method was likely due to embryo removal from the chemical exposure at 24 hpf and singulation into chemical-free medium in plate wells.

The concentration uptake ratio of naphthalene (nmol/embryo per nmol available in EM) at 24 hpf in the vial exposures was 40x higher than in the 96-well exposure method 0.1303 vs. 0.0032 (Fig. S4). By 48 hpf the uptake ratio was slightly higher than the plate exposures of Geier et al. (2018b) (0.0013) but lower than the plate exposures of Rude et al. (2025) (0.0032) despite the large reduction between 24 and 48 hpf in the vial exposures (Fig. S4) (Geier et al., 2018b; Rude et al., 2025). This demonstrated the effectiveness of the vial exposure method in reducing chemical volatilization and increasing exposure concentrations.

While all chemicals followed the same trend over time, the body burden concentrations were highly variable at 24 hpf ranging from 2.10–10.6 nmol/mg bw. 1,7-dimethylnaphthalene had the highest concentration at 24 hpf followed by 1,5- and 2,7-dimethylnaphthalene, then 2,4,5-trimethylnaphthalene and naphthalene (Fig. 4B). Using body burden concentration instead of logK_ow_ as a proxy for membrane partitioning, we further tested the validity of the narcosis mode of action hypothesis. Body burden concentration at 24 hpf did not correlate with the estimated percent Any effect at 20 μM (from concentration-response curves described in Section 2.3.2; Spearman’s rank correlation p-value = 0.23 and ρ value of 0.7. The body burden of 2,4,5-trimethylnaphthalene was less than half that of 1,5- and 2,7-dimethylnaphthalene but had a much higher percent effect at 20 μM than either of the latter two chemicals (Fig. 4B). 2,4,5-trimethylnaphthalene had the lowest percent reduction in body burden between 24 and 48 hpf at 46 %. The moderate increase in 2,4,5-trimethylnaphthalene tissue persistence may have increased its toxicity but does not explain its toxicity relative to its reduced body burden. While 1,5- and 2,7-dimethylnaphthalene did have lower body burdens than 1,7-dimethylnaphthalene, this did not correspond to the differences in estimated percent Any effect (13 % and 12 % vs 52 % respectively). Incongruence between body burden and estimated percent Any effect further supported a non-narcotic mode of action.

Targeted transcriptomics

3.4.

Transcriptomic data showed that separation along exposure/extraction day and batch correction, described in section 2.5.4.1, was successful in removing the apparent batch-based separation (Fig. S1). Based on PCA of batch corrected and normalized read counts, the overall data variability was small. Principal components one and two explained only 4.84 % of the variability. This suggested that few transcriptomic changes were captured in association with the chemical exposures. This may be due to the 20 μM concentration chosen to ensure sufficient larvae at 48 hpf. Based on our concentration-response curves, we predict a maximum effect of ~ 50 % for the most potent chemical at 120 hpf at 20 μM. When sampling RNA at only one timepoint there is the possibility of missing the strongest transcriptomic responses. It is also possible that disrupted pathways may not be represented in the genes included in the targeted transcriptomics platform.

The number of DEGs (p-value ≤ 0.05 and log_2_FC ≥ |1|) identified for each chemical did not correlate with potency as determined by BMC_50_ (Fig. 5A). Low overall DEGs could explain this result. However, the discordance is not unique to this dataset; similar discord was observed for a diverse set of polyfluorinated alkyl substances (Rericha et al., 2024). Notably, 1,7-dimethylnaphthalene had substantially more DEGs than the rest of the chemicals tested even compared to chemicals with a similar BMC_50_ like 1,4,5,7-tetramethylnaphthalene (Fig. 5A). This is especially interesting due to the morphologically distinct effects of 1,7-dimethylnaphthalene, which together may indicate a distinct mode of action for this chemical.

Similarities in the structures of these chemicals suggested that at least some might have the same mode(s) of action. To investigate this, we looked for common DEGs across the chemicals. Five genes were differentially expressed in association with six or more chemicals (ctss2.1, fgbp2b, klf9, klf11a, and fkbp5) (Fig. 5B). Fkbp5 was the most frequent DEG, associated with 18 of the 22 chemicals and with an average log_2_FC of −3.57. Additionally, fkbp5 expression was negatively correlated with increasing potency determined by Spearman’s rank correlation analysis (Table S4). Fkbp5 expression was strongly correlated with BMC_50_ with a rho of 0.63 and an adjusted p-value of 0.01. Fkbp5 is part of the FK506 binding protein family and an immunophilin with peptidylprolyl-isomerase activity. It is known to play critical rolls in protein folding, immune regulation and is central to glucocorticoid receptor activity (Eachus et al., 2023; Rein, 2020). Specifically, fkpb5 is a co-chaperone protein for the glucocorticoid receptor, modulating the affinity of cortisol for this receptor (Schiene-Fischer & Yu, 2001). In zebrafish fkbp5 is downregulated by fluoxetine and sertraline (selective serotonin reuptake inhibitors), supporting its roll in glucocorticoid receptor activity in zebrafish (Park et al., 2012). Sertraline treatment of zebrafish leads to hypoactivity in the dark phase of a photomotor assay, like that observed for some chemicals in our study (Yang et al., 2021). Additionally, FKBP5 is clinically associated with atrial fibrillation (AF) and Fkbp5 knock out (KO) mice exhibit an increase in AF occurrence (X. Wang et al., 2023). FKBP5 is a co-chaperone of HSP90′s roll in the degradation of HIF-1a and down regulation of fkbp5 may stabilize HIF-1a in cardiomyocytes resulting in downstream dysregulation of Ca^2+^ waves.

The next two most frequent DEGs were both in the KLF family (klf9 and klf11a) which are key regulators in embryogenesis and highly conserved between zebrafish and mammals (Nagai et al., 2009; Xue et al., 2015). KLF9 and 11 are members of the same group within the KLF family that interact with the sin3A protein to remodel chromatin and modify histones (Nagai et al., 2009; Xue et al., 2015). Both genes had negative log_2_FCs when differentially expressed and klf9 had the largest average log_2_FC of the frequent DEGs at −4.3. Klf11 is downregulated in cancers and KO mice exhibit downregulation of oxidative stress scavengers (Fernandez-Zapico et al., 2003). In zebrafish klf9 expression fluctuates diurnally and KO was found to increase liver size, decrease oxygen consumption rate (OCR), and upregulate glycolytic genes (Drepanos et al., 2023). Notably, Klf9 plays a key role in downstream transcriptomic response from glucocorticoid receptor activity and down regulates the expression of fkbp5 in zebrafish (Gans et al., 2021).

Both fgfbp2b and ctss2.1 were differentially expressed by six of the chemicals. Fgfbp2b encodes for one of five zebrafish fibroblast growth factor binding proteins critical for growth factor signaling during embryogenesis and development (Li et al., 2018). Ctss2.1 encodes a lysosomal cysteine protein important in immune function and removal of damaged proteins under oxidative stress (Fu et al., 2020; Lalmanach et al., 2020; Lu et al., 2024).

Investigation of frequent DEGs revealed commonly shared transcriptomic responses across many of the alkylated naphthalenes. Correlation between potency and fold change of expression of fkbp5 suggested a relationship between expression of this gene and alkyl substituted naphthalene toxicity. The key involvement of klf9 and fkbp5 in glucocorticoid receptor activity suggested disruption of glucocorticoid signaling by alkyl naphthalene exposure.

Molecular docking

3.5.

Zebrafish phenotypes and transcriptomics provide important evidence for understanding how alkyl naphthalenes affect development but are insufficient evidence to generate molecular initiation event hypotheses. To this end, we used a reverse molecular docking screen of our alkyl naphthalene chemical structures against a large library predicted human protein binding pockets (Sim et al., 2023). Given that alkyl-substituted naphthalenes may undergo metabolic activation molecular docking of their metabolites would be highly informative. While the measurement or prediction of the metabolites for the 24 chemicals included was outside the scope of this study, this would be a valuable future direction for this work. Using binding pockets predicted from human sequences versus zebrafish may limit interpretation but offsets the considerable uncertainty in structural predictions of zebrafish proteins due to the limited experimentally derived protein structures and the lack of an available database of predicted binding domains for zebrafish protein structures. Given that over 70 % of human protein coding genes have an easily identifiable ortholog in the zebrafish genome, predicted human binding pockets should be reasonably similar to zebrafish binding pocket structures (Howe et al., 2013). Incorporation of human protein binding domains supports our objective to better understand the impact of alkyl naphthalene exposure on human biology and provides a valuable contribution to the field on its own. Follow-up docking was conducted using a zebrafish sequence-derived protein structure prediction to assess whether binding affinity trends were consistent across species for the protein of interest. NT5E docking using experimentally derived human protein structures was also used to validate the results from the predicted structures in the screen.

Reverse molecular docking screen

3.5.1.

Naphthalene and the 24 alkyl-substituted naphthalenes were docked against ~ 30,000 predicted binding pockets from 16,954 AlphaFold predicted protein structures (Sim et al., 2023). The heatmap in Fig. 6 shows trends in binding affinities across chemicals and proteins. The predicted binding affinity ranks were similar across chemicals apart from 2,7-di-tert-butylnaphthalene, which was associated with many high-ranking proteins in the top 10 of few or no other chemicals, like OPN5 and CD1B.

Of all the CYP proteins included in the screen, CYP1A2 showed the highest binding affinity for most of the alkyl substituted naphthalenes. Previously, in vitro studies found that CYP1A2 had the highest affinity of the CYPs for naphthalene (Cho et al., 2006). Recapitulation of this through reverse molecular docking demonstrated the potential to describe biologically significant interactions. While toxic metabolites from CYP activity are hypothesized to be the cause of naphthalene and monomethylated naphthalene toxicity, previous in vivo work demonstrated that naphthalene carcinogenicity was not due to AHR inducible CYPs specifically like CYP1A2 as naphthalenes are poor AHR ligands (Genter et al., 2006). This suggested that although naphthalenes may have a high affinity for binding CYP1A2 may not be induced in vivo and thus not contribute to the observed toxicity.

All three Folate acid receptors (FOLR1, FOLR2, and FOLR3) were frequently in the top binding affinities. FOLR1 specifically was in the top 10 predicted binding affinities for all 25 chemicals. Folate receptors are primarily responsible for the transport of folate (Vitamin B) into cells. Folate is essential for cell growth and normal development. In humans, FOLR1 is the main folate receptor associated with this activity. In zebrafish there is one folate receptor (Folr) which is most similar in structure and activity to the human FOLR1 protein (Jones et al., 2017). The folate pathway is tightly regulated during development and its disruption is linked to numerous developmental defects in zebrafish and humans (M. S. Lee et al., 2012; Rosenquist & Finnell, 2001). If alkyl-substituted naphthalenes can bind Folr in vivo, it could inhibit cross membrane folate transport and be a potential toxic mechanism that warrants further investigation.

NT5E was among the top three predicted binding affinities for 22 of the chemicals and was predicted to have the highest binding affinity for 15. NT5E encodes protein CD73, an extracellular membrane dimerized protein that converts extracellular AMP to adenosine (Jackson, 2011; Minor et al., 2019). NT5E has been implicated in heart failure and as an immune suppressor with a role in cancer progression (Kordaß et al., 2018; H. Wang et al., 2024). Additionally, CD73 plays an important role in the cAMP pathway linked to the regulation of glucocorticoid receptor signaling through extracellular signal-regulated kinase (ERK) (Jackson, 2011; Kang et al., 2023). CD73 has a zebrafish ortholog by the same name, expressed at 24 hpf in developing zebrafish embryos, making it plausible that these chemicals may interact with it during the exposure window (Thisse & Thisse, 2004). We saw nt5e expression, at low read counts, in many of the chemical exposures in our targeted transcriptomics, however, it was detected in only one control sample (SI-T11). While nt5e was not significantly differentially expressed in these exposures, likely due to the low read count, we did see decreased expression of another 5′-nucleotidase, nt5dc2, in four chemical exposures (1,7 and 2,7-dimethylnaphthalene, 2-ethylnaphthalene, and 1,4,5-trimethylnaphthalene). NT5DC2 is a less well described protein than NT5E, but it also contains a 5′-nucleotidase sequence and is suspected to be involved in adenosine production and/or degradation (Cros-Perrial & Jordheim, 2025). This suggested a potential for disruption of purinergic signaling by our exposures. This, as well as the predicted high affinity for the NT5E protein in silico, supported the binding of nt5e and disruption of purinergic signaling by alkyl-substituted naphthalenes as a potential mechanism of the observed glucocorticoid related transcriptomic response.

Follow-up molecular docking

3.5.2.

Based on our findings from the reverse molecular docking screen using predicted human protein binding sites, we conducted follow-up docking studies on alkyl-substituted naphthalenes and NT5E. All chemicals were docked against the predicted structure for NT5E based on the zebrafish sequence (zf-Nt5e) as well as previously published crystal structures for the open and closed forms of the NT5E (7P9N and 4H2I respectively). Zf-Nt5e was determined to be in the closed conformation based on overlay with both crystal structures. The predicted binding pocket was similarly located to the binding pocket for the closed conformation crystal structure 4H2I. Adenosine and AMP were docked to both the screening binding pocket and follow-up structures to compare the predicted binding affinities of the alkyl-substituted naphthalenes to those of the endogenous ligands. Results from this follow-up docking are shown in Fig. 7.

The follow-up docking scores for all structures recapitulated the trends seen in the screening data across chemicals, with a few notable exceptions; 2,7-ditertbutylnaphthalene, 1,4-dimethylnaphthalene, and 2,6-dimethylnaphthalene. All three of these chemicals had much lower binding affinities in the screening than they did in the follow-up docking. Most of the chemicals tested had higher binding affinities for all NT5E structures than for adenosine and AMP. This was more pronounced in the open conformation human crystal structure 7P9N, where the endogenous ligands had lower binding affinities. The two AlphaFold predicted structures had higher binding affinities across all chemicals, including the endogenous ligands, when compared to the experimentally derived structures, apart from the three chemicals previously mentioned. These results further confirm the findings from the reverse molecular docking screen which suggested that NT5E was a target for alkyl-substituted naphthalenes that could outcompete its endogenous ligand. Collectively, while these findings suggest that alkyl-substituted naphthalenes have a strong affinity for NT5E, additional enzymatic and genetic experiments are required to confirm its involvement in the toxicity of these chemicals in zebrafish and human model systems.

Conclusions

The vial exposure method developed in this study captured morphological and behavioral effects from developmental exposure to alkyl-substituted naphthalenes by overcoming volatility losses inherent in with 96-well plates. The improved exposure protocol enables more accurate hazard assessment of semi-volatile contaminants, with implications for improved risk assessment. Our results showed high variability in activity across structurally similar compounds and 1,7-dimethylnaphthalene as having distinct morphological effects. A primary hypothesis for naphthalene toxicity involves oxidative damage induced by toxic metabolites requiring aromatic ring oxidation products. Contrary to a previous hypothesis that increased alkyl substitution would lead to decreased toxicity due to decreased ring oxidation, naphthalene and its monomethylated derivatives were not overtly toxic while most of the more highly alkyl substituted derivatives were. Metabolite analysis could confirm the extent of in vivo alkyl substituted naphthalene ring oxidation, while measurement of oxidative stress would further clarify the role of this mechanism in the observed toxicity. Conducting these experiments across time, concentration, and structural variants would provide the most insight. This approach would enable the identification of temporal trends in metabolite formation and stress responses, link morphological effects to metabolites and oxidative stress, and reveal the structural features that drive metabolite formation. Variable potency of the alkyl substituted naphthalenes was not explained by logK_ow_ predictions or body burden concentrations, suggesting a non-narcotic mode of action.

Despite low variability in the transcriptomic responses across chemicals, targeted transcriptomics suggested a shared mode of action by many of the chemicals involving glucocorticoid signaling. Measuring cortisol concentrations across time, concentration, and chemical potency could confirm these impacts to glucocorticoid signaling. Reverse molecular docking identified potential protein targets and molecular initiating events. Predicted high-affinity interactions with proteins such as CYP1A2, NT5E, and FOLR1 suggested possible mechanisms underlying the transcriptomic and phenotypic outcomes. Recapitulation of binding affinity trends for NT5E in follow-up docking further supported the hypothesis that alkyl substituted naphthalenes act through disruption of purinergic signaling. Future work will focus on validating the predicted interaction between alkyl substituted naphthalenes and NT5E using a cross-species approach. High-throughput approaches will be used to directly evaluate the binding affinity and functional effects of each alkyl-substituted naphthalene on human NT5E and zebrafish Nt5e (Freundlieb et al., 2014; Senger et al., 2004), and complementary genetic experiments in zebrafish will definitively determine the role of Nt5e in mediating the toxicity of these compounds.

Supplementary Material

1

2