mRNA Transcript Analysis of Hypothesis-driven Pathways as Known Responders to Organophosphate Exposure: Rhinella Arenarum Larvae Transcriptome Study


 Transcriptional analysis of the network of transcription regulators and target pathways in exposed organisms may be a hard task when their genome remains unknown. We used a whole transcriptome study on Rhinella arenarum larvae exposed to the organophosphorus pesticides azinphos-methyl and chlorpyrifos to evaluate transcriptional effects on a priori selected groups of genes. This approach allowed us to evaluate the effects on hypothesis-selected pathways such as target esterases, detoxifying enzymes, polyamine metabolism and signaling and regulatory pathways modulating them. We could then compare the responses at the transcriptional level with previously described effects at the enzymatic or metabolic levels to obtain global insight into toxicity-response mechanisms. The effects of both pesticides on the transcript levels of these pathways could be considered moderate, while the responses elicited by chlorpyrifos were more potent and earlier than those elicited by azinphos-methyl. Finally, we infer a prevailing downregulation effect of pesticides on signaling pathways and transcription factor transcripts encoding products that modulate/control the polyamine and antioxidant response pathways. We additionally tested and selected potential housekeeping genes based on those reported for other species. These results allow us to go through future confirmatory studies on pesticide gene expression modulation in toad larvae.


Introduction
The amphibian Rhinella arenarum (Hensel 1867) is widely distributed throughout Argentina and partially in South America 1 . Its life cycle includes two fundamental stages before reaching the adult stage: the embryonic and larval stages. The physiological reproduction of the species in the Upper Valley of Rio Negro and Neuquén (North Patagonia) occurs once a year during the spring months in backwaters and irrigation channels. The reproductive season coincides with the period of greatest pesticide application to protect fruit production (INTA, 1993). Various kinds of pesticides have been detected in super cial waters where amphibian reproduction occurs 2 . This situation implies that both R. arenarum embryos and larvae are potentially exposed, at least temporarily, to high concentrations of these toxicants, posing a hazard to this and other species that inhabit the area 3,4,5,6,7,8,9,10 .
When evaluating the risks of exposure to a contaminant for any species, it is desirable to count earlyresponse biomarkers capable of anticipating irreversible damage that may occur at later stages of development 11 . Molecular targets, effectors or modulators of toxicant effects are among early-response biomarkers. However, their development in autochthonous, nonmodel species such as R. arenarum is di cult due to the lack of both sequenced genomes that allow primer design and transcript analysis, as well as speci c antibodies necessary to support protein expression analysis 3 . In these cases, toxicogenomic and particularly transcriptomic analyses through bulk RNA sequencing (RNA-Seq) are powerful tools that enable screening for effects and responses to a toxicant and therefore allow further molecular studies 12,13 .
RNA-Seq is a highly sensitive and accurate tool for measuring expression across the transcriptome, enabling researchers to detect changes that otherwise would go unnoticed. It enables the detection of both known and new features in a single assay, uncovering isoforms, fused genes, single-nucleotide polymorphisms (SNPs), and other functional features without the limitation of prior knowledge.
Transcriptome sequencing makes it possible to obtain a genic expression pro le under a given stimulus at a given moment of time. It also makes it possible to quantify the levels of abundance or relative changes for each transcript during a speci c developmental stage or under a speci c treatment 14 . In the case of nonmodel organisms such as R. arenarum, RNA-Seq technology allows the obtention of genomic information through a relatively accessible process from an economic point of view and considering the cost-bene t relationship 13 .
In the context of RNA-Seq and transcriptomic analysis performed on R. arenarum larvae exposed to the organophosphorus (OP) pesticides azinphos-methyl (AZM) and chlorpyrifos (CPF) 13 , we performed hypothesis-driven data analysis processing on a priori selected mRNA transcripts. The aims of the study were a) to compare transcript expression levels in genes from pathways that have been previously recognized as impacted by OP pesticides and that have also been studied at the biochemical, metabolite or physiological level in R. arenarum to obtain a more comprehensive mechanism of toxicity and response and b) to test a set of potentially adequate housekeeping genes in R. arenarum larval stages for future quantitative PCR studies.

Chemicals
High purity-certi ed standards of azinphos-methyl (98.3% AZM) and chlorpyrifos (CPF; 99.5% purity) were purchased from Chem Service Inc. (West Chester, PA, USA). Standard solutions of 18 g L -1 AZM and 1 g L -1 CPF were prepared by dissolving the pesticide standards in acetone. The exact concentrations of AZM and CPF in the standard solution were checked by capillary gas chromatography coupled to a nitrogenphosphorus detector (GC-NPD).

Biological material
Adult females and males of the South American common toad (Rhinella arenarum) were collected in reference areas, free of pesticide application, in agreement with the corresponding collection permission 040/2020 from the Environment Secretary of Río Negro Province, Argentina. Animals used in this study were maintained and treated with regard to the alleviation of suffering according to recommendations of the Guide for the Care and Use of Laboratory Animals (National Research Council 2011) 55 . The animals were kept in captivity outdoors for 24-48 h until their use. Female ovulation was induced by intraperitoneal injection of 2500 international units (IU) of human chorionic gonadotropin (ELEA Laboratory, Buenos Aires, Argentina.) and embryos were obtained by in vitro fertilization 5 . Embryos were maintained until they reached the complete operculum (CO) stage (stage 25, according to 56 ). Ten days after reaching the CO stage, the larvae were used for acute toxicity assays.

Acute toxicity assays
The whole protocol was approved by the Faculty Committee for Care and Use of Experimental Animals (CICUAL-Facultad de Ciencias Agrarias Universidad Nacional del Comahue 01/13-7 -2020). Larvae were randomly collected to perform the assays. Sublethal concentrations of AZM (0.5 mg L -1 , 1/20 96 h-LC50; 4 ) and CPF (0.1 mg L -1 ), 1/15 96 h-LC50; 8 ) were selected to carry out exposures for up to 24 h in glass dishes, maintaining a ratio of 1 larva/10 mL in amphibian Ringer's solution with 0.3% acetone ( nal v/v). The exact pesticide concentrations were checked by gas chromatography and nitrogen-phosphorus detection. Control acetone treatment was included to discard possible solvent effects. The treatments were carried out in duplicate, and R. arenarum larvae were grown in 10 different glass receptacles to perform AZM/CPF 6 h exposures, AZM/CPF 24 h exposures, and control treatments. From each receptacle, fteen random larvae were collected and pooled at the corresponding times. Larvae were washed three times with cold Ringer's solution, placed in 1.5 mL tubes with RNALater® (Thermo Fisher Scienti c Inc.) and stored at -20 ° C until processed.
RNA extraction, cDNA library generation and massive parallel sequencing RNA extraction, cDNA library generation and massive parallel sequencing were carried out as described by Ceschin et al. (2020) 13 . Brie y, total RNA of each sample was extracted, and the cDNA library for transcriptome analysis was prepared. The ten library samples were normalized to 10 nM cDNA to be sequenced on a HiSeq 1500 Illumina platform, generating nonstrand speci c "paired-ends" (PE) 2 × 100 bp readings.

Statistical analysis of expression levels
The readings obtained by massive sequencing of the R. arenarum transcriptome were aligned with Bowtie2 v2.3.5 57 against TSA: GHCG00000000.1 (BioProject PRJNA485066), and transcript expression quanti cation was performed using RSEM v1.3.0 58, 59 . Expression values were normalized by the TMM method using the R and EdgeR packages 60,61 . Once the adequate HK transcripts were identi ed, the TMM values were standardized both by the average of the selected reference genes and by the respective control values. Finally, a nonparametric analysis was performed by the median and Kruskal-Wallis tests to assess signi cant differences or tendencies using the exact p-values. Raw data, standardization steps and statistical analyses are available in Supplementary le II.

Results
Gene selection from annotated transcripts for OP effect analysis.
We used a list of the available gene transcripts in R. arenarum published by us in the database from massive RNA transcript sequencing and gene annotation (http://rhinella.uncoma.edu.ar/; 13 ). We selected a group of 14 potential housekeeping genes, considering those currently proposed for mRNA expression normalization in Xenopus laevis 15,16,17 (Table S1, GROUP A). We found a considerable number of annotated genes that could be included in the following pathways, which were selected on the basis of available information about OP effects at the biochemical level in amphibians 6,8,18 and other vertebrates: polyamine metabolism genes, (Group B, 14 genes); antioxidant response genes (Group C, 16 genes); OPmetabolizing, primary-and secondary-target enzyme genes (Group D, 16 genes); and gene expression regulators, transcription factors and phosphorylation cascade effectors corresponding to the Mitogen-Activated protein kinase pathway, Transcription factor AP-1, Aryl hydrocarbon receptor (AHR) pathway and Nuclear factor erythroid 2-related factor 2 antioxidant response pathway (Group E, 17 genes) ( Figure   1). The complete list of genes is detailed in Supplementary Table S1.

Corroboration of gene annotation for selected transcripts
Alignments were carried out using the tools available in BLAST (BLASTN, BLASTX and BLASTP) to verify that the sequences of the selected transcripts effectively corresponded to the annotated genes. Those sequences that yielded a match greater than 50% in the vertebrate databases were selected for further analysis. A summary of this sequence validation is shown in Figure 1, and the complete analysis is available in Supplementary File I. From a total of 225 transcript sequences that were compared to the corresponding annotated genes, 93.8% could be effectively con rmed for further analysis; the best correspondence was obtained in group E covering transcription factors and signaling pathways with a 100%, and the lowest percentage was approximately 90% for the group of antioxidant stress transcripts (Group C).

Validation of adequate transcription levels for ltered transcripts
Previously validated transcripts were analyzed to determine if the transcription levels of treatment replicates were appropriate for statistical analysis. These levels resulted from the massive RNAseq ampli cation and quantitation of each transcript fragment in the different treatments. As 'not appropriate' levels, we considered fragments with zero expression levels in any sample and/or replicates with very low and erratic values. The summary of transcripts accepted in each group of selected genes is shown in Figure 1. The raw data corresponding to all transcript fragments are available in Supplementary File II.

Housekeeping genes stability
We carried out a transcription stability analysis of the selected potential housekeeping genes, comparing their variability within and between the treatments. The respective TMM means, minimum and maximum expression values, standard deviations, median expressions and percentual coe cients of variation (CV%) are presented in Table 1. We assumed a maximum CV% of 20% in the TMM as an acceptable limit for a transcript to be considered housekeeping. In this way, we were able to select 9 out of 12 transcripts as appropriate reference genes in the rst larval stage of R. arenarum. The genes with the least and acceptable variations in group A were EF1A0, EF1GA, EF1D, EF1B, TBA, TBB, TBB4B, ACTB, and RL8.
The transcripts belonging to tubulin TBA1 and glyceraldehyde 3P dehydrogenase (G3P) showed a variability of approximately 22% in their expression levels, roughly in the exclusion limits we considered, and might be considered housekeeping genes if a re ned and speci c analysis proves better results. In turn, ACT3, belonging to sarcomeric actin, showed the greatest variation in expression at 83%, clearly indicating that it cannot be considered a housekeeping gene in R. arenarum larvae, at least for OP pesticide studies. The TMM values of the selected HKs, rated with respect to their control values, were averaged for each treatment and used to standardize the TMM values of the transcripts corresponding to the genes of groups B, C, D and E (raw and standardized data shown in Supplementary le II).
Effects of OP on the transcription of polyamine metabolism genes OP mainly decreased the expression of genes related to polyamine synthesis ( Figure 2A). Although some differences could be noted between CPF and AZM and the exposure times, decreases reaching 30-40% in the expression levels were found in ornithine decarboxylase (DCOR1), S-adenosylmethionine decarboxylase proenzyme (AMD, transcript 1-b), spermidine synthase (SPEE), and ornithine decarboxylase antizyme (OAZ1). Only AMD transcript 1-a showed relevant increases of nearly 3 times the control values when exposed to both OP pesticides. In turn, the ornithine decarboxylase regulator OAZ2 and the antizyme inhibitors AZIN1 and 2 showed no signi cant or no relevant changes (Supplementary Table S2).
In turn, polyamine-degrading enzymes showed more variable responses to OP pesticides ( Figure 2B). AOC1 transcript expression, corresponding to diamine oxidase, showed a signi cant inhibition of 60-75% by AZM treatment and 50-70% by CPF treatment. The AOC2 transcript showed a similar tendency with even deeper inhibition responses, but the difference was not signi cant (Table S2). AOC3 and 4 transcript isoforms showed variable and nonsigni cant responses to OP pesticide exposure. Acetyl spermidine/spermine oxidase (PAOX) expression showed an induction tendency of 40% with both OPs at 6 h of exposure, while spermine acetylase transcript expression (SAT1, 2-a, 2-b) showed an inhibitory trend with CPF up to 40%. Spermine oxidase expression was scarcely inhibited (20%) by AZM and induced by CPF (20%, 6 h) ( Figure 2B).

Effects of OP on oxidative stress response -antioxidant enzyme genes
Superoxide dismutase transcripts displayed differential responses to OP treatments ( Figure 3A). SODC expression showed a decrease in response to AZM exposure at 24 h and to CPF at 6 and 24 h of exposure (up to 30% inhibition) with respect to the control. SOD3 was less affected, as its transcription was signi cantly inhibited by approximately 25% by AZM at 24 h and CPF at 6 h of exposure. On the other hand, SODE expression was signi cantly induced in 60% by CPF at 24 h. Catalase transcription (CATA) was in turn slightly inhibited by AZM exposure at 24 h and by CPF up to 25%.
Glutathione-dependent peroxidase transcripts were mainly induced by OP pesticides. GPX1 transcription was considerably induced by AZM at 24 h (3X) and by CPF at both times (up to 2.7X). GPX3 transcription was induced to a lesser extent but in a signi cant way by CPF (1.5X), while GPX4 also showed an induction of approximately 30% by CPF exposure. In turn, the transcript corresponding to GPX8B presented a tendency to decrease in samples exposed to AZM at 24 and to CPF at 6 and 24 h, up to 60% with respect to the control. Two other GPX transcripts, GPX2 and 7, as well as one GSH-reductase transcript (GSHR), showed decreasing but not signi cant trends due to OP pesticide exposure (Supplementary Table S2). Finally, one transcript for the enzyme glutathione synthetase (GSHB) was scarcely downregulated by OP exposure (AZM at 24 h; CPF at both times, approximately 25%) ( Figure  3B).

OP pesticide targets and detoxifying enzymes
Although some fragments corresponding to cholinesterase transcripts were detected (BCHE), their levels were very low and showed erratic responses (raw data in Supplementary le II). Among the esterase group, carboxylesterase (EST5A) transcript showed a relevant decrease in larvae exposed to AZM at both 6 and 24 h, while in larvae exposed to CPF, a decrease at 6 h (50% inhibition) was observed, returning to control values at 24 h ( Figure 4A). A similar pattern was observed for another carboxylesterase transcript (EST3B), with no signi cant effects (Supplementary Table S2). The serum paraoxonase/arylesterase-2 transcript (PON2) showed a small but signi cant decrease of approximately 20% only in larvae exposed to CPF at 6 h with respect to the controls. For cytochrome P450 enzymes, only the transcripts putatively related to OP metabolization were analyzed; CYP1A1 transcript expression was signi cantly reduced by AZM at 24 h and by CPF exposure, although the effects were minor (20% inhibition). Cytochrome CYP2C19 expression was more markedly inhibited by both OP pesticides after 6 h of treatment, reaching 50-60% inhibition with respect to controls.
Glutathione-S-transferase (GST) transcripts were among the most abundant groups, as 15 different isoforms could be identi ed. The GST isoforms alpha (GSTA3), pi (GSTP1 and 2), and theta (GSTT3) were downregulated by both OP pesticides, with transcription inhibitions ranging from 20 to 50% with respect to controls ( Figure 4B). Isoform GSTP1 showed the highest inhibition and in a signi cant way. Although GST Mu (GSTM1) showed a decreasing trend, the effects were not signi cant (Table S2). On the other hand, GST Kappa 1 transcript (GSTK1) was signi cantly induced by CPF (45%), while GSTM3 showed the highest transcription induction values by CPF and by AZM at 24 h, reaching 2 X -3.3 X expression values with respect to controls, but with a p-value of 0.09. The GST omega 1 transcript (GSTO1) also showed a signi cant increase of 45% in its expression for AZM at 6 h and CPF at 24 h.
Microsomal GST transcripts were also analyzed, nding increases in the expression levels of isoforms 1, 2 and 3 (M-GST1, M-GST2, M-GST3, respectively, Figure 4C). The M-GST1 transcript was one of the most affected transcripts, with increases of 40% (AZM-24 h; CPF-6 h) to 70% (CPF-24 h). One M-GST3 transcript showed roughly similar effects (30-60% increase), while a second transcript showed a moderate increase (20%) only for CPF at 24 h of exposure. For the M-GST2 isoform, we also identi ed two transcripts, one of which showed a moderate increase (up to 35%) after CPF exposure, and the other showed slight decreases of approximately 20% in its expression with respect to the controls.

OP effects on transcription factor-and signaling pathway-related genes
The MAP kinase phosphorylation pathway seemed to be affected solely at the rst level of regulation, as mitogen-activated protein kinase kinases (MP2K) were downregulated early by CPF at 6 h and late by AZM at 24 h ( Figure 5); the MP2K2 transcript was diminished to 70% by both OP pesticides, while the MP2K1 transcript was affected only by CPF but not at a relevant level (Supplementary Table S2). Similarly, the mitogen-activated protein kinase p38 transcript was scarcely affected by both OP pesticides (approximately 12%), while the JUNK transcript remained unaffected. None of the transcription factor JUN transcripts detected (JUN-1, JUN-B, JUN-D1) were signi cantly affected by OP exposure (Table S2). In turn, the transcription factor FOS, which associates with JUN as a heterodimer partner in activator protein-1 (AP1), showed a severe reduction in its transcription levels by AZM at 6 h of exposure (to 35% of control levels) but recovered at 24 h and surpassed controls by 25%; CPF caused an inverse pattern, inducing transcription of FOS at 6 h (40%) but inhibiting it at 24 h (30%) ( Figure 5).
On the other hand, nuclear factor erythroid 2-related factor 2 (NRF2), linked to the antioxidant response through the antioxidant response element (ARE) pathway, showed a decrease in transcription of approximately 30% in the larvae exposed to AZM for 24 h and in those exposed to CPF for 6 and 24 h compared to controls ( Figure 5).
Within the aryl hydrocarbon receptor (AhR) pathway, some differences were found for the nuclear translocator (ARNT) and the aryl receptor repressor (AHRR). Two transcripts for ARNT were analyzed. The rst ARNT transcript showed a decrease of nearly 30% in larvae exposed to AZM at 24 h and to CPF at 6 and 24 h compared to the controls. The second transcript, ARNT2, showed a relevant induction in transcription of nearly 2X after 6 h of exposure to CPF. The AHRR transcript showed an increase in larvae after 24 h of exposure to AZM (80%) and CPF (40%) at 24 h with respect to the controls. The AhR transcript itself showed a decreasing but statistically nonsigni cant pattern after exposure to CPF (Table  S2). The transcript corresponding to HSP90AB1, the chaperone that binds AhR in the cytosol, showed a decrease of 20% after 24 h of exposure to CPF, while the transcript corresponding to the cochaperone aryl-hydrocarbon-interacting protein-like-1, AIPL1, was downregulated by AZM (approximately 50%) but increased by CPF (up to 75%). Another isoform, AIP, showed no effects, as well as the third cochaperone prostaglandin E synthase-3, PTGES3.

Discussion
We succeeded in applying a hypothesis-driven transcript analysis from a transcriptome database on R. arenarum larvae exposed to two OP pesticides, with two purposes into mind: 1) to test HK genes currently used in model species and 2) to compare the effects of AZM and CPF on transcript levels for pathways previously described as impacted at biochemical or protein expression levels. These goals in a native, nonmodel species such as the amphibian R. arenarum are a very good example of what transcriptomics analysis can do, solve or imply for the advance in molecular toxicology when other tools are not readily available 12,19,20 . Considering the results obtained, it is evident that the transcriptomic expression data allow a screening analysis into selected pathways to decide for further qPCR studies if necessary. Furthermore, we can con rm that transcriptomic information is an essential tool when developing molecular biology assays in nonmodel organisms.
One of the problems in gene transcription analysis is being able to accurately determine the ampli cation degree of the nucleic acid present in the reaction. This is commonly achieved using HK genes, that is, those genes constitutively expressed and consequently not affected by the exposure conditions. In a comprehensive study conducted with Xenopus laevis, RL8 and GAPDH were stable as reference genes during the rst embryonic stages, ODC1 was stable in the initial and nal stages, and H4 was adequate throughout embryonic development 15 . In another study in X. laevis, eEF1A1 and SUB1 L were identi ed as the genes whose expression remained more stable, which would allow their use as reference genes 16 .
In this work, we were able to analyze the expression levels of potential reference genes from the R. arenarum transcriptome data considering those proposed for X. laevis, selecting a list of 9 genes that would t to the different requirements in future qPCR studies. We con rm that RNA-Seq data have the potential to identify genes with less variation in their expression, as suggested by other authors 21,22,23 .
By analyzing the transcriptomic data of selected pathways of R. arenarum larvae exposed to the OP pesticides AZM and CPF, we obtained an overview of their regulatory effects on gene transcription at previously known targets. A qualitative representation of a reduced heatmap is shown in Figure 6. This method of representing transcription level data enables a rapid comparison among the different treatments or states 12 . The rst conclusion to be highlighted is that, in general, the exposure of R. arenarum larvae to AZM or CPF at sublethal concentrations and up to 24 h did not cause notable changes in gene expression in the selected pathways. This is very interesting, considering that target esterase genes, detoxifying genes and oxidative stress response genes might be a priori expected to be upregulated, as their protein products are either inactivated by OP pesticides or their activities are engaged in xenobiotic transformation and detoxi cation 24,25,26,27 . Most of the transcripts showing changes in their levels after OP exposure did so in a moderate way, and only a few showed fold-changes higher than 2X. Thus, moderate changes in mRNA expression in these genes would be enough for ample response changes in their product activities. Another observation is that exposure to AZM clearly provokes a relatively poor effect on the genes analyzed at the short exposure time of 6 h; AZM effects are mainly evident after 24 h of exposure. This may be seen in all the selected groups, except for polyamine metabolism genes, where most of their expression levels were altered. In turn, the CPF effects on R. arenarum were more potent, affecting most of the gene expression early at 6 h and sustaining the effects after 24 h of exposure. Another interesting feature to be taken into consideration is that downregulation cases notably exceed induction cases, in a proportion of 3 to 1.
The polyamine metabolism pathway was probably the most affected group of those analyzed in this work, as shown in Figure 6. The general pattern of the effects of both OP pesticides on polyamine synthesis suggests downregulation through successive involved genes. Ornithine decarboxylase, Sadenosyl methionine decarboxylase precursor and spermidine-spermine synthase expression were mainly inhibited, except for one AMD-1 transcript that showed an induction reaching 2-fold or more. Antizyme expression is also inhibited, but this may be a feedback effect due to ODC inhibition itself.
Complementary to these effects on polyamine synthesis, degradation genes tend to be downregulated by OP pesticides, probably as a cellular response intended to avoid a lethal drop in polyamine levels 28,29,30,31 . Nevertheless, some enzymes related to spermidine and spermine degradation show increasing expression in larvae exposed to CPF, which seems to elicit stronger effects than AZM. A decrease in putrescine and spermine was observed in R. arenarum embryos of the complete operculum stage exposed to AZM 32 . Similarly, CPF downregulated ODC activity and decreased putrescine and spermidine levels in early R. arenarum embryos, correlated with the percentage of embryonic developmental arrest 18 . The opposite effect was reported in R. arenarum late embryos exposed to AZM, where ODC activity and putrescine content were increased 31 . We also found concordance between the transcript expression levels and some of the polyamine-degrading enzymes and previously reported activity values in R. arenarum embryos exposed to AZM or CPF, i.e. the inhibition of DAO and SMOX activities and the increase in PAOX activity 32 . Therefore, there is a good correspondence between the expression levels found at the transcriptomic level for the regulatory enzymes, some of the polyaminemetabolizing enzyme activities and the polyamine contents in R. arenarum. On the other hand, spermidine did not show variation throughout embryonic development for AZM-exposed embryos 32 . These ndings are related to the ne regulation of polyamine metabolism, both at the transcriptional level and at the enzymatic level. According to studies performed in the rst embryonic stages in Xenopus, AMD overexpression leads to the interruption of normal embryonic development, causing cell apoptosis 33 .
Correspondingly, the increased expression of the AMD1 transcript elicited by exposure to AZM and CPF in our work may be related to an apoptotic mechanism triggered to eliminate abnormal and damaged cells caused by pesticide exposure. On the other hand, the increase in AMD1 transcription could also be due to the need to synthesize a greater amount of spermidine. Spermidine acts as a donor of the aminobutyl group for posttranslational modi cations of the lysine residue in the transcription factor eIF5A. This modi cation is essential for the activity of the transcription factor, whose function is to contribute to transcription, nucleus-cytoplasm transport, mRNA renewal and apoptosis 34 .
Oxidative stress has been reported in different organisms exposed to OP pesticides. In particular, the effects of AZM and CPF on oxidative stress and antioxidant responses at the metabolite and enzymatic levels have been studied in R. arenarum development. Indeed, we have proposed that there is a link between OP effects on polyamine metabolism and levels, oxidative stress and teratogenic effects in R.
arenarum embryos and larvae 5,6,7,8,18,24,31,32,35,36 . Smirnova et al. (2012) 37 analyzed oxidative stress as the cause of alterations in polyamine metabolism due to the dysregulation of ODC and SSAT; human hepatoma cells chemically induced to increase ROS production showed overexpression of ODC and SSAT, which are transcriptionally regulated by NRF2 through a speci c recognition site 38 . The downregulation of NRF2 mRNA expression with elevated Nrf2 protein levels has been reported in liver pathologies 39 . We found that OP exposure in R. arenarum larvae caused a downregulation of NRF2 mRNA expression, and the accompanying downregulation at the transcriptional level of ODC, SSAT and several antioxidant-detoxifying enzyme genes containing ARE sequences (such as GSTA and GSTP) suggests that Nrf2 protein or activity might also be downregulated. Ma et al. (2018) 40 refer to the alterations in transcription levels of detoxifying and oxidative stress-related enzymes in the GST and CYP groups in a transcriptomic analysis performed on trichlorfon-exposed Rana chensinensis. This OP pesticide upregulated CYP2C transcripts but downregulated CYP3A and GSTK1. We report an increase in GSTK1 expression because of CPF exposure in R. arenarum, as well as for other GST isoforms, such as mu and the microsomal isoforms. This upregulation in several GST transcripts is in concordance with the reported increase in GST activity using CDNB as a substrate in R. arenarum embryos and larvae after exposure to both AZM and CPF, among other OP pesticides 7,35,41 . The effect on GST activity was associated with an increase in GST-Pi1 protein in R. arenarum larvae exposed to arsenic 11 . Other genes under Nrf2 regulation are those belonging to antioxidant defenses, such as SOD3, glutathione peroxidases (GPx2, 3, 6 and 8), and glutathione reductases (GSHR1) 40,42 . Accordingly, we observed downregulated levels of SOD3 and GPX8 accompanied by NRF2 downregulation. Other SOD transcripts, glutathione synthase GSHB, and CATA, also seemed downregulated, while SODE, GPX1, 3 and 4 were induced, mainly by CPF. In our experience, the antioxidant enzymatic activity response varies greatly in R. arenarum embryos and larvae exposed to OP pesticides, showing cycles of induction followed by inhibition attributed to catalytic site inactivation due to ROS attack 8, 24,35,41,43 . In the response of R. arenarum larvae to arsenic-induced oxidative stress, we observed an increase in Cu/Zn SOD protein, followed by a decrease together with a decrease in CAT protein, with similar enzymatic activity patterns 11,44 .
We also report here the downregulation of CYP1A1 and CYP2A19 for both AZM and CPF exposures, probably as a negative regulation after a previous detoxi cation response. CYP1A1/2 are among the genes regulated by the AhR pathway by the binding of activated and nuclear-translocated AhR-Arnt heterodimeric transcription factor to DRE/XREs in their regulatory sequences. The AhR protein is currently stabilized in the cytosol by the binding of the heat-shock protein HP90AB1 as a chaperone, together with the X-associated protein AIP and the prostaglandin E3 synthase PTGES3; once activated by planar aromatic hydrocarbon binding, it translocates to the nucleus, where it binds to the nuclear translocator Arnt. There is also an AhR repression mechanism involving protein-protein interactions with AHRR 45,46 . Studies carried out with Sparus aurata exposed to PCBs have shown a different expression pattern between AHRR and AhR-ARNT, in accordance with those studies that suggest that multiple mechanisms can contribute to the downregulation of AhR 47 . Our ndings coincide with these studies, as we observed induced levels of AHRR, while the levels of ARNT were decreased in R. arenarum larvae exposed to both CPF and AZM. Furthermore, the decreased levels of AhR-ARNT are consistent with the decreased expression found in CYP1A, since the entire mechanism is downregulated. A decrease in HSP90 and AIP transcripts, supposing reduced chaperone levels, would also contribute to lower AhR protein levels 48 . The AhR pathway is also involved in paraoxonase-1 expression regulation, according to results on exposure to polyphenols and speci c inducing ligands 49 . In turn, the PON-2 gene has at least a CRE regulatory sequence for AP-1 regulation linked to the oxidative response and JNK activation and a PAPR-regulated site related to polyphenol activation and MEK pathway downregulation by phosphorylation 50,51 . Our results agree with downregulated MEK (MP2Ks), AP-1 (cFOS) and AHR pathways causing PON2 transcript repression in R. arenarum larvae exposed mainly to CPF. Similarly, the xenobiotic/endobiotic detoxifying carboxylesterase family (CES) is transcriptionally regulated by a series of nuclear factors and pathways: AhR, constitutive androstane receptor (CAR), pregnane X receptor (PXR) and Nrf2 are involved mostly in the upregulation of some of these families 52 . Although we were not able to annotate putative transcripts corresponding to CAR and PXR pathways in our R. arenarum transcriptome, the downregulation trends followed by AhR and Nrf2 pathways in larvae exposed to AZM and CPF are in line with CES5A transcript downregulation. This may be surprising, in the sense that carboxylesterases are recognized suicide "buffer" enzymes that irreversibly react with OP pesticides to protect the primary target acetylcholinesterase in the nervous system, decreasing their activities. This is in fact corroborated in toad embryos and larvae exposed to different OP pesticides, including AZM and CPF, in most of our reports in R. arenarum 3,7,8,35,53 . Thus, an induction of carboxylesterase and cholinesterase mRNA and/or protein synthesis would be expected to reestablish normal activity levels. Finally, remarking on the complexity of the responses and crosstalk between pathways, we have reported for different developmental states in R. arenarum and cell cultures exposed to CPF, other OP pesticides or arsenic, an increase in MEK1/2 and ERK1/2 proteins, their translocation to the nucleus, and ERK phosphorylation, considering that the MAP kinase pathway regulates the Nrf2-mediated response. Additionally, cFOS and cJUN proteins are increased after oxidative stress in R. arenarum embryos and larvae 43,44,54 .

Concluding remarks
In conclusion, we acknowledge the power of a hypothesis-driven transcriptomic data analysis on selected pathways. First, we identi ed an appropriate battery of potential housekeeping genes in R. arenarum larvae to further analyze gene expression by conventional or quantitative PCR after pesticide exposure. Second, we were able to visualize the effects of two OP pesticides in a priori selected metabolic and signaling pathways and compare them to previous metabolite and enzyme activity analyses. In our analysis, we infer a gradient of effects since CPF is more potent than AZM and acts earlier on gene transcription. Finally, we found a prevailing downregulation in signaling cascades and transcription factors that act upstream of the polyamine metabolism pathway and antioxidant responses, which partially coincides with previously well-characterized responses at the protein and/or activity and metabolite levels.

Declarations
Gene ID Data calculated from massive RNA-seq and transcriptomic analysis are expressed as TMM values. A limit of 20% for the CV is proposed to consider a transcript as an HK gene in expression analysis. Figure 1 Gene selection process by metabolic pathway.

Figures
Hypothesis-driven selection into pathways recognized as affected by organophosphorus pesticides at enzyme activity or metabolic product levels. HK: housekeeping genes; PM: polyamine metabolism; AS: antioxidant system; DS: detoxifying systems; TF: transcription factors, signaling pathways. Once selected in a rst instance, transcripts were checked by nucleotide and predicted amino acid sequence identity; afterwards, the transcription levels were checked to discard erratic low values.

Figure 2
Expression levels of transcripts corresponding to polyamine metabolism genes in R. arenarum larvae. A, Genes related to polyamine synthesis and its regulation. B, genes related to polyamine degradation.
Signi cance levels for Kruskal-Wallis and median tests, * p=0.09; ¤ p=0.08; † p=0.04; III p=0.0001. AZM: azinphosmethyl, CPF: chlorpyrifos, at 6 h and 24 h exposures; gene codes are detailed in the text.   Expression levels of transcripts corresponding to transcription factors and signaling pathway genes in R. arenarum larvae. Signi cance levels for Kruskal-Wallis and median tests, * p=0.09; ¤ p=0.08. AZM: azinphosmethyl, CPF: chlorpyrifos, at 6 h and 24 h exposures; gene codes are detailed in the text. Heatmap representation of OP effects on selected gene expression in R. arenarum larvae.
Data are presented for four groups of selected genes whose products or activities are known or suspected targets of OP pesticides. Larvae were exposed to azinphos-methyl (AZM) and chlorpyrifos (CPF) for 6 and 24 h. Gene codes are detailed in the text.