Transcriptome Analysis of Thiram-Treated Zebrafish (Danio rerio) Embryos Reveals Disruption of Reproduction Signaling Pathways

Simple Summary Thiram is a widely used fungicide, yet studies have shown that it has hazardous effects on animals, including fish. We performed developmental and transcriptome analyses of thiram-exposed zebrafish embryos to further understand thiram toxicity. Thiram-exposed embryos exhibited several developmental abnormalities. Transcriptome analysis revealed 1754 differentially expressed genes associated with thiram exposure when compared to controls. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of thiram-treated samples highlighted reproduction-related pathways. The current findings shed light on the molecular mechanisms through which thiram exerts developmental and reproductive toxicity in zebrafish. Abstract Thiram, a dithiocarbamate fungicide, is used for the treatment of various fungal infections affecting crops and ornamentals. However, thiram-associated toxicity has been reported in animals, including fish, and the underlying molecular mechanisms are unclear. Herein, we employed zebrafish (ZF) to gain further insights into thiram toxicity-associated molecular mechanisms. We studied developmental abnormalities and performed whole-transcriptome analysis of ZF embryos exposed to thiram for 96 h. Embryos exposed to 4.0 μg/L thiram exhibited several phenotypic abnormalities, including bradycardia, spinal curvature, hatching arrest, and growth retardation. Whole-transcriptome analysis revealed 1754 differentially expressed genes (DEGs), with 512 upregulated and 1242 downregulated DEGs. The majority of biological processes affected by thiram were metabolic. Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis yielded terms related to reproduction, such as steroid biosynthesis and steroid hormone biosynthesis. Quantitative real-time polymerase chain reaction validation results were in line with sequencing data for ten DEGs. The study results improve our current understanding of the effects of thiram exposure in ZF.


Introduction
Pesticide toxicity is a major threat to people and ecosystems worldwide. Pesticide usage has considerably increased in the previous decade in a bid to improve agricultural productivity and quality through the control of pests, fungi, and weeds [1,2]. In particular, fungicides are essential for crop protection, as they prevent fungal growth [3]. Their increasing use has led to increased environmental contamination and concerns regarding their effects on the environment as well as on human health [2][3][4]. Dithiocarbamate (DTC) fungicides were shown to be teratogenic, mutagenic, carcinogenic, and neurotoxic in animal and fish studies, possibly sharing a similar mechanism of toxicity [5][6][7][8]. Thiram is a DTC fungicide that is widely used to protect fruits, vegetables, ornamentals, and turf from pests [9]. It is also used as an accelerator and vulcanizing agent in the rubber industry [10]. hatching, with no embryos hatching at 72 h post-fertilization (hpf) and less than 5% embryos hatching at 96 hpf (Table 1), was observed. The developmental anomalies observed at 24,48,72, and 96 hpf are indicated in Figure 1.   ZF exhibit average heartbeats of 150-200 beats per min (bpm) between 30 and 60 hpf at 28 • C [25]. Thiram affected ZF cardiac function at a dose of 4.0 µg/L. DMZF had mean heart rates of 170.4 ± 7.1 (48 hpf), 185.2 ± 9.8 (72 hpf), and 194.2 ± 11.3 (96 hpf) bpm ( Figure 2). LTZF exhibited similar mean heart rates of 172.6 ± 11.8 (48 hpf), 182.8 ± 10.2 (72 hpf), and 193.6 ± 12.4 (96 hpf) bpm. However, HTZF showed significantly lower heart rates than did DTZF (Video S1), at 145.0 ± 16.1 (48 hpf), 157.2 ± 8.4 (72 hpf), and 166.2 ± 11.5 (96 hpf) bpm. There was a high incidence of PE in the HTZF group (Table 1). At 24 hpf, 12.8% of ZF had edema symptoms, with PE increasing to 32.2% at 72 hpf. Hyperemia ( Figure 1) and blood flow abnormalities (data not shown) were also evident in the HTZF group.
Thiram affected the overall growth of ZF ( Figure 3). The mean growth of larvae in the LTZF group was similar to that of the DMZF group larvae, while HTZF had a shorter body length than the DMZF (27.74% shorter). Mean body lengths in the DMZF, LTZF, and HTZF groups at 144 hpf were 3.91 ± 0.39, 3.88 ± 0.29, and 2.82 ± 0.53 mm, respectively. Normal spines were noted in both the DMZF and LTZF groups, as opposed to spinal curvatures in the HTZF group (Table 1 and Figure S1), with wavy distortions of the notochords (Table 1 and Figure 4). Thiram affected the overall growth of ZF ( Figure 3). The mean growth of larvae in the LTZF group was similar to that of the DMZF group larvae, while HTZF had a shorter body length than the DMZF (27.74% shorter). Mean body lengths in the DMZF, LTZF, and HTZF groups at 144 hpf were 3.91 ± 0.39, 3.88 ± 0.29, and 2.82 ± 0.53 mm, respectively. Normal spines were noted in both the DMZF and LTZF groups, as opposed to spinal curvatures in the HTZF group (Table 1 and Figure S1), with wavy distortions of the notochords (Table 1 and Figure 4). Results are presented as mean ± standard deviation (n = 3). Unpaired Student's t-tests were performed to assess significance. * statistically significant (p < 0.05).   Bar graph showing mean heartbeats per minute observed at 48, 72, and 96 hpf. Results are presented as mean ± standard deviation (n = 3). Unpaired Student's t-test was performed for statistical analysis. * statistically significant (p < 0.05).
Thiram affected the overall growth of ZF ( Figure 3). The mean growth of larvae in the LTZF group was similar to that of the DMZF group larvae, while HTZF had a shorter body length than the DMZF (27.74% shorter). Mean body lengths in the DMZF, LTZF, and HTZF groups at 144 hpf were 3.91 ± 0.39, 3.88 ± 0.29, and 2.82 ± 0.53 mm, respectively. Normal spines were noted in both the DMZF and LTZF groups, as opposed to spinal curvatures in the HTZF group (Table 1 and Figure S1), with wavy distortions of the notochords (Table 1 and Figure 4). Results are presented as mean ± standard deviation (n = 3). Unpaired Student's t-tests were performed to assess significance. * statistically significant (p < 0.05). Results are presented as mean ± standard deviation (n = 3). Unpaired Student's t-tests were performed to assess significance. * statistically significant (p < 0.05).

Whole-Transcriptome Analysis
Transcriptome analysis was performed (DMSO (n = 3), thiram 0.4 μg/L (n = 3), and thiram 4.0 μg/L (n = 3), 96 hpf) to comprehensively assess the molecular response to thiram exposure in ZF. Table 2 shows sequencing data from Illumina. Illumina-based sequencing data returned an average of 71,928,169 clean reads for DMSO-treated samples, 78,526,392 clean reads for 0.4 μg/L-treated samples, and 75,297,512 clean reads for 4.0 μg/L-treated samples. Over 98% of clean reads satisfied the Q20 threshold, and over 94% satisfied the Q30 threshold. An average of 86.95% clean reads were unequivocally mapped to the ZF reference genome (GRCz11). A heat map of two-way hierarchical clustering analysis with the z-score obtained from the normalized value (log2-based) is illustrated in Figure 5A. Using the p-value and fold change of each transcript, volcano graphs were constructed to show the overall distribution of transcripts as well as to determine differentially expressed genes (DEGs) under thiram treatment ( Figure 5B,C). Cut-off values for DEG selection were a 2-fold change and a corrected false discovery rate (FDR) p-value of 0.05. A total of

Gene Ontology (GO) Enrichment Analysis
Significant DEGs were classified into three functional categories as per GO enrichment analysis, including biological process (BP), cellular component (CC), and molecular function (MF). Figure 6 depicts the top ten enriched GO terms in each category following high-dose thiram therapy (HTZF group). Metabolic processes, biosynthetic processes, and catabolic processes were among the 10 most enriched terms for "BP" ( Figure 6A). With regard to the top 10 enriched terms for "CC", most were related to the membrane and extracellular  Figure 6B). Most MF terms were associated with oxidoreductase, ion binding, and catalytic activity ( Figure 6C).

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis
DEGs were subjected to KEGG pathway enrichment analysis to determine signaling pathways regulated in response to thiram treatment. A total of 77 signaling pathways were significantly enriched in the HTZF group (p < 0.05). These were categorized into six subgroups. The metabolism subcategory was at the top, with 52 signaling pathways, followed by cellular processes, with 11 signaling pathways. Organismal systems and environmental processing included six pathways each, while human diseases and genetic processing consisted of one pathway each ( Figure S3). A total of 70 signaling pathways were enriched in the LTZF group (p < 0.05), grouped into six categories. The metabolism subcategory was once again the largest, with 44 signaling pathways, followed by cellular processes with eight signaling pathways. The organismal systems, environmental processing, genetic information processing, and human diseases subcategories included 6, 5, 1, and 1 signaling pathways, respectively ( Figure S3).

Quantitative Real-Time Polymerase Chain Reaction (qPCR) Validation of Transcriptome Data
The transcript levels of ten randomly selected DEGs with up-or downregulation tendencies were examined using qPCR. All genes presented the same expression patterns in agreement with RNA-seq data (Table 4).     Figure 6 depicts the top ten enriched GO terms in each category following high-dose thiram therapy (HTZF group). Metabolic processes, biosynthetic processes, and catabolic processes were among the 10 most enriched terms for "BP" ( Figure 6A). With regard to the top 10 enriched terms for "CC", most were related to the membrane and extracellular region ( Figure 6B). Most MF terms were associated with oxidoreductase, ion binding, and catalytic activity ( Figure 6C).

Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis
DEGs were subjected to KEGG pathway enrichment analysis to determine signaling pathways regulated in response to thiram treatment. A total of 77 signaling pathways were significantly enriched in the HTZF group (p < 0.05). These were categorized into six subgroups. The metabolism subcategory was at the top, with 52 signaling pathways, followed by cellular processes, with 11 signaling pathways. Organismal systems and environmental processing included six pathways each, while human diseases and genetic processing consisted of one pathway each ( Figure S3). A total of 70 signaling pathways were enriched in the LTZF group (p < 0.05), grouped into six categories. The metabolism subcategory was once again the largest, with 44 signaling pathways, followed by cellular processes with eight signaling pathways. The organismal systems, environmental processing, genetic information processing, and human diseases subcategories included 6, 5, 1, and 1 signaling pathways, respectively ( Figure S3).

Discussion
Given the extreme toxicity of pesticides, it is crucial to obtain a comprehensive understanding of the molecular mechanisms activated in response to pesticide exposure in off-target organisms. Pesticides exert considerable toxicity, capable of disrupting cellular homeostasis [26]. Herein, we employed next-generation RNA-sequencing technology to explore the response to thiram exposure in ZF at the molecular level. To this end, developmental deformities were scored, and ZF samples were subjected to whole-transcriptome analysis. A wide spectrum of developmental deformities were noted in the HTZF group ( Figure 1), with 1754 DEGs, including 512 upregulated and 1242 downregulated DEGs ( Figure S1). Thus, 4.0 µg/L of thiram considerably disrupted cellular homeostasis, impacting the normal growth and development of ZF.
Thiram has been included in the European Union's Category 1 list of endocrine disruptors based on extensive scientific evidence. Female rats treated with thiram exhibited a disrupted hormonal control of ovulation [27]. Thiram inhibited the rise in luteinizing hormone (LH) at a dose of 50 mg/kg in all rats tested [27]. Chen et al. recently reported that early thiram exposure disrupts endocrine function in ZF [19]. Furthermore, results from our laboratory showed that thiram affected the fecundity of ZF (data not shown). Interestingly, KEGG pathway enrichment analysis of DEGs obtained from the HTZF group revealed a significant enrichment of hormone-related signaling pathways, including steroid biosynthesis, steroid hormone biosynthesis, and retinol metabolism (Figure 7). "Steroid hormone biosynthesis" was the top-ranked enriched pathway as per KEGG results for the HTZF group (Figure 7), indicating that the disruption of steroid hormone homeostasis is the major underlying mechanism of thiram toxicity in ZF. Steroid hormones are required for the initiation and maintenance of secondary sex characteristics, reproduction, and other bodily functions [28]. In the steroid hormone biosynthesis pathway, several UDP-glucuronosyltransferase (UGT) (21 genes) and cytochrome P450 (CYP) (five genes) genes were downregulated in response to thiram. UGT and CYP enzymes both have detoxification activity [29]. More specifically, UGT glucuronidates a variety of endogenous and exogenous substrates, while CYPs oxidize. In addition, UGTs are primarily involved in androgen signaling [30] and CYPs in estrogen signaling [31]. Several important genes, such as srd5a2a (steroid 5-alpha-reductase, alpha-polypeptide 2a), had11b1la (hydroxysteroid (11-beta) dehydrogenase 1-like a), and comta (catechol O methyltransferase a) were also downregulated in the HTZF group. SRD5A2 converts testosterone into the more potent dihydrotestosterone (DHT), an essential hormone in testicular development and spermatogenesis [32]. 11β-hydroxysteroid dehydrogenase is essential in fish spermatogenesis, as it is required for the conversion of 11β-hydroxytestosterone to 11-ketotestosterone (11-KT) and the conversion of cortisol to cortisone [33]. Catechol O methyl transferase is essential for reproduction and metabolizing estrogen as well as epinephrine, catecholamines, and norepinephrine [34]. Sex steroids permanently constitute the brain circuits responsible for reproductive behavior during fetal development, stimulating these specialized pathways when released from the gonads during puberty [35]. For example, the ZF male holding water contains 17 beta-diol, 5 alpha-androstane-3 alpha, and cholesterol, which act as ovulation inducers [36].
Another KEGG pathway enriched for the HTZF group was "steroid biosynthesis" (Figure 7). Human fertility is negatively impacted by genetic abnormalities that manifest as functional changes in proteins of the steroidogenic pathway [37]. Four genes known to play an essential role in cholesterol biosynthesis, namely, methyl sterol monooxygenase 1 (msmo1), lanosterol synthase (lss), emopamil-binding protein (ebp), and 7-dehydrocholesterol-∆7 reductase (DHCR7), were downregulated in the HTZF group (Table 3). Msmo1 is the key enzyme in the post-squalene cholesterol biosynthetic pathway and removes a methyl group from C4-methyl sterols [38]. Lanosterol synthase is required for the cyclization of the first sterol, an essential step in cholesterol biosynthesis [39]. Emopamil-binding protein acts as a sterol isomerase in the pathway [40]. DHCR7 is required to complete the conversion of 7-dehydrocholesterol to cholesterol, as the final step [41]. Furthermore, cholesterol serves as a critical precursor for the biosynthesis of steroid hormones, which in turn regulate reproductive pathways and the development of secondary sexual characteristics, and de novo defects in cholesterol synthesis impact fertility and embryo viability [37].
"Retinol metabolism" was also enriched in the KEGG pathway analysis of HTZF (Figure 7). Retinoids play an essential role in various physiological processes, including reproduction, growth and development, vision, and the immune system [42,43]. Exposure to diethylaminobenzaldehyde (DEAB), an inhibitor of retinoid synthesis, or intake of a low-retinoid diet resulted in ZF producing 70 to 90% fewer eggs [42]. Retinoic acid deprivation at the embryonic stages of ZF development resulted in microphthalmia before birth [42]. Embryos deprived of retinoic acid also developed severe PE, smaller eyes, yolk sac edema, and visual disturbances [42]. As the HTZF group exhibited similar developmental anomalies (Figure 1), it is conceivable that some of these were due to defects in retinol metabolism.
The next protein downregulated in response to thiram was Zondhesin (ZAN). ZAN is a multi-domain sperm protein essential in the binding interaction of sperm and eggs during reproduction [44]. Angiotensin I-converting enzyme 2 (ACE2) was downregulated about 16-fold in the HTZF group. According to Keskus et al., ace2 expression is crucial for ZF embryonic development, with in silico analysis of transcriptome datasets revealing that ace2 expression predominantly increases from 3 days post fertilization (dpf) to 4 dpf [45]; this result provides support for its role in ZF embryonic development.
All fish in the HTZF group displayed wavy notochord distortions and skeletal deformities ( Figure 4). The results are in line with those of Chen et al., who found that ZF treated with thiram exhibited wavy distortions of the notochords [19]. The findings are consistent with those in previous reports of body axis curvatures and notochord distortions in thiram-treated ZF [20]. Inhibition of copper-dependent lysyl oxidases was previously shown to lead to abnormal notochords in ZF embryos [46], while the addition of copper restored the wavy notochord phenotypes in ZF [47]. It is therefore possible that the wavy notochord phenotype seen after thiram treatment ( Figure 5) is due to the synergistic effects of thiram-mediated lysyl oxidase suppression (Table 3) and copper ion chelation. Downregulated mucins (mucin 2.1, 16.76-fold, and mucin 13b, 5.66-fold) could equally explain the skeletal deformities observed in HTZF, since mucin knockout mice have a decreased bone mass and increased osteoblast mineralization activity [48]. The expression of three genes involved in cholesterol biosynthesis (msmo1, lss, and ebp) was downregulated in the HTZF group (Table 3), with deficiencies in genes required for cholesterol biosynthesis inducing bone defects. Moreover, msmo1 and lss are essential in bone formation, and mutations of these genes in prehypertrophic chondrocytes cause severe abnormalities in bone formation due to cholesterol deficiency and intermediate sterol accumulation [49]. Mutation of the ebp gene in humans is associated with craniofacial and skeletal abnormalities, as well as ocular or visual defects [50]. Undoubtedly, further research is needed to better understand the causes of these abnormalities in relation to thiram.
Another group of enzymes that were significantly downregulated in the HTZF group includes chitinase, acidic.1 (chia.1), chitinase, acidic.2 (chia.2), chitinase, acidic.5 (chia.5), and chitin synthase 1 (chs1). Chitin is a structural polysaccharide that is endogenously synthesized in fish, particularly during the early phases of development. Knockdown of chs1 in developing ZF embryos resulted in a significant reduction in chitin signals within the gut [51]. Mammalian chitinases/chitinase-like proteins have been linked to various innate and acquired immunological responses, as well as tissue remodeling and fibrosis induced by both chitinous and non-chitinous stimuli. Chitinase(s) are also known to be involved in normal development of the trunk and tail in ZF [52].
There was a 17-fold downregulation in the expression of meprin A, alpha, tandem duplicate 1 (mep1a.1) in response to thiram. Merpins are expressed in a variety of tissues in ZF, including the gills, gut, liver, skin, kidneys, brain, and heart [53]. According to Schütte et al., meprins plays essential roles in embryonic cell differentiation and proliferation as well as angiogenesis [54]. Meprin 1 knockdown ZF larvae displayed a dilated pericardium or deformed tissues in the trunk and tail due to issues with cell differentiation. Therefore, cardiac abnormalities observed with thiram-treated ZF are possibly due to the reduction of meprin levels.
There was a more than 5-fold downregulation in the expression of nine tRNAs (tRNA-Asp, tRNA-Glu, tRNA-Pro, tRNA-Ala, tRNA-Leu, tRNA-Trp, tRNA-Val, tRNA-Thr, and tRNA-Cys) in the HTZF group compared with that in the DMZF group. tRNAs are assumed to be exclusively involved in translation, as they are originally known for carrying amino acids during protein synthesis. However, recent studies have demonstrated that tRNAs are essential in numerous other processes, such as the control of transcription and translation, post-translational modifications, stress reactions, and diseases [55]. For example, Revendro et al. reported that tRNA dysregulation in ZF causes reduced embryo viability and increased malformations, and postulated that tRNA dysregulation promotes protein aggregation and genomic instability [56]. Reactive oxygen species accumulation, nuclear and mitochondrial DNA damage, and mitochondrial network alterations have also been observed in ZF embryos with tRNA mutations resulting in the incorrect incorporation of serine (Ser) [56]. The expression of several factors known to mediate global protein biosynthesis [57,58], such as translation initiation factors (eIF4A1a and IF21) and termination factors (eTF1a), was also downregulated in the HTZF group (File S1). It is possible that disruptions in protein synthesis were responsible for the growth retardation and other abnormalities seen in the HTZF group ( Figure 3). Interestingly, transcriptome results from the present study (thiram treatment) and our earlier study (phosmet treatment) suggested the disruption in protein synthesis as one of the main mechanisms of toxicity [23]. However, while phosmet decreased the expression of several aminoacyl synthases that ligate tRNAs to their corresponding amino acids, thiram downregulated the expression of several tRNAs, suggesting that pesticides dysregulate protein synthesis through different mechanisms.
The expression of some miRNAs, such as miR-126a (8.4 fold), miR-181a-1, (8.37 fold), and miR-133b (5.03 fold) was downregulated in the HTZF group. MiR-126 is highly enriched in vascular endothelium and is crucial for vasculogenesis, angiogenesis, and endothelial inflammation [59]. MiR-126a silencing impaired the formation of parachordal lymphangiogenic budding from the posterior cardinal vein and thoracic duct of ZF embryos [60]. There is evidence that miR-181 is essential in embryonic development, cell proliferation, autophagy, apoptosis, immune response, and mitochondrial function [61]. The expression of miR-133, crucial in cardiac diseases and development, was decreased during cardiac regeneration, enhancing cardiomyocyte regeneration [62]. In ZF, downregulation of miR-133 expression alters the muscle gene expression pattern and leads to actin disorganization during sarcomere assembly [63].
Early ZF embryonic development is exclusively dependent on the egg yolk for nutrition [64]. Thus, hatching delay or inhibition is deleterious [65]. ZF embryos start hatching at 35 hpf, and almost all hatch at 96 hpf under ideal conditions [66]. Hatch failure was the most significant deformity observed in ZF after treatment with thiram (Table 1). In the HTZF group, almost all embryos remained in their chirons. These results are consistent with findings from previous studies [19,20]. As ZF embryos did not hatch in the presence of other DTCs, such as maneb [67] and ziram [68], it is conceivable that the inhibition of hatching is among the most important toxic effects of DTCs on ZF. The potential to prevent hatching could be explained by either a decrease in the larva's ability to destroy the chorion or by the suppression of proteins required for chorion digestion. The downregulation of expression of tetraspanins (tspan 34 and tspan 13a) (Table 3), required for chorionic digestion in ZF [69], and/or morphological defects, such as tail deformities and spinal curvature, impair the embryo's capacity to rupture the chorion [70,71], possibly causing the hatching failure observed in the HTZF. Because hatching is a critical step of reproduction, hatching delays have a considerable negative impact on populations and ecology [72].
Transcriptome analysis of the LTZF group revealed 138 DEGs compared to the DTZF group ( Figure S2), suggesting that thiram elicits molecular responses even at low concentrations (0.4 µg/L). Despite the lack of an apparent phenotype in the LTZF group (Table 1; Figures 1-4), we cannot completely rule out the possibility that lower-dose exposure may have long-term consequences on ZF development.
Understanding the mechanisms activated during pesticide poisoning in ZF is expected to provide clues for the analogous responses in humans owing to the high degree of genetic homology [73]. The ZF model allows for faster in vivo toxicity evaluation than mammalian tests [74]. Animal toxicological studies often reveal consequences that require additional, costly, and time-consuming research into the exact cause. Combined with transcriptome analysis, toxicity screening approaches employed in ZF allow for the identification of off-target effects as well as for the study of putative underlying mechanisms. ZF models are essential for the identification of the safest candidates for pesticide production and for determining optimal approaches to treat pesticide poisoning. The current work elucidated the molecular response to pesticide toxicity, demonstrating the feasibility of state-of-the-art methods for this purpose.

Ethics Statement
All animal tests were performed in compliance with the standards of the care and use of laboratory animals and were approved by the Animal Ethics Committee of the National Institute of Agricultural Sciences, Rural Development Administration, Republic of Korea (NAS-202102).

ZF Husbandry and Embryo Selection
ZF (Danio rerio, AB strain) were cultured at the fish breeding facility of the National Institute of Agricultural Sciences in South Korea in order to obtain embryos. Fish cultures were maintained as previously described [70,71]. They were kept in a glass tank (50 L) with dechlorinated water and a continuously circulating filter method. The breeding facility was kept at a constant temperature of 25 ± 1 • C and a photoperiod of 14 h light/10 h dark. The fish received one of three diets, including bloodworm (Hikari Bio-pure, Himeji, Japan), live Artemia (INVE Aquakultur, Dendermonde, Belgium), and dry flake diet (Top Meal, Tabia, Korea), at least twice a day. Collection and washing techniques for ZF embryos were previously described [20,23]. Males and females in a 1:1 ratio were mated in the dark for 10 h. Once the light was turned on for 30 min, the embryos were harvested and washed repeatedly with E3 medium to completely remove debris from their surface. Fertilized embryos were selected at 2 hpf and used for the pesticide treatment.

Thiram Treatment and Scoring of Developmental Deformities
Thiram (98.4%, Sigma-Aldrich, St. Louis, MO, USA) was dissolved in DMSO (Sigma-Aldrich). Thiram stocks were aliquoted and stored at -20 • C until use. Test doses were determined according to a previous study [20]. Two test doses were chosen, with the higher thiram dose of 4.0 µg/L being the average EC 50 value for ZF, which would ensure that the most deformities were induced in the collected embryos. The other test dose, 0.4 µg/L, was a 10-fold dilution of the higher test dose. The DMSO concentration in all treatment groups was 0.01%. For the developmental toxicity assay, two 24-well plates, each containing 20 embryos, were used for each test condition. The experiments were repeated three times. Chemical exposure was started after 2.0 hpf, and the experimental conditions for chemical treatment, deformity assessment, and the calculation of percentage deformed embryos were similar to those previously described [20].

Heartbeat Survey
Heartbeat studies were performed under a stereomicroscope (Stemi 508, Zeiss, Germany) at three time intervals (48,72, and 96 hpf) at 26 ± 1 • C. Heartbeats were counted for 20 s, converted to beats per min, and reported. The experiment was performed three times, and each time the heartbeats of 10 embryos were evaluated (n = 30 embryos).

Touch-Evoked Escape Response (TEER)
As previously reported [23], TEER was performed at 144 hpf. The TEER was measured by carefully touching the ZF's head and tail with a thin (1 mm) flexible wire. The experiment was performed three times, and each time, 10 embryos were evaluated (n = 30 embryos).

Body Length Survey
Body length studies were performed at 144 hpf using a stereomicroscope to take pictures. OptiView 3.7 software (Korealabtech, Seongnam, Republic of Korea) was used to calculate body length (from the tip of the mouth to the end of the caudal fin). The experiment was performed three times, and each time, the body length of 10 embryos was evaluated (n = 30 embryos).

Statistical Analysis
An unpaired t-test (Prism 5.0, GraphPad, San Diego, CA, USA) was used to determine significant differences between samples in the analyses of developmental deformities, heart rate, and body length. The p values ≤ 0.05 were considered as statistically significant.

RNA-Sequencing
Samples were collected in triplicate for transcriptome analysis and qPCR validation. Samples containing at least 30 embryos of the 4.0 µg/L thiram-treated, 0.4 µg/L thiramtreated, and DMSO-treated groups were snap-frozen in liquid nitrogen and stored at −80 • C until use. Total RNA extraction, RNA quantification, RNA integrity assessment, RNA library preparation, and cDNA library preparation were performed according to previously described protocols [23]. Paired-end (2 × 100 bp) sequencing was performed on the NovaSeq platform (Illumina, Inc., San Diego, CA, USA).

DEG Selection
By comparing DMSO-and thiram-treated groups, DEGs were selected using nbi-nomWaldTest and DESeq2. The DEG selection criteria were a p-value of 0.05 and a cut-off of the absolute log2-fold change of 2. The acquired DEGs were subjected to further analyses.

Quantitative Reverse Transcription-Polymerase Chain Reaction Validation
The total RNA content was extracted from three samples per test condition using the TRIzol reagent (Bioneer, Daejon, South Korea). The Nanodrop 2000 was used to examine the RNA (Thermo Fisher Scientific, MA, USA). To reverse-transcribe RNA into cDNA, the ReverTra Ace™ qPCR RT master mix (Toyobo, Osaka, Japan) was used. TOPreal TM qPCR 2X premix (Enzynomics, Daejon, Republic of Korea) and CFX96 Dx real-time PCR detection equipment (Biorad, CA, USA) were utilized to perform qPCR. The Ct values of all samples were obtained in triplicates. Average Ct values were used to calculate the gene expression fold change. As a housekeeping gene, β-actin was employed. For the expression analysis, the 2 −∆∆CT approach was used. The oligoprimers utilized for qPCR are listed in Table S2.

Conclusions
As thiram is widely used in agriculture, it is important to understand its toxicity to aquatic life. The current results show that thiram extensively influences gene expression, affecting the growth, reproduction, and behavior of the fish. More specifically, thiram disrupts ZF reproduction by interfering with pathways involved in steroid hormone biosynthesis. These findings also provide substantive insights into the molecular pathways altered in response to thiram and thus causing developmental toxicity in ZF. Further research into thiram disposal in aquatic environments is required.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biology12020156/s1, Figure S1: Representative images showing spinal curves observed in high-dose (4.0 µg/L) thiram-treated zebrafish (HTZF); Scale = 0.5 mm. Figure S2: Numbers of up-and downregulated differentially expressed genes. Yellow bar, upregulated DEGs; blue bar, downregulated DEGs. Figure S3: List of enriched pathways in KEGG pathway enrichment analysis. Table S1: List of top upregulated (FC > 5) and top downregulated (FC < −10) genes after 4.0 µg/L thiram treatment. Table S2: List of primers used for qPCR. Supplementary file 1: Total list of DEGs obtained in high-dose (4.0 µg/L) thiram-treated zebrafish (HTZF) (FC > 2 and FC < −2). Video S1: video showing the heartbeats of DMSO-treated zebrafish embryos and 4.0 µg/L thiram-treated zebrafish embryos. Videos were captured under a stereomicroscope (Stemi 508, Zeiss, Germany) at 96 hpf.  Institutional Review Board Statement: The study was conducted in accordance with guidelines for care and use of laboratory animals as approved by the Animal Ethics Committee of the National Institute of Agricultural Sciences, Rural Development Administration, Republic of Korea (protocol code NAS-202202 (19 January 2022)).