Transcriptome-Based Identification of Genes Responding to the Organophosphate Pesticide Phosmet in Danio rerio

Organophosphate pesticides (OPPs) are one of the most widely used insecticides. OPPs exert their neurotoxic effects by inhibiting acetylcholine esterase (AChE). Most of the gross developmental abnormalities observed in OPP-treated fish, on the other hand, may not be explained solely by AChE inhibition. To understand the overall molecular mechanisms involved in OPP toxicity, we used the zebrafish (ZF) model. We exposed ZF embryos to an OPP, phosmet, for 96 h, and then analyzed developmental abnormalities and performed whole transcriptome analysis. Phenotypic abnormalities, such as bradycardia, spine curvature, and growth retardation, were observed in phosmet-treated ZF (PTZF). Whole transcriptome analysis revealed 2190 differentially expressed genes (DEGs), with 822 and 1368 significantly up-and downregulated genes, respectively. System process and sensory and visual perception were among the top biological pathways affected by phosmet toxicity. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed significant enrichment of metabolic pathways, calcium signaling pathway, regulation of actin cytoskeleton, cardiac muscle contraction, drug metabolism–other enzymes, and phototransduction. Quantitative real-time PCR results of six DEGs agreed with the sequencing data expression profile trend. Our findings provide insights into the consequences of phosmet exposure in ZF, as well as an estimate of the potential risk of OPPs to off-target species.


Introduction
Pesticides improve crop productivity and quality by killing insects, fungi, and weeds; however, their off-target effects on aquatic flora and fauna have raised concerns about pesticide-related human health and environmental risks [1].
Organophosphates are one of the most widely used insecticides. In addition to their use in agriculture, organophosphates are used as chemical warfare agents, and flame retardants [2,3]. In fact, organophosphate pesticides (OPPs) are extremely toxic to humans and account for most self-poisoning deaths in developing countries [4,5]. The consequences of OPPs on humans are primarily linked to neurological disorders [6]. Children were reported to exhibit neurological and developmental abnormalities after being exposed to indoor OPPs [7]. OPPs were reported to induce developmental and behavioral abnormalities in animals, birds, reptiles, amphibians, and fish [8]. Because OPPs are increasingly being used in developing countries [9], it is essential to better understand their off-target effects.
Phosmet is a non-systemic, phthalimide-derived OPP used to control moths, aphids, mites, suckers, and fruit flies in plants and animals [10]. It is widely used on fruit crops, ornamental, cattle, sheep, goats, and pigs [10]. According to the California Department of pesticide regulation, an average of 0.4 million pounds of phosmet was applied in California based on 2000-2018 statistics [11]. Oral ingestion or inhalation of phosmet has been

Effects of Phosmet on ZF Development
Phosmet has been reported to cause various developmental abnormalities in ZF embryos and larvae [19]. Figure 1 shows the developmental deformities of PTZF at 24,48,72, and 96 hpf. A wide range of developmental defects, including abnormal somites, reduced retinal pigmentation, edema, body-axis curvature, and sensory abnormalities, was observed in PTZF at a dose of 8.0 mg/L (Table 1).  The cardiac defects in PTZF ranged from an extremely slow heart rate to subtle changes in blood flow ( Figure 2 and Table 1). The average heart rates of the controls were 163.6 ± 12.2, 185.2 ± 6.3, and 193.2 ± 11.5 beats per minute (bpm) at 48, 72, and 96 hpf, respectively ( Figure 2). PTZF, with 127.8 ± 22.0, 142.8 ± 12.6, and 165.2 ± 11.1 bpm at 48, 72, and 96 hpf, respectively, had a lower heart rate than non-treated ZF (NTZF). Pericardial edema (PE) was also observed in PTZF (Table 1). Approximately 20% of PTZF showed symptoms of PE at 24 hpf, and by 72 hpf, the incidence of PE increased to 62%. In ZF with PEs, reduced blood flow and hyperemia were also frequently observed (Table 1). These findings suggest that phosmet affects cardiac development and function in ZF. NTZF quickly escaped upon touch, whereas PTZF displayed an abnormal touchevoked escape response (TEER) ( Table 1). When PTZF were gently touched on the head or tail, they reacted in a mild to non-reactive manner. Most NTZF quickly identified the object placed in their field of binocular vision and escaped quickly to the opposite corners of the Petri dish. However, PTZF failed to identify the object and the subsequent vision-mediated escape response (VMER) ( Table 1), suggesting abnormal sensory responses in PTZF.
Phosmet decreased the overall growth of ZF (Figure 3). At 144 hpf, PTZF displayed a 20% shorter body length than NTZF. NTZF and PTZF had an average body length of 4.4 ± 0.1 mm and 3.5 ± 0.1 mm, respectively. The graph represents the average body length at 144 hpf. Representative images show the body length at 144 hpf. NTZF, non-treated zebrafish; PTZF, phosmet-treated zebrafish. Scale = 1.0 mm. The results are mean ± SD (n = 4). Statistical significance was evaluated using an unpaired t-test. * (p < 0.001).

Gene Ontology Term Enrichment
Next, we performed gene ontology (GO) enrichment analysis of significant DEGs to functionally classify the genes in three subcategories: cellular component (CC), biological process (BP), and molecular function (MF). Figure 5 shows the top 10 GO terms in each category after phosmet treatment. Significantly enriched terms for "BP" were associated with system process, visual perception, nervous system processes, and sensory perception ( Figure 5A); for "CC," with membrane and extracellular region ( Figure 5B); for "MF," with structural molecule activity, transporter activity, and ion binding ( Figure 5C).

Discussion
Given the high toxicity of OPPs, understanding the molecular mechanisms affecting their non-specific targets is critical. Pesticide-induced toxicity is undoubtedly complicated and can disrupt cellular homeostasis [30]. In this study, we used next-generation RNA sequencing, which is one of the promising techniques for comprehensive transcriptome analysis [31], to translate and understand phosmet-induced toxicity at the molecular level in the ZF model. When comparing the transcriptome of PTZF to that of NTZF, 2,190 DEGs were obtained, with 822 and 1368 significantly up-and downregulated DEGs, respectively (File S1). The identified DEGs were enriched in pathways related to metabolism, muscle development and contraction, ion transport and signaling, and neuron development and signaling ( Figures 5 and 6).
AChE suppression is a major toxic mechanism of OPPs, and AChE measurement is the most frequently used biomarker for environmental pollution with OPPs [30]. AChE inhibition leads to ACh buildup at the synapse, which induces a wide range of behavioral defects in fish, including abnormal swimming and a lack of responsiveness to touch and visual cues [22]. Consistently, decreased AChE expression was observed in PTZF (File S1) and in treatments with other OPPs, including envoy 50 SC [32], chlorpyrifos [33], and diazinon [24], all of which have been reported to induce behavioral and developmental abnormalities in ZF. By contrast, some OPPs, such as dichlorvos, malathion, and methylparathion, did not suppress AChE, but induced deformities in ZF [21]. Richendrfer and Créton revealed that relatively low quantities of OPPs can influence ZF behavior without altering AChE activity [34]. These findings imply that in addition to, or instead of, suppressing AChE, OPPs affect other molecular mechanisms that regulate development and behavior. Based on the obtained transcriptome data (Figures 5 and 6), it is reasonable to speculate that the phosmet can affect many pathways in addition to AChE inhibition. Further research is required to distinguish AChE inhibition from the pathways described.
Gross developmental abnormalities such as curved spines, edema formation, and growth retardation were the most significant in PTZF (Figures 1 and 3, and Table 1). The OPP, dichlorvos, which induces developmental abnormalities in ZF [20], altered several genes involved in energy metabolism and oxidative stress responses [35]. In ZF, metabolism is primarily responsible for development into adult fish [36]; therefore, maintaining metabolic homeostasis is critical. With 151 DEGs, "metabolic pathways" (dre01100) was the top enriched term in KEGG pathway analysis (Figure 6), suggesting that phosmet exposure can affect multiple metabolic pathways in early developmental stages.
The calcium signaling pathway (dre04020) was found to be significantly enriched in PTZF (Figure 6), implying a disruption of Ca 2+ homeostasis. Being a secondary messenger, Ca 2+ plays vital roles in somite formation, eye development, brain partitioning, and heart formation, thus the establishment of the basic vertebrate body plan of ZF [42]. Calcium signaling has also been implicated in early muscle development, neuronal induction, regulation of neural tube closure, cardiogenesis, and phagocytosis at wound sites [43]. Inhibiting calcium signaling during embryonic development induced several developmental abnormalities in ZF, including cyclopia and tail deformities [44]. Consistently, in this study, PTZF exhibited gross developmental and behavioral abnormalities, such as increased incidence of edema and eye defects and decreased response to touch (Table 1). PTZF failed to establish cardiac function (Table 1 and Figure 2), which is consistent with our previous results [19], where phosmet induced a dose-dependent increase in PE, bradycardia, and abnormalities in blood flow. OPPs, such as diazinon [45], dichlorvos [21], and sumithion [20], were also reported as cardiac toxicants. Cardiac problems could arise from defects in heart development and abnormalities in the conduction system or the contractile apparatus [46]. PE is a sign of abnormal cardiac development, whereas heart rate variations are indicators of abnormal cardiac function [47]. Alterations in the rate or pattern of blood flow are an additional sign of abnormal heart function. Our functional analysis revealed enrichment of numerous cardiac-related pathways (Figures 5 and 6) that are crucial for maintaining the structural and functional integrity of the heart, indicating that phosmet exposure affected both heart development and function. Cardiac muscle contraction (dre04260), ECM-receptor interaction (dre04512), adrenergic signaling in cardiomyocytes (dre04261), calcium signaling (dre04020), regulation of actin cytoskeleton (dre04810), and cytokine-cytokine receptor interactions (dre04060) are a few indicators of such alterations.
The top downregulated gene in PTZF, mir1-1, was implicated in the regulation of cardiac contraction and embryonic angiogenesis in ZF [48]. In ZF, the downregulation of mir-1 is associated with the dysregulation of muscle gene expression, thereby affecting sarcomeric actin assembly. Knocking out mir-1 led to the failure of embryonic blood vessel formation. In fact, when observed under the microscope, PTZF showed malformed blood vessels, resulting in abnormal blood flow (Table 1). RNA sequencing analysis also showed that the expression of tnnt2, ttn, heg1, mybpc2, trdn, and itga10 were downregulated in PTZF compared with NTZF (Table S2). Silent heart mutants lacking the cardiac troponin T gene (tnnt2) failed to form myofibrils [49], and in the absence of titin, nascent myofibrils developed normally, but higher-order sarcomeric structures were absent [50]. Moreover, heart development protein with EGF-like domains 1 (heg1) mutants exhibited atrial and ventricular enlargement, reduced heart rate and blood flow, and PE [51]. Myosin binding protein C (encoded by mybpc) is critical for the maintenance of sarcomere integrity; the knockdown of mybpc caused cardiac hypertrophy, diastolic heart failure, and PE [52]. In mice, Trdn knockout led to impaired heart muscle excitation-contraction coupling [53]. Integrins (itga10) are cell surface receptors that play a significant role in mechanotransduction, and therefore the heart rate [54]. In comparison to NTZF, PTZF exhibited a 10-fold overexpression of cyp1a (Table S2). In the ZF heart failure model, Cyp1 suppression reduced doxorubicin-induced cardiomyopathy [55]. These findings, consistent with deformity data such as pericardial edema symptoms and a decreased heart rate ( Table 1 and Figure 2), suggest that phosmet can interfere with cardiogenesis and the cardiac conduction system.
The cytoplasmic Ca 2+ in the heart area increases several folds during cardiogenesis [56], and therefore, perturbation of calcium signaling interferes with normal cardiogenesis and subsequently, cardiac function. For instance, the downregulation of pkd2, as observed in PTZF (Table S2), encodes polycystin-2, a non-selective calcium-regulated cation channel expressed on the endoplasmic/sarcoplasmic membrane, resulting in reduced cardiac contractile function [57]. There is evidence that mitochondrial dysfunction is linked to cardiovascular manifestations [58] and indeed, mitochondrial inhibitors induced PE, bradycardia, and arrhythmia in the ZF model [59]. Thus, the cardiovascular abnormalities of PTZF can be attributed in part to mitochondrial dysfunction and calcium signaling perturbations during development.
According to our KEGG analysis, the aminoacyl-tRNA biosynthesis pathway was highly enriched in PTZF (Figure 6), indicating that protein biosynthesis was disrupted. Most genes encoding aminoacyl tRNA synthetase (ARS), which ligate tRNAs to their corresponding amino acids, were downregulated (Table 2). Thus, the growth retardation observed in PTZF (Figure 3) is attributable to abnormal protein synthesis. In addition to their important involvement in protein biosynthesis, ARSs also required in other biological processes. Recently, Inoue et al. reported that leucyl-tRNA synthetase deficiency induces autophagy in ZF [65]. Loss of leucyl-tRNA synthetases leads to liver failure in ZF [66]. Waldron et al. reported that histidyl-tRNA synthetases are essential for the proliferation and survival of neuronal progenitors during development [67].
Considering that ZF shares a high degree of gene homology with humans [15], understanding the mechanism of pesticide toxicity in ZF can provide insights into its consequences in humans. The ZF model can produce in vivo toxicity assessments within weeks, which is much faster than that for mammalian testing [73]. Animal toxicology studies frequently uncover effects that demand further investigation to elucidate the causes, which is both expensive and time consuming. The existing screening methods in ZF enable understanding very early details of off-target effects, and when combined with transcriptome analysis, understanding the potential mechanism and crosstalk between mechanisms. In addition, the impact of pesticides on the cardiovascular system, central nervous system, digestive tract, and visual functions and their proconvulsant potential can be assessed at the molecular level. Therefore, ZF-based technology should be recognized as an important prefilter to aid in the identification of the safest candidates for pesticide development and the determination of pesticide poisoning treatments. This study uncovered various molecular mechanisms affected by pesticide toxicity, thereby demonstrating the feasibility of using novel techniques for cumulative risk assessment.

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

ZF Maintenance and Embryo Collection
ZF (AB strain) were maintained in a glass aquarium filled with dechlorinated tap water and equipped with a recirculating filtration system. The photoperiod was 14 h light/10 h dark and the temperature was 25 ± 1 • C. ZF were fed 2-3 times a day with live brine shrimp (INVE aquaculture, Dendermonde, Belgium), blood worms (Hikari Bio-pure, Himeji, Japan), and dry flake food (Top meal, Tabia, Korea). The protocols used for ZF embryo collection and cleaning were as described [70,74]. Briefly, five females and five males were transferred to a nylon mesh, which was installed in a plastic aquarium filled with water, and covered with a transparent lid. On the next day, the fertilized eggs were collected 30 min after the light was turned on and washed three times in E3 medium (290 mg of sodium chloride, 8.3 mg of potassium chloride, 48 mg of calcium chloride, 81.5 mg of magnesium chloride, and 100 µL of 1% methylene blue per 1 L of distilled water).

Pesticide Treatment and Deformities Scoring
The final test concentration of phosmet [2-dimethoxyphosphinothioylthiomethyl) isoindoline 1,3-dione] (#36195, purity 98%, Sigma Aldrich, St. Louis, MO, USA) was prepared in E3 buffer. The test dose of 8.0 mg/L was chosen based on our prior work [19] as having twice the EC 50 value (4.38 ± 0.18 mg/L) of phosmet to ZF, ensuring that most scored abnormalities are present in all collected embryos (RNA sequencing samples). Phosmet and E3 medium exposures were started at approximately 2.0 hpf in an incubator with a constant temperature of 26 ± 1 • C and under dark conditions. The test solutions were replaced every 24 h. The experiments were conducted in 24-well plates and repeated twice with two 24-well plates (each plate containing 20 embryos) per test condition. Deformities were scored under a stereoscopic microscope (Stemi 508, Zeiss, Germany). The methods of deformities scoring and percentage calculation were as described [18,74]. Abnormalities in the somites and tail detachment, as well as symptoms of edema, were scored at 24 hpf. Accumulation of retina pigmentation, abnormal blood flow (caudal region), and hyperemia were scored at 48 hpf. At 72 hpf, hatching, pericardial, and yolk sac edema were scored. Body curvature and hatching were scored at 96 hpf. The deformity percentage was calculated as the ratio of embryos that showed deformity overall live embryos at the time point. Statistical significance was evaluated using an unpaired t-test.

Heartbeat Survey
A heartbeat survey was conducted at 48, 72, and 96 hpf. Embryos or larvae were anesthetized with a final concentration of 10 mg/L tricaine (ethyl 3-aminobenzoate methanesulfonate) (Sigma-Aldrich, St. Louis, MO, USA) just before being observed. The heartbeats of embryos or larvae were counted for 20 s under a stereoscopic microscope at room temperature of 26 ± 1 • C. The counts were adjusted to bpm and presented. A total of 10 embryos from each test condition at each time point were analyzed, and the experiment was performed 3 times. Statistical significance was evaluated using an unpaired t-test.

TEER and VMER
TEER was monitored at 96 hpf by gently touching the head and tail of ZF with a flexible plastic wire (1 mm). VMER was monitored at 96 hpf by gently moving a 1 mm diameter plastic tip in the field of binocular vision. The experiments were repeated three times, each time with ten embryos per test condition. Statistical significance was evaluated using an unpaired t-test.

Body Length Survey
A body length survey was conducted at 144 hpf. Pictures were captured with a stereomicroscope. Body length (from the mouth tip to the end of the tail fin) was obtained using OptiView 3.7 software (Korealabtech, Seongnam, Korea). The experiments were repeated twice with 12 embryos per test condition. Statistical significance was evaluated using an unpaired t-test.

RNA Sequencing Preparations
Three phosmet-treated and three non-treated ZF samples were collected by snap freezing in liquid nitrogen and were stored at −80 • C until use. Each sample contained approximately 60 embryos or larvae. Total RNA was extracted using the QIAzol®Lysis Reagent (QIAGEN, Cat. No. 79306, Hilden, Germany) and RNeasy®Mini Kit (QIAGEN, Cat. No. 74106, Hilden, Germany). Total RNA concentration was quantified using Quant-IT RiboGreen (Invitrogen, #R11490, Waltham, MA, USA), and RNA integrity was assessed using the TapeStation RNA screen tape (Agilent, Santa Clara, CA, USA). RNA preparations with a RIN greater than 7.0 were used to construct the RNA libraries. Illumina TruSeq Stranded Total RNA Library Prep Gold Kit (Illumina, Inc., San Diego, CA, USA, #20020599) was used to prepare the independent libraries using 0.5 µg of total RNA from each sample. The Ribo-Zero rRNA Removal Kit (Human/Mouse/Rat Gold) (Illumina, Inc., San Diego, CA, USA) was used to remove rRNA from total RNA. Then, the remaining mRNA was fragmented into small pieces using divalent cations under high temperatures. Using SuperScript II reverse transcriptase (Invitrogen, #18064014, Waltham, MA, USA) and random primers, the cleaved RNA fragments were transcribed into the first strand of cDNA. Second-strand cDNA synthesis was then performed using DNA polymerase I, RNase H, and dUTP. These cDNA fragments were then subjected to end repair, followed by the addition of a single A base and adapter ligation. After purifying the products, the final cDNA libraries were prepared by PCR. The KAPA Library Quantification Kits for Illumina Sequencing Platforms (KAPA BIOSYSTEMS, #KK4854, Bath, UK) were used to quantify the libraries, and TapeStation D1000 ScreenTape (Agilent Technologies, # 5067-5582, Santa Clara, CA, United States) was used to qualify libraries. The indexed libraries were subjected to paired-end (2 × 100 bp) sequencing by Macrogen Incorporated on an Illumina NovaSeq platform (Illumina, Inc., San Diego, CA, USA).

Transcriptome Alignment and Analysis of DEGs
DEGs were selected (nbinomWaldTest using DESeq2) with a cut-off of absolute log2 fold change ≥ 2 and a p-value of 0.05 by comparing PTZF and NTZF. Further analyses were performed by using the obtained DEGs.

Quantitative Real-Time PCR
Total RNA was extracted from three phosmet-treated and three non-treated ZF samples using the RNeasy Mini Kit (Bioneer, Daejon, South Korea), and analyzed with the Nanodrop 2000 (Thermo Fisher Scientific, MA, USA). RNA was reverse transcribed into cDNA using the ReverTra Ace™ qPCR RT master mix with gDNA remover (Toyobo, Osaka, Japan). Quantitative real-time PCR was performed on the CFX96 Dx real-time PCR detection system (Biorad, CA, USA) using the TOPreal™ qPCR 2X premix (Enzynomics, Daejon, Korea). Each sample was analyzed in triplicate and the average Ct values were used for fold change calculations. β-actin was used as a housekeeping gene. The 2 −∆∆CT method was applied for the expression analysis. The list of primers used for qPCR is presented in Supplementary Table S1.

Conclusions
In conclusion, these results demonstrate that phosmet alters the expression of a plethora of genes, impacting fish overall growth and behavior. These findings primarily provided detailed information on the altered molecular pathways to phosmet that causes developmental toxicity in ZF.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/genes12111738/s1, Figure S1: List of enriched pathways after phosmet treatment, Table S1: List of primers used for, Table S2: qPCRList of top upregulated and downregulated genes after phosmet treatment. File S1: List of DEGs after phosmet treatment.