Zebrafish (Danio rerio) was raised according to standard protocols. The following zebrafish lines were used: AB wild-type (wt) strain, and transgenic Tg(eef1a1l1:EGFP) expressing enhanced green fluorescent protein (GFP) in neuron cells. Zebrafish embryos were obtained by natural spawning, collected within 30 min after fertilization, and cultured at 28.5°C (Tian et al., 2019). To evaluate the toxicity of ESIONPs, a dose–response analysis was carried out to determine the median lethal dose (LC50). The zebrafish embryos at 4 h post-fertilization (hpf) were distributed in 6-well plates (30 for each group), and exposed to ESIONPs suspensions at different concentrations (0 mg/L, 10 mg/L, 20 mg/L, 30 mg/L, 40 mg/L, 60 mg/L, 80 mg/L, and 100 mg/L). The medium was changed every 24 h. The survival rate was determined every day by counting the embryos that survived. The exposed embryos were collected at indicated stages for different analysis. Embryos from each group were observed and photographed taken an SMZ25 stereomicroscope with a DS-Ri2 digital camera (Nikon, Japan). All experimental procedures on zebrafish were approved by the Institutional Animal Care and Use Committee of Northwest University and carried out in accordance with the approved guidelines (NWU-AWC-20190103Z).

Embryos in 40 mg/L ESIONP-exposed group and control group at 48 hpf, 72 hpf, 96hpf, and 120 hpf were collected and used for total RNA extraction. mRNA was isolated using the NEBNext PolyA mRNA Magnetic Isolation Module (New England Biolabs, Ipswich, MA, United States). Libraries were prepared with the NEB Next Ultra Directional RNA Library Prep Kit (New England Biolabs, United States), and subjected to Illumina sequencing with paired end 2 × 150 as the sequencing mode. The clean reads were mapped to reference genome (D. rerio: NCBI_GRCz11). Gene expression levels were estimated using FPKM (fragments per kilobase of exon per million fragments mapped) by StringTie v1.3.4d (Pertea et al., 2015). Differential expressed genes (DEG) were measured using R package, edgeR v3.24.2 (Robinson et al., 2010). The false discovery rate (FDR) was used to calculate the adjusted p-value in multiple testing in order to evaluate the significance of the differences. Here, only gene with an adjusted q-value < 0.05 and |log2FC| ≥ 1 were used for subsequent analysis. The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (Chen et al., 2021) in National Genomics Data Center (CNCB-NGDC Members and Partners, 2024), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA016266) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.

A weighted gene co-expression network analysis was performed using the WCGNA package in R (Langfelder and Horvath, 2008). Samples were clustered by hclust to filter outliers (h > 15). In order to construct scale-free network, the optimal soft-thresholding power β was defined by picking Soft Threshold function (β = 6, R ² ≥ 0.8). Based on pairwise correlations between genes, genes with similar expression patterns were clustered into a group through a TOM clustering tree according to the dynamic tree cut method, and similar groups were combined into one module. The key hub genes, which were the node of co-expression network, were defined based on the connectivity by the CytoHubba plugin in Cytoscape v3.9.1.

Gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were carried on DEGs (Gene Ontology Consortium, 2021; Kanehisa et al., 2023). Enrichment analysis was performed using the R package “clusterProfiler” (Wu et al., 2021). GO terms and KEGG analysis with corrected p-value < 0.05 were considered to be significantly enriched (Wang et al., 2021).

Total RNA was isolated from zebrafish embryos in each group using TRIzol™ reagent (Ambion, United States). The cDNA was synthesized using the SuperScriptIII (Invitrogen, United States), according to the manufacturer’s protocol. qPCR was conducted using SYBR FAST Universal qPCR kit (KAPA, Germany) and ViiA 7 Real-Time PCR System (ABI, United States), as described previously (Wang et al., 2021). Primer sequences were shown in Supplementary Table S1.

Zebrafish embryos at 96 hpf were anesthetized with tricaine (0.04%) and embedded in 5% methylcellulose to score the skeletal muscle lesions. Birefringence was imaged under SMZ25 stereo microscope equipped with a DS-Ri2 digital camera (Nikon, Japan), as previously described (Lu et al., 2021). For quantification analysis, 10 somites between the levels of somite 5 to 15 were imaged per embryo.

For locomotion tracking, single zebrafish larvae (treated with or without ESIONP) developed to 120 hpf was placed in individual wells of 24-well cell culture plate containing approximately 500 µL embryo medium. Swimming behavior was monitored at room temperature using a DanioVision system and EthoVision XT 11.5 locomotion tracking software (Noldus Information, Netherlands), according previously described (Lu et al., 2021).

The oxidative stress and damage caused by ESIONP was detected by measuring ROS production in zebrafish embryo (Zhu et al., 2022). Briefly, embryos at 72 hpf treated with or without ESIONPs were stained with an oxidation-sensitive fluorescent probe dye, dichloro-dihydro-fluorescein diacetate (DCFH-DA) (Beyotime, China) at a final concentration of 20 μg/mL. Stained embryos were incubated at 28°C for 1 h and then washed with PBS. The photos were taken under a fluorescence microscope with a DS-Ri2 digital camera (Nikon, Japan). The fluorescence intensity of embryos was quantified using ImageJ software (NIH, United States).

To detect apoptotic cells in zebrafish embryos, acridine orange (AO), a fluorescent dye was used. Zebrafish larvae developed to 96 hpf were incubated with 10 μg/mL AO staining solution (Beyotime, China) at 28.5°C for 30 min in the dark, and rinsed with PBS (Zhu et al., 2022). The zebrafish embryos were observed and recorded under a fluorescence microscope with a DS-Ri2 digital camera (Nikon, Japan). The intensity of the fluorescence signal was measured and analyzed using ImageJ software (NIH, United States).

Each experiment was repeated at least three times. All data were presented as the mean ± SD. Student’s t-test was applied for comparisons among different groups. p-value < 0.05 was considered significant.

To assess the toxicity of ESIONPs, zebrafish embryos were exposed to different concentrations of ESIONPs (0, 10, 20, 30, 40, 60, 80, and 100 mg/L), the survival rate was counted at 24, 48, 72, 96, and 120 hpf (Supplementary Figure S1A). The LC50 of ESIONPs was determined at 72 hpf (Supplementary Figure S1B). The degrees of malformations are defined in 4 classes, including dorsal bending, shortened body length, yolk sac swelling, cardiac edema, smaller eyes, and head (Supplementary Figure S1C). A concentration of 40 mg/L was chosen for subsequent experiments, to ensure the maximum survival rate of embryos while including diverse morphological abnormalities.

In order to explore the mechanism of toxicity of ESIONPs on zebrafish embryonic development, the DEGs of the control groups and the ESIONPs-exposed groups at 48 hpf, 72 hpf, 96 hpf, 120 hpf were analyzed by RNA-Seq. To determine the relationship between replicates, the samples were clustered using Principal Component Analysis (PCA) and correlation analysis. The PCA showed high repeatability in duplicate samples (Figure 1A), and the Pearson correlation coefficient for each group of replicates also indicated a high repeatability (|r| ≥ 0.8 for all) (Figure 1B). A total of 32,057 transcripts from 16 samples were fused to construct the co-expression network which was constructed by R packages of WGCNA. Firstly, samples were clustered by hclust (h > 5,500), no outlier samples were found in the hierarchical clustering (Figure 1C). In order to build a network with scale-free distribution and preserve the information of DEGs as much as possible, we found the best soft-thresholding powers β (β = 6). The connectivity between genes in the network is relatively high (β = 6, R ² = 0.813), indicating that the network was scale-free (Figure 1D). Genes were divided into 22 modules (Figure 1E; Supplementary Table S2). In order to study the mechanism of toxicity of ESIONPs at each development stages, the module with the highest correlation of changes was selected at each stage. The yellow module (|R ²| = 0.62), magenta module (|R ²| = 0.3765), lightgreen module (|R ²| = 0.725), and brown module (|R ²| = 0.6) showed the highest correlation of changes at 48 hpf, 72 hpf, 96 hpf, and 120 hpf, which were used for further research (Figure 1F).

Genes in key modules were screened out according to the eigengene-based connectivity (kME) values (|kME| > 0.9). There were 61 hub genes in the yellow module, 128 hub genes in the magenta module, 54 hub genes in the lightgreen module, and 1,157 hub genes in the brown module (Supplementary Table S3).

In the enrichment analysis of the GO pathways, we presented top 10 GO annotations (Supplementary Figures S2A–D; Supplementary Tables S4–S7). The important terms of enrichment in the yellow module (48 hpf) were related to neuron system and muscle development (Supplementary Figure S2A). The terms of related neuron development also mainly were enriched in the magenta module (72 hpf) (Supplementary Figure S2B). The terms of inflammatory were enriched in the lightgreen module (96 hpf, mainly including inflammatory response and chemotaxis of immunocyte (Supplementary Figure S2C). The terms of neuronal signal transmission were enriched in the brown module (120 hpf), including (Supplementary Figure S2D).

In enrichment analysis of KEGG pathways (Supplementary Figures S2E–H; Supplementary Tables S4–S7), yellow module genes (48 hpf) were mainly enriched in embryonic development key signaling pathways, including notch signaling pathway, wnt signaling pathway, hedgehog signaling pathway (Supplementary Figure S2E). The terms of hormone secretion and metabolism were enriched in the magenta module (72 hpf), including pentose phosphate pathway, ubiquitin mediated proteolysis, protein processing in endoplasmic reticulum, GnRH signaling pathway, Endocytosis (Supplementary Figure S2F). Lightgreen module genes (96 hpf) were enriched in metabolism pathways (Biosynthesis of nucleotide sugars, arachidonic acid metabolism, amino sugar and nucleotide sugar metabolism, glycerophospholipid metabolism), necroptosis, C-type lectin receptor signaling pathway, and ferroptosis (Supplementary Figure S2G). Brown module genes (120 hpf) were enriched in pathways related to neuronal signal transmission such as calcium signaling pathway, neuroactive ligand-receptor interaction, cell adhesion molecules, and oxidative phosphorylation (Supplementary Figure S2H).

The enrichment analysis of modules corresponding to each stage of embryonic development by GO and KEGG indicated that ESIONPs might mainly be toxic to the nervous system development, neural conduction, and motor system of zebrafish, and might induce inflammation and ferroptosis in zebrafish embryos.

To investigate the mechanism of toxicity of ESIONPs on zebrafish embryos, we next filtered out the hub genes affected by ESIONPs. Networks were constructed to explore relationships among hub genes, which were used as nodes of the scale-free network and had the highest correlation. The top hub genes as the most important nodes in each module were identified and highlighted (Figure 2). There were 3 hub genes in the yellow module: kif15, ncapd2, smarcc1a (Figure 2A), 5 hub genes in the magenta module: zc2hc1a, hunk, pbx1a, lypd6, si:ch211-235e15.1 (Figure 2B), 5 hub genes in the lightgreen module: irg1l, si:ch211-153b23.3, si:ch211-153b23.5, zgc:100868, noxo1a (Figure 2C), and 3 hub genes in the brown module: chrna6, slc4a10a, slc6a11b (Figure 2D). Although some hub genes (si:ch211-235e15.1, si:ch211-153b23.3, si:ch211-153b23.5, zgc:100868) of the magenta and lightgreen modules remain unannotated, the other hub genes of each module could reflect the toxicity of ESIONPs to embryos. For example, genes involved in neuron development (kif15, ncapd2, smarcc1a, hunk, pbx1a, lypd6), neurotransmission (chrna6, slc4a10a, slc6a11b), immune system regulation (irg1l), and oxygen stress (noxo1a). For each module, two hub genes were selected for validation by qRT-PCR (Figure 2F), which exhibited similar expression trends to the RNA-Seq profiles (Figure 2E). The WGCNA analysis indicated that ESIONPs might have neurotoxicity, which could damage neuron development, nerve conduction, and synaptic transmission. In addition, ESIONPs might also cause ferroptosis in zebrafish embryos.

In Tg (eef1a1l1:EGFP) transgenic zebrafish embryos, ESIONPs-exposed embryos exhibited significant abnormal development in the nervous system at 72 hpf (Figure 3A), compared with the control embryos. Furthermore, the expression of neuron developmental markers (pax2a, neurog1, axin2) was significantly downregulated in ESIONPs-exposed embryos (Figure 3B). The analysis of the movement track at 120 hpf showed that the movement ability of the ESIONPs-exposed group was significantly weakened (Supplementary Figure S3A), and the expression of the neuromuscular junction and synapse markers (lrp4, musk, mpz) was also downregulated in ESIONPs-exposed embryos (Supplementary Figure S3B). Moreover, the muscle polarization of zebrafish larvae exposed to ESIONPs was significantly reduced (Figure 3C), and the expression of muscle markers (acta2, ttn, lpin1) was also downregulated (Figure 3D). These results showed that ESIONPs not only impaired neuron development, synaptic signal transmission, and neuromuscular junction signal transmission, but also reduced muscle development.

Next, we determined whether ferroptosis occurred by detecting oxidative stress and apoptosis in zebrafish embryos (Dixon et al., 2012). The oxidative stress was significantly increased in the ESIONPs-exposed embryos (Figure 3E), and the expression of oxidative stress markers (cybb, nox1, rac2) was significantly upregulated in ESIONPs-exposed embryos (Figure 3G). The ESIONPs-exposed embryos also displayed significant cell apoptosis (Figure 3F), and the expression of apoptotic markers (jnk1, bcl2a, tp53) was significantly upregulated (Figure 3H). In addition, the expression of inflammatory markers (il1b, il6, tnfa) was also increased significantly (Supplementary Figure S3C). These results revealed that ESIONPs could induce ferroptosis, resulting in oxidative stress, cell apoptosis, and inflammatory response in zebrafish embryos.