- Email: [email protected]
Transcriptome analysis of Rana chensinensis liver under trichlorfon stress
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Contents lists available at ScienceDirect
Ecotoxicology and Environmental Safety journal homepage: www.elsevier.com/locate/ecoenv
Transcriptome analysis of Rana chensinensis liver under trichlorfon stress Yu Ma a b c
a,b
a,b
, Bo Li
b
c
, Yang Ke , Yongan Zhang , Yuhui Zhang
a,⁎
MARK
College of Life Science, Shaanxi Normal University, Xi’an 710062, China Shaanxi Microbiology Institute, Xi’an 710043, China Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan 430072, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Organophosphate insecticides Trichlorfon Differentially expressed genes Rana chensinensis RNA-seq
Trichlorfon is a selective organophosphate insecticide that is widely applied in aquaculture and agriculture for control of various parasites. However, repeated and excess applications of trichlorfon often lead to water pollution and threaten non-targeted species. Our previous studies showed that trichlorfon could cause oxidative stress, lipid peroxidation and hepatic lesions in the liver of Rana chensinensis, but the related molecular mechanisms remain unclear. To explore the interference of trichlorfon in gene transcription, the differentially expressed genes in the liver of R. chensinensis exposed to trichlorfon were characterized using the RNA-seq platform. A search of all unigenes against non-redundant protein sequence (Nr), non-redundant nucleotide (Nt), Swiss-Prot, Kyoto Encyclopaedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG) and Gene Ontology (GO) databases resulted in 22,888, 21,719, 20,934, 16,923, 7375 and 15,631 annotations, respectively, and provided a total of 27,781 annotated unigenes. Among the annotated unigenes, 16,923 were mapped to 257 signalling pathways. A set of 3329 differentially expressed unigenes was identified by comparison of the two groups in liver. Notably, relative expression of metabolism-related genes, including both upand down-regulated genes, were also validated by qPCR. The present study depicts the high degree of transcriptional complexity in R. chensinensis under trichlorfon stress and provides new insights into the molecular mechanisms of organophosphate insecticide toxicology. Some of these metabolism-responsive genes could be useful for understanding the toxicological mechanism of trichlorfon on non-target aquatic organisms and will contribute to the conservation of aquatic life.
1. Introduction Organophosphates have been used as insecticide, acaricide and anthelmintic compounds for several decades. Trichlorfon is a selective organophosphate insecticide used for control of various parasitic infestations in fish and agriculture. The dosage of trichlorfon applied to aquaculture ponds to eradicate ectoparasites varies from 0.1 to 1.0 mg/ L (Herwig, 1979; WHO, 1992). Trichlorfon has been widely used in aquaculture, however, the trichlorfon-polluted waters from the fish and agriculture industries are potential threats to non-targeted species, such as fish, crabs, and shrimps (Soumis et al., 2003; Chang et al., 2013). Therefore, the expansion of aquaculture worldwide has generated concerns about its impact on both public and environmental health (Coelho et al., 2011). Because the effects of environmental concentrations very rarely cause lethality, the risk assessment of trichlorfon as a long-term pollutant in aquatic ecosystems remains to be comprehensively investigated in a specific species. The liver plays a critical role in the detoxification of xenobiotics.
⁎
Biochemical and histopathological data have shown the hepatotoxic tendency of trichlorfon in aquatic organisms (Mataqueiro et al., 2009; Xu et al., 2012). The metabolic changes observed in aquatic animals exposed to trichlorfon can create widespread disturbances in enzymatic activities and oxidative metabolism in general physiological processes (Thomaz et al., 2009; Coelho et al., 2011; Xu et al., 2012; Chang et al., 2013). These studies suggest that the metabolism of xenobiotics is a key event in trichlorfon-exposed organisms. However, the molecular mechanism of the metabolic pathways involved in xenobiotic detoxification induced by trichlorfon is still unclear. Previous studies have found that superoxide dismutase (SOD) gene expression levels were up-regulated after exposure to 0.4 mg/L trichlorfon for 6 h in the giant freshwater prawn (Macrobrachium rosenbergii) (Chang et al., 2013). Indeed, many more genes may participate in this process and need to be identified in more species. The development of techniques such as RNA-seq have provided a better understanding of the molecular mechanisms of toxicological responses in non-model organisms. RNA-seq is a new method for both
Corresponding author. E-mail address: [email protected] (Y. Zhang).
http://dx.doi.org/10.1016/j.ecoenv.2017.09.016 Received 8 March 2017; Received in revised form 5 September 2017; Accepted 7 September 2017 0147-6513/ © 2017 Elsevier Inc. All rights reserved.
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
exposure. Then, the tissue samples were immediately collected in cryotubes and stored in liquid nitrogen for later RNA extraction and qPCR. All procedures involving the handling and treatment of frogs used in this study were conducted in accordance with the approval of the Animal Care and Use Committee of Shaanxi Normal University.
mapping and quantifying the transcriptome; it has clear advantages over the existing approaches and improves our understanding of the extent and complexity of the eukaryotic transcriptome including nonmodel organisms that lack completely sequenced genomes (Wang et al., 2009a, 2009b; Ekblom and Galindo, 2011). RNA-seq has been recently used in some amphibian species, including Odorrana margaretae (Qiao et al., 2013) and Rana chensinensis (Zhang et al., 2013), but the use of the high-throughput RNA-seq technology for studying environmental stressors is still limited. This study aims to understand the molecular mechanism of trichlorfon impact on liver function in amphibians. Amphibians spend the first phase of the life cycle in aquatic environments and their highly permeable skin makes them particularly sensitive to contaminants. Furthermore, their low mobility and high diversity in their reproductive modes through critical hormone-related developmental stages in aquatic environments render them vulnerable to environmental changes. The Chinese brown frog (Rana chensinensis) is a unique amphibian species that is widely distributed in northern China. This species is considered a potential candidate indicator due to its high sensitivity to chemical pollution in an aquatic environment; thus, it can be conveniently used as a native species to research ecotoxicology in local amphibians (Lundberg et al., 2007; Li et al., 2014, 2016). Our previous studies showed that trichlorfon can cause oxidative stress, lipid peroxidation and hepatic lesions in the liver of R. chensinensis (Li et al., 2017). These findings suggest that the liver is the major organ targeted by trichlorfon. However, the related molecular mechanisms remain unclear. The present study aims to obtain transcriptome information for R. chensinensis exposed to trichlorfon by using newly developed RNA-seq technology. The mRNA expression levels of four genes involved in metabolism and detoxification (eNOS, Cyt 2c, Gst kappa I,and Cyt 3a) were examined using quantitative real-time PCR (qPCR). The resulting data provides a preliminary research direction for the molecular mechanisms underlying environmental exposure, enriches the information currently in the cDNA database and identifies the genes involved in the metabolism and detoxification of trichlorfon toxicity-related mechanisms in amphibians.
2.2. RNA isolation and deep sequencing Total RNA was isolated from liver tissues from trichlorfon-exposed and unexposed frogs at 4 weeks using a TRIzol Kit (Invitrogen, CA, USA) according to the manufacturer's instructions. To obtain as many gene transcripts as possible and avoid the negative effects of individual variation, equal volumes of high-quality RNA from 12 individuals from each group were then pooled together for transcriptome sequencing. The RNA integrity was checked using an Agilent 2100 BioAnalyzer (Agilent Technologies, CA, USA), and the mRNA was purified from the total RNA using the PolyATract mRNA Isolation Systems Kit (Promega) following the manufacturer's instructions. Fragmentation was carried using divalent cations under elevated temperature in NEB Next First Strand Synthesis Reaction Buffer. First stand cDNA was synthesized using random hexamer primer and M-MuLV reverse transcriptase (RNase H-), and then the second strand cDNA synthesis was performed using DNA polymerase I and RNase H. In order to select cDNA fragments that were preferentially 150–200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). The adaptor modified fragments were selected by gel purification and amplified through PCR to create the final cDNA library. Transcriptome sequencing was conducted using an Illumina Hiseq™ 2000 according to the manufacturer's instructions. The invalid reads which contain adapters or unknown or low-quality bases, were abandoned to avoid negatively affecting the bioinformatics analysis. De novo transcriptome assembly was carried out using the short-read assembly program Trinity software (Grabherr et al., 2011). 2.3. Bioinformatic analysis
2. Materials and methods
Unigenes from the de novo assembly were further processed by sequence annotation using the NCBI (http://blast.ncbi.nlm.nih.gov/) non-redundant protein sequence (Nr) and non-redundant nucleotide (Nt) (Pruitt et al., 2007), Swiss-Prot (http://www.uniprot.org/) (Boeckmann et al., 2003), Clusters of Orthologous Groups (COG) (http://www.ncbi.nlm.nih.gov/COG/) (Koonin et al., 2004), Gene Ontology database (GO) (http://www.geneontology.org/) (Götz et al., 2008), and Kyoto Encyclopaedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/pathway.html) databases (Kanehisa and Goto, 2000). Gene names were assigned to each protein sequence based on the best Blast hit from all the Blast results.
2.1. Animals and treatments Adult Chinese brown frogs, Rana chensinensis, were obtained from the rivulet in unpolluted areas in Qinling Mountains (Shaanxi, China; N34°0′30.95′′, E109°06′48.64′′). The frogs were acclimatized in aquariums for 7 days prior to the experiment. After that, a total of 48 healthy frogs of both sexes, weighing approximately 8.967 ± 1.383 g, were randomly sampled and divided into experimental and control groups. Based on the preliminary experiment and the applied dose of trichlorfon in aquaculture, the exposure trichlorfon dose at 0.1 mg/L (90% purity, Hubei, China) was selected, and a blank control group was set. Each aquarium (40 cm × 20 cm × 20 cm) had 8 individuals with 2 L of test solution for static-renewal assay, and there were three replicates for each treatment. During the acclimation and experimental periods, the frogs were fed worms (Tenebrio molitor) three times a week. The tanks were placed on a slant to provide the option of both aqueous and dry environments. Trichlorfon-treated water was fresh prepared and renewed every 48 h to keep the trichlorfon concentration constant, and the tank was cleaned thoroughly. Throughout the experiment, the water temperature was maintained at 20 ± 2 °C, the pH values ranged from 6.9 to 7.1, the dissolved oxygen concentration remained at 7.0 ± 0.1 mg/L, and the photoperiod was set as 12:12 light:dark. GDYS-201 M multi parameter water quality analyzer (Little Swan, China) was used for the measurement of pH, dissolved oxygen and total hardness of water. The livers from 24 individuals of experimental and control groups were dissected out separately after ether anesthetized at 4 weeks after
2.4. Identification of differentially expressed genes (DEGs) The fragments per kilobase per million reads (FPKM) algorithm was used to quantify the transcript abundance in the differentially expressed genes (DEGs) analyses (Mortazavi et al., 2008). The DEGs were identified using an algorithm developed by Audic and Claverie (1997). A false discovery rate (FDR) ≤ 0.001 was used as the threshold for the pvalue in multiple tests to judge the significance of the gene expression differences. When the p-value ≤ 0.001 and a greater than two-fold change (absolute value of log2 ratio > 1) was observed in the gene expression, the transcript differences were considered significant. 2.5. Validation of RNA-seq profiles by qPCR To see the dynamic expression process of the different genes during trichlorfon stress, four differentially transcribed genes related to metabolic pathways and metabolism of xenobiotics were examined, 488
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
including induced and repressed genes. The four candidate genes were measured in the same RNA samples for validation by qPCR. Total RNA was extracted from the target hepatic tissues from the control and treated groups using the Trizol method (Invitrogen, CA, USA), as mentioned above. First-stand cDNA was synthesized using the PrimeScript RT reagent Kit with DNA eraser according to the manufacturer's protocol (Takara, Japan). The candidate gene-specific primers were designed with Primer Premier 6.0 (Supplementary file 1). The rpl8 gene was used as an internal control to normalize the targets for quantification. The qPCR analysis was conducted on an IQ5 Real Time PCR Detection system (Bio-Rad Laboratories, Richmond, CA) using SYBR Premix Ex Taq (Takara, Japan) according to the manufacturer's recommendations. The reaction mixtures were 20 μL in total volume and included 10 μL 2 × SYBR Permix Ex Taq, 0.4 μL each of the 10 μM forward and reverse primers, 2 μL cDNA, and RNase-free ddH2O for a volume of 20 μL. Three technical replicates and three biological replicates were conducted in each assay. The mRNA relative expression levels were calculated using the comparative CT method for relative quantification. Relative gene expression levels were evaluated using the 2-△△Ct method (Livak and Schmittgen, 2001). All data are expressed as the mean ± SD. The values obtained were compared based on TukeyKramer HSD test using Graph Pad Prism 5.0 (GraphPad Software Inc., USA). p < 0.05 was regarded as statistically significant.
Table 1 Summary of all unigenes annotated in the liver from R. chensinensis.
All unigenes Annotated using Nt database Annotated using Nr database Annotated using Swiss-Prot database Annotated using COG database Annotated using GO database Annotated using KEGG database All annotated unigenes
Number
Ratio
55,235 21,719 22,888 20,934 7375 15,631 16,923 27,781
100% 39.32% 41.44% 37.90% 13.35% 28.30% 30.64% 50.30%
3. Results 3.1. Sequencing and de novo assembly A total of 37,803,562 raw reads were produced from the liver in the control group, compared with 37,836,552 raw nucleotide reads from the treated group. After quality-trimming and adaptor-clipping, 37,015,892 and 37,178,282 clean reads remained for the control and treated groups, respectively. The Q20 and GC contents for the control and treated livers were 96.24%, 96.33%, 44.67%, and 45.51%, respectively. The unknown bases were 0.00%. Through assembly and elimination of redundancy, 105,047 contigs with an average length of 339 bp (N50 = 559 bp) for a total length of 35,600,355 nt were generated from the control group, whereas 103,523 contigs with an average length of 343 bp (N50 = 576) for a total length of 35,508,761 nt were created for the trichlorfon stress group. After splicing and redundancy removal, the contigs were further assembled into 55,235 unigenes with an average length of 756 bp (N50 = 1245 bp). 3.2. Functional annotation of unigenes The unigenes were annotated based on the identity of the translation frame, and 55,235 unigenes were predicted to play roles in mapping protein-coding sequences. Unigene sequences were aligned by Blast-n to the nt nucleotide database (E-value < 1e-5), and annotated by Blast-x to protein databases, such as Nr, Swiss-Prot, COG, GO and KEEG (E-value < 1e-5). A total of 21,719 (39.32%) unigenes were unambiguously aligned to the reference when Blast-n used NCBI's nt database, whereas Blast-x against NCBI's nr database or the SwissProt database returned 22,888 (41.44%) and 20,934 (37.90%) hits, respectively. Specifically, the size of the unigenes has a marked impact on the annotation efficiency (Table 1). Additionally, the threshold E-value of the annotated unigenes against the Nr database was 1e-5 (Fig. 1A). Only 15.13% of the unigenes have strong similarity (less than 1e-100) in the Nr database, whereas 84.87% of the unigenes ranged from 1e-5 to 1e-100. The similarity distributions demonstrated that approximately 40.28% of the sequences have a similarity greater than 80% and that the remaining 59.72% of the sequences have a similarity ranging from 17% to 80% (Fig. 1B). For species distributions matched against the Nr database, 53.16% of the matched unigenes showed similarities with Xenopus tropicalis, followed by Xenopus laevis (19.38%) and Rana catesbeiana
Fig. 1. Characteristics from the homology search with the Illumina sequences against the nr database. A, E-value distribution. B, similarity distribution. C, species distribution.
(2.63%) (Fig. 1C). More unigenes were similar to X. tropicalis and X. laevis because of their closer phylogenetic relationship and their abundant genomic information. Furthermore, the remaining 24.83% of matched unigenes showed similarities to other species due to the absence of genome information in amphibians. For this reason, the unigenes in the liver from R. chensinensis in this study should be further annotated with the flurry of published crustaceans genes and with more available gene background information. COG assignments were performed to predict and classify the possible unigene functions (Fig. 2). Analysis of the COG annotations showed that 7375 (13.35%) hits were assigned to the 55,235 unigenes. The hits from the COG prediction were functionally classified into 26 489
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
Fig. 2. Distribution of all unigenes annotated with COG terms.
proportion of the molecular functions (Fig. 3).
categories, in which the most enriched terms were general function (2937 unigenes), followed by replication, recombination and repair (1323 unigenes), and translation, ribosomal structure and biogenesis (1252 unigenes). As mentioned above, GO analysis was performed by Blast2GO on the putative proteins (Götz et al., 2008). After parsing the GO annotations output, 15,631 (28.3%) GO terms were finally obtained, with 52.18% for biological process, 34.83% for cellular components and 12.99% for molecular functions. The three main GO categories were classified into 58 subcategories. The distribution of the GO terms showed that cellular process, metabolic process and single-organism process were significantly enriched terms among biological processes, that the cell and cell part were the most represented terms in cellular components and that binding and catalytic activity formed the largest
3.3. KEGG pathway analysis In addition to GO and COG analyses, KEGG pathway paring was carried out. The annotations by KEGG provided various molecular pathways (Table 2). A total of 16,923 (30.64%) unigenes were mapped to 257 KEGG signalling pathways, including metabolic pathways, endocytosis, MAPK signalling pathway, and PPAR signalling pathway. The well-represented top three pathways with the greatest number of annotated sequences were metabolic pathways, pathway in cancer and regulation of actin cytoskeleton (Table 2). The number of sequences annotated in KEGG pathways ranged from 1 to 2087. These annotations provided valuable information for studying the specific biological and
Fig. 3. Gene Ontology of the unigenes annotated with the GO terms.
490
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
Table 2 The top 20 pathways with the most annotated sequences in the liver from R. chensinensis. Number
Pathways
All genes with pathway annotation (16,923)
Ratio
Pathway ID
1 2 3
Metabolic pathways Pathways in cancer Regulation of actin cytoskeleton Influenza A Focal adhesion Epstein-Barr virus infection Purine metabolism Herpes simplex infection HTLV-I infection Endocytosis RNA transport Ubiquitin mediated proteolysis Phagosome Tuberculosis Huntington's disease MAPK signalling pathway Tight junction Protein processing in endoplasmic reticulum Transcriptional misregulation in cancer Fc gamma R-mediated phagocytosis
2087 755 577
12.33% 4.46% 3.41%
ko01100 ko05200 ko04810
560 553 543
3.31% 3.27% 3.21%
ko05164 ko04510 ko05169
520 490
3.07% 2.9%
ko00230 ko05168
485 467 444 433
2.87% 2.76% 2.62% 2.56%
ko05166 ko04144 ko03013 ko04120
421 420 415 404
2.49% 2.48% 2.45% 2.39%
ko04145 ko05152 ko05016 ko04010
400 377
2.36% 2.23%
ko04530 ko04141
373
2.2%
ko05202
368
2.17%
ko04666
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Fig. 5. The relative expression of four candidate genes involved in metabolic pathways by qPCR. These data are expressed as the mean ± SD relative to the control. *p < 0.05 and ***p < 0.001.
3.5. Validation of RNA-seq profiles by qPCR To test and verify the credibility of the RNA-seq results, qPCR experiments were performed for validation. Four candidate genes involved in metabolic pathways and metabolism of xenobiotics, which included both induced and repressed genes, were selected to perform qPCR and be compared with the RNA-seq data at the gene transcript level. The results showed that the expression patterns for the four DEGs corroborated with the RNA-seq-based regulation analysis (Fig. 5). During 4 weeks of trichlorfon treatment, unigene 22540 and CL307.contig2 were significantly up-regulated, implying that these genes were involved in the response to trichlorfon. Similarly, a significant decrease in the expression of the encoding genes CL1105.contig1 and CL3031.contig2 were shown after 4 weeks of trichlorfon treatment.
metabolic process, functions and molecular mechanisms in the followup studies. 3.4. Identification of DEGs induced by trichlorfon stress To identify genes that displayed significant changes in expression during trichlorfon stress, DEGs were analysed by comparing the treated and control libraries. In all, the abundance of 3329 filtered unigenes (FDR ≤ 0.001 and log2Ratio ≥ 1) of varying overall length showed evidence of differential transcription. A higher number of genes were down-regulated than up-regulated in response to trichlorfon stress. From the data sets, a total of 1444 up-regulated and 1885 down-regulated DEGs were detected after exposure to trichlorfon for 4 weeks (Fig. 4).
4. Discussion RNA-seq technology is an ideal method for high-throughput sequence determination, which has provided efficiency and speed in gene discovery. In this study, RNA-seq was used to evaluate the transcriptional responses in the livers of R. chensinensis exposed to trichlorfon. A
Fig. 4. DEGs in the liver from R. chensinensis under trichlorfon stress. A, The number of DEGs identified in each library contrast by applying a threshold of the ratio change ≥ 2 and a q-value of < 0.05. The red and green columns represent genes up- and down-regulated by trichlorfon stress, respectively. B, Scatter plot of DEGs (FDR ≤ 0.001 and |log2Ratio| ≥ 1) illustrating the full set of genes samples. Red points are up-regulated genes, blue points are non-DEGs, and green points are down-regulated genes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
491
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
number of genes that encode P450 enzymes, which are involved in the detoxification of exogenous chemicals and the metabolism of endogenous substrates (Nebert et al., 2004). The CYP 1–3 families are involved in xenobiotic metabolism and are expressed primarily in the liver (Kawai et al., 2013). CYP gene over-expression has been attributed to pesticide resistance in several insects (Daborn et al., 2001; Yang et al., 2006; Komagata et al., 2010). These results revealed that trichlorfon increased the CL307.contig2 (CYP2C) expression levels in the liver, which could lead to enhanced trichlorfon metabolism. Meanwhile, there is strong evidence suggesting that the activation of CYP2C epoxygenase in endothelial cells is an essential step in the NO-independent vasodilatation of several vascular beds (Maier and Roman, 2001). This was identical to the increase in eNOS expression in the liver, which highlights the role of CYP2C in vascular homeostasis. In addition to the increase in CYP2C levels in the presence of trichlorfon, CL1105.contig1 (CYP3A) mRNA expression was decreased after exposure to trichlorfon for 4 weeks. The decreased levels of CYP3A may hinder the metabolism of trichlorfon, resulting in the accumulation of the compound and its toxicity. This result is consistent with the report that hepatic CYP3A1 mRNA levels were significantly decreased after renal ischemia/reperfusion (Wang et al., 2009a, 2009b). Xenobiotics should regulate/inhibit some key roles of CYP3A. For example, changes in the serotonin levels in the brain and catecholaminergic systems or in important transcriptional factors, such as the glucocorticoid receptor or the constitutive and rostane receptor, have been shown to regulate CYP genes (Rysz et al., 2015). GST constitutes an important supergene family involved in protecting against the deleterious effects of oxidative stress and xenobiotics. Down-regulation of CL3031.contig2 (glutathione S-transferase kappa 1) gene expression highlighted the fact that this enzyme probably plays divergent physiological roles during trichlorfon detoxification in R. chensinensis. In a previous study, the GST mRNA expression decreased over time with 500 μg L−1 triclosan for 168 h (Ku et al., 2014). The following are possible explanations for these phenomena: (1) different xenobiotics might affect transcription factors and then further affect the transcript of GST expression (Jacobs et al., 2005) and (2) the incomplete degradation of xenobiotics absorbed by the organisms. The results demonstrate that trichlorfon stress induces the expression of selected metabolism and antioxidant defence system-related genes in the liver, which could be involved in the response to trichlorfon stress. Furthermore, the results also validated the Illumina sequencing and assembly results. Although the specific functions of these genes remain ambiguous, our goal in the present study was to gain an overall understanding of R. chensinensis responses to trichlorfon stress and provide elementary insights into crucial signalling pathways processes. Follow-up studies will utilize the results obtained here as a foundation for more targeted studies, such as cloning the genes involved in metabolism and antioxidant defence systems, functional analysis of these genes and the regulation of associated signalling pathways.
total of 55,235 unigenes in the hepatic tissue from R. chensinensis were yielded from Illumina sequencing, 41.44% of which were aligned to the Nr database. The huge amount of sequencing data would be a stepping stone for the use of amphibian models in molecular toxicology. Furthermore, it supports the assumption that RNA-seq could generate more novel transcripts compared with traditional methods. 4.1. Functional annotations The functional annotations of the unigenes according to the COG, GO and KEGG databases provided ample numbers of candidate genes and more information about their biological features in R. chensinensis under trichlorfon stress. Therein, unigene 22540, CL307. Contig2, and CL1105. Contig1 were classified into metabolic pathways, whereas CL3031. Contig2, CL307. Contig2, and CL1105. Contig1 were classified as the metabolism of xenobiotics and drugs. Additionally, antioxidant genes including endothelial nitric oxide synthase (eNOS) and glutathione S-transferase (GST) were also identified in this transcriptome database. Previous research also confirmed that aqueous trichlorfon induced oxidative stress, affected intermediary metabolism, and subsequently increased susceptibility to pathogen infections (Venturini et al., 2015). All of this information suggests that metabolism and antioxidant defence systems should be activated to protect amphibians against pathogen infections under trichlorfon stress. We also obtained an abundance of transcripts related to detoxification. Cytochrome p450 monooxygenases (P450s) and GSTs are protein families that metabolize a wide range of xenobiotic (such as pesticides, toxic chemicals and hydrocarbons) and endogenous compounds, which may be drugs or products of oxidative stress (Ranson et al., 2002; Contreras-Vergara et al., 2004). In this study, our transcriptome data identified 295 hits coding for P450s, which were the most abundant sequences among the detoxification-related genes. There were 155 hits encoding P450s in metabolism of xenobiotics and 140 hits encoding P450s related to drug metabolism. Further, 89 hits encoded other enzymes in drug metabolism and 20 hits encoded GSTs. The results showed a consistency between our study and a transcriptome study in Nilaparvata lugens, as these proteins have been found in virtually all organisms (Contreras-Vergara et al., 2004). There were also 19 hits encoding eNOS genes, which are associated with the production of reactive oxygen species (ROS) during oxidative stress. This suggests that general physiological processes, such as enzymatic activities and oxidative metabolism, are likely disturbed in aquatic organisms exposed to organophosphorus pesticides. 4.2. Validation of RNA-seq profiles via qPCR The four unigenes subjected to qPCR analysis were selected based on the results from the annotations. The mRNA expression levels of these genes involved in metabolism and antioxidant defence systems were all significantly changed after exposure to trichlorfon for 4 weeks. Overall, the similar gene expression profiles obtained from the two distinct experiments, using different methods, showed that trichlorfon stresses induced a reproducible response in the liver of R. chensinensis. These results also further confirmed the reliability of the RNA-seq data and the accuracy of the Trinity assembly. Unigene 22540 is a highly regulated eNOS that derives NO in the liver and is essential for regulating normal hepatic blood flow (Palmer et al., 1987). The bioavailability of nitric oxide (NO) can be affected by cellular metabolism, including increases in the generation of reactive oxygen species, which will reduce NO bioavailability and in turn affect eNOS expression. Altered eNOS expression in the liver is an important factor in the progression of liver fibrogenesis to stress (Leung et al., 2008). The increased eNOS expression in this observation is in agreement with previous findings in brain tissues from rats after treatment with nebivolol (Heeba and El-Hanafy, 2012). The cytochrome p450 (CYP) gene superfamily consists of a large
5. Conclusion In this report, we have determined the transcriptomic response of R. chensinensis under trichlorfon stress with RNA-seq to examine the liver transcriptome profile and identify DEGs. Metabolic pathways-related transcripts were detected in this dataset. The results indicated that these DEGs are implicated in response to trichlorfon stress. Furthermore, the bioinformatics analysis of the transcriptome in this study generated a large number of novel genes, which will be of great value for understanding trichlorfon's mechanisms of toxicity in nonmodel species. Conflict of interests The first two authors made equivalent contributions to this work. 492
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
All of the authors read and approved the final manuscript. The authors declare that there are no conflicts of interest regarding the publication of this paper.
detoxification system in the yellow catfish (Pelteobagrus fulvidraco): expressions of CYP and GST genes and corresponding enzyme activity in phase I, II and antioxidant system. Comp. Biochem. Phys. C. 166, 105–114. Leung, T.M., Tipoe, G.L., Liong, E.C., Lau, T.Y., Fung, M.L., Nanji, A.A., 2008. Endothelial nitric oxide synthase is a critical factor in experimental liver fibrosis. Int. J. Exp. Pathol. 89 (4), 241–250. Li, B., Ma, Y., Zhang, Y.H., 2017. Oxidative stress and hepatotoxicity in the frog, Rana chensinensis, when exposed to low doses of trichlorfon. J. Environ. Sci. Health B 52 (7), 476–482. Li, X.Y., Xiao, N., Zhang, Y.H., 2014. Toxic effects of octylphenol on the expression of genes in liver identified by suppression subtractive hybridization of Rana chensinensis. Ecotoxicol 23 (1), 1–10. Li, X.Y., Liu, J., Zhang, Y.H., 2016. Octylphenol induced gene expression in testes of Frog, Rana chensinensis. Ecotoxicol. Environ. Saf. 128, 75–82. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25 (4), 402–408. Lundberg, R., Jenssen, B.M., Leiva-Presa, A., Rönn, M., Hernhag, C., Wejheden, C., Larsson, S., Orberg, J., Lind, P.M., 2007. Effects of short-term exposure to the DDT metabolite p,p′-DDE on bone tissue in male common frog (Rana temporaria). J. Toxicol. Environ. Heal. A. 70, 614–619. Maier, K.G., Roman, R.J., 2001. Cytochrome P450 metabolites of arachidonic acid in the control of renal function. Curr. Opin. Nephrol. Hypertens. 10 (1), 81–87. Mataqueiro, M.I., Nakaghi, L.S., Souza, J.P., Cruz, C.D., Oliveira, G.H., Urbinati, E.C., 2009. Histopathological changes in the gill, liver and kidney of pacu (Piaractus mesopotamicus, Holmberg, 1887) exposed to various concentrations of trichlorfon. J. Appl. Ichthyol. 25 (1), 124–127. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5 (7), 621–628. Nebert, D.W., Dalton, T.P., Okey, A.B., Gonzalez, F.J., 2004. Role of aryl hydrocarbon receptor-mediated induction of the CYP1 enzymes in environmental toxicity and cancer. J. Biol. Chem. 279 (23), 23847–23850. Palmer, R.M., Ferrige, A.G., Moncada, S., 1987. Nitric oxide release accounts for the biological activity of endothelium-derived relaxing factor. Nature 327 (6122), 524–526. Pruitt, K.D., Tatusova, T., Maglott, D.R., 2007. NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucl. Acids Res. 35, D61–D65. Qiao, L., Yang, W., Fu, J., Song, Z., 2013. Transcriptome profile of the green odorous frog (Odorrana margaretae). PLoS One 8 (9), e75211. Ranson, H., Claudianos, C., Ortelli, F., Abgrall, C., Hemingway, J., Sharakhova, M.V., Unger, M.F., Collins, F.H., Feyereisen, R., 2002. Evolution of supergene families associated with insecticide resistance. Science 298 (5591), 179–181. Rysz, M., Bromek, E., Daniel, W.A., 2015. Activation of the brain serotonergic system decreases the expression and activity of liver cytochrome P450 isoforms CYP1A, 2C11 and 3A. Eur. Neuropsychopharm 25 (S2), 209–210. Soumis, N., Lucotte, M., Sampaio, D., Almeida, D.C., Giroux, D., Morais, S., Pichet, P., 2003. Presence of organophosphate insecticides in fish of the Amazon River. Acta Amazon. 33, 325–337. Thomaz, J.M., Martins, N.D., Monteiro, D.A., Rantin, F.T., Kalinin, A.L., 2009. Cardiorespiratory function and oxidative stress biomarkers in Nile tilapia exposed to the organophosphate insecticide trichlorfon (NEGUVON). Ecotoxicol. Environ. Saf. 72 (5), 1413–1424. Venturini, F.P., Moraes, F.D., Cortella, L.R., Rossi, P.A., Cruz, C., Moraes, G., 2015. Metabolic effects of trichlorfon (Masoten®) on the neotropical freshwater fish pacu (Piaractus mesopotamicus). Fish. Physiol. Biochem. 41 (1), 299–309. Wang, B.Y., Li, Q.X., Li, J., Xie, X.F., Ao, Y., Ai, Y.X., 2009a. Hepatotoxicity and gene expression down-regulation of CYP isozymes caused by renal ischemia/reperfusion in the rat. Exp. Toxicol. Pathol. 61 (2), 169–176. Wang, Z., Gerstein, M., Snyder, M., 2009b. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10 (1), 57–63. World Health Organization (WHO), 1992. Environmental Health Criteria 132: Trichlorfon, International Program on Chemical Safety. WHO, Geneva. Xu, W., Liu, W., Lu, K., Jiang, Y., Li, G., 2012. Effect of trichlorfon on oxidative stress and hepatocyte apoptosis of Carassius auratus gibelio in vivo. Fish. Physiol. Biochem. 38 (3), 769–775. Yang, Y., Chen, S., Wu, S., Yue, L., Wu, Y., 2006. Constitutive overexpression of multiple cytochrome P450 genes associated with pyrethroid resistance in Helicoverpa armigera. J. Econ. Entomol. 99 (5), 1784–1789. Zhang, M., Li, Y., Yao, B., Sun, M., Wang, Z., Zhao, Y., 2013. Transcriptome sequencing and de novo analysis for Oviductus Ranae of Rana chensinensis using illumina RNA-Seq technology. J. Genet. Genom. 40 (3), 137–140.
Acknowledgements This study was supported by the CAS “Light of West China” visiting scholar training project and the Foundation of the Shaanxi Academy of Science (2017k-05, 2017k-27). Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ecoenv.2017.09.016. References Audic, S., Claverie, J.M., 1997. The significance of digital gene expression profiles. Genome Res. 7 (10), 986–995. Boeckmann, B., Bairoch, A., Apweiler, R., Blatter, M.C., Estreicher, A., Gasteiger, E., Martin, M.J., Michoud, K., O'Donovan, C., Phan, I., Pilbout, S., Schneider, M., 2003. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucl. Acids Res. 31 (1), 365–370. Chang, C.C., Rahmawaty, A., Chang, Z.W., 2013. Molecular and immunological responses of the giant freshwater prawn, Macrobrachium rosenbergii, to the organophosphorus insecticide, trichlorfon. Aquat. Toxicol. 130–131, 18–26. Coelho, S., Oliveira, R., Pereira, S., Musso, C., Domingues, I., Bhujel, R.C., Soares, A.M., Nogueira, A.J., 2011. Assessing lethal and sub-lethal effects of trichlorfon on different trophic levels. Aquat. Toxicol. 103 (3–4), 191–198. Contreras-Vergara, C.A., Harris-Valle, C., Sotelo-Mundo, R.R., Yepiz-Plascencia, G., 2004. A mu-class glutathione S-transferase from the marine shrimp Litopenaeus vannamei: molecular cloning and active-site structural modeling. J. Biochem. Mol. Toxicol. 18 (5), 245–552. Daborn, P., Boundy, S., Yen, J., Pittendrigh, B., ffrench-Constant, R., 2001. DDT resistance in Drosophila correlates with Cyp6g1 over-expression and confers cross-resistance to the neonicotinoid imidacloprid. Mol. Genet. Genom. 266 (4), 556–563. Ekblom, R., Galindo, J., 2011. Applications of next generation sequencing in molecular ecology of non-model organisms. Heredity 107 (1), 1–15. Götz, S., García-Gómez, J.M., Terol, J., Williams, T.D., Nagaraj, S.H., Nueda, M.J., Robles, M., Talón, M., Dopazo, J., Conesa, A., 2008. High-throughput functional annotation and data mining with the Blast2GO suite. Nucl. Acids Res. 36 (10), 3420–3435. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29 (7), 644–652. Heeba, G.H., El-Hanafy, A.A., 2012. Nebivolol regulates eNOS and iNOS expressions and alleviates oxidative stress in cerebral ischemia/reperfusion injury in rats. Life Sci. 90 (11–12), 388–395. Herwig, N., 1979. Handbook of drugs and chemicals used in the treatment of fish diseases. In: Charles, C. (Ed.), A Manual of Fish Pharmacology and Material Medica. Springfield, Illinois, USA. Jacobs, M.N., Nolan, G.T., Hood, S.R., 2005. Lignans, bacteriocides and organochlorine compounds activate the human pregnane X receptor (PXR). Toxicol. Appl. Pharm. 209 (2), 123–133. Kanehisa, M., Goto, S., 2000. KEGG: kyoto encyclopedia of genes and genomes. Nucl. Acids Res. 28 (1), 27–30. Kawai, Y.K., Watanabe, K.P., Ishii, A., Ohnuma, A., Sawa, H., Ikenaka, Y., Ishizuka, M., 2013. De novo sequence analysis of cytochrome P450 1-3 genes expressed in ostrich liver with highest expression of CYP2G19. Comp. Biochem. Phys. D 8 (3), 201–208. Komagata, O., Kasai, S., Tomita, T., 2010. Overexpression of cytochrome P450 genes in pyrethroid-resistant Culex quinquefasciatus. Insect Biochem. Mol. 40 (2), 146–152. Koonin, E.V., Fedorova, N.D., Jackson, J.D., Jacobs, A.R., Krylov, D.M., Makarova, K.S., Mazumder, R., Mekhedov, S.L., Nikolskaya, A.N., Rao, B.S., Rogozin, I.B., Smirnov, S., Sorokin, A.V., Sverdlov, A.V., Vasudevan, S., Wolf, Y.I., Yin, J.J., Natale, D.A., 2004. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol. 5 (2), R7. Ku, P., Wu, X., Nie, X., Ou, R., Wang, L., Su, T., Li, Y., 2014. Effects of triclosan on the
493
Contents lists available at ScienceDirect
Ecotoxicology and Environmental Safety journal homepage: www.elsevier.com/locate/ecoenv
Transcriptome analysis of Rana chensinensis liver under trichlorfon stress Yu Ma a b c
a,b
a,b
, Bo Li
b
c
, Yang Ke , Yongan Zhang , Yuhui Zhang
a,⁎
MARK
College of Life Science, Shaanxi Normal University, Xi’an 710062, China Shaanxi Microbiology Institute, Xi’an 710043, China Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan 430072, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Organophosphate insecticides Trichlorfon Differentially expressed genes Rana chensinensis RNA-seq
Trichlorfon is a selective organophosphate insecticide that is widely applied in aquaculture and agriculture for control of various parasites. However, repeated and excess applications of trichlorfon often lead to water pollution and threaten non-targeted species. Our previous studies showed that trichlorfon could cause oxidative stress, lipid peroxidation and hepatic lesions in the liver of Rana chensinensis, but the related molecular mechanisms remain unclear. To explore the interference of trichlorfon in gene transcription, the differentially expressed genes in the liver of R. chensinensis exposed to trichlorfon were characterized using the RNA-seq platform. A search of all unigenes against non-redundant protein sequence (Nr), non-redundant nucleotide (Nt), Swiss-Prot, Kyoto Encyclopaedia of Genes and Genomes (KEGG), Clusters of Orthologous Groups (COG) and Gene Ontology (GO) databases resulted in 22,888, 21,719, 20,934, 16,923, 7375 and 15,631 annotations, respectively, and provided a total of 27,781 annotated unigenes. Among the annotated unigenes, 16,923 were mapped to 257 signalling pathways. A set of 3329 differentially expressed unigenes was identified by comparison of the two groups in liver. Notably, relative expression of metabolism-related genes, including both upand down-regulated genes, were also validated by qPCR. The present study depicts the high degree of transcriptional complexity in R. chensinensis under trichlorfon stress and provides new insights into the molecular mechanisms of organophosphate insecticide toxicology. Some of these metabolism-responsive genes could be useful for understanding the toxicological mechanism of trichlorfon on non-target aquatic organisms and will contribute to the conservation of aquatic life.
1. Introduction Organophosphates have been used as insecticide, acaricide and anthelmintic compounds for several decades. Trichlorfon is a selective organophosphate insecticide used for control of various parasitic infestations in fish and agriculture. The dosage of trichlorfon applied to aquaculture ponds to eradicate ectoparasites varies from 0.1 to 1.0 mg/ L (Herwig, 1979; WHO, 1992). Trichlorfon has been widely used in aquaculture, however, the trichlorfon-polluted waters from the fish and agriculture industries are potential threats to non-targeted species, such as fish, crabs, and shrimps (Soumis et al., 2003; Chang et al., 2013). Therefore, the expansion of aquaculture worldwide has generated concerns about its impact on both public and environmental health (Coelho et al., 2011). Because the effects of environmental concentrations very rarely cause lethality, the risk assessment of trichlorfon as a long-term pollutant in aquatic ecosystems remains to be comprehensively investigated in a specific species. The liver plays a critical role in the detoxification of xenobiotics.
⁎
Biochemical and histopathological data have shown the hepatotoxic tendency of trichlorfon in aquatic organisms (Mataqueiro et al., 2009; Xu et al., 2012). The metabolic changes observed in aquatic animals exposed to trichlorfon can create widespread disturbances in enzymatic activities and oxidative metabolism in general physiological processes (Thomaz et al., 2009; Coelho et al., 2011; Xu et al., 2012; Chang et al., 2013). These studies suggest that the metabolism of xenobiotics is a key event in trichlorfon-exposed organisms. However, the molecular mechanism of the metabolic pathways involved in xenobiotic detoxification induced by trichlorfon is still unclear. Previous studies have found that superoxide dismutase (SOD) gene expression levels were up-regulated after exposure to 0.4 mg/L trichlorfon for 6 h in the giant freshwater prawn (Macrobrachium rosenbergii) (Chang et al., 2013). Indeed, many more genes may participate in this process and need to be identified in more species. The development of techniques such as RNA-seq have provided a better understanding of the molecular mechanisms of toxicological responses in non-model organisms. RNA-seq is a new method for both
Corresponding author. E-mail address: [email protected] (Y. Zhang).
http://dx.doi.org/10.1016/j.ecoenv.2017.09.016 Received 8 March 2017; Received in revised form 5 September 2017; Accepted 7 September 2017 0147-6513/ © 2017 Elsevier Inc. All rights reserved.
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
exposure. Then, the tissue samples were immediately collected in cryotubes and stored in liquid nitrogen for later RNA extraction and qPCR. All procedures involving the handling and treatment of frogs used in this study were conducted in accordance with the approval of the Animal Care and Use Committee of Shaanxi Normal University.
mapping and quantifying the transcriptome; it has clear advantages over the existing approaches and improves our understanding of the extent and complexity of the eukaryotic transcriptome including nonmodel organisms that lack completely sequenced genomes (Wang et al., 2009a, 2009b; Ekblom and Galindo, 2011). RNA-seq has been recently used in some amphibian species, including Odorrana margaretae (Qiao et al., 2013) and Rana chensinensis (Zhang et al., 2013), but the use of the high-throughput RNA-seq technology for studying environmental stressors is still limited. This study aims to understand the molecular mechanism of trichlorfon impact on liver function in amphibians. Amphibians spend the first phase of the life cycle in aquatic environments and their highly permeable skin makes them particularly sensitive to contaminants. Furthermore, their low mobility and high diversity in their reproductive modes through critical hormone-related developmental stages in aquatic environments render them vulnerable to environmental changes. The Chinese brown frog (Rana chensinensis) is a unique amphibian species that is widely distributed in northern China. This species is considered a potential candidate indicator due to its high sensitivity to chemical pollution in an aquatic environment; thus, it can be conveniently used as a native species to research ecotoxicology in local amphibians (Lundberg et al., 2007; Li et al., 2014, 2016). Our previous studies showed that trichlorfon can cause oxidative stress, lipid peroxidation and hepatic lesions in the liver of R. chensinensis (Li et al., 2017). These findings suggest that the liver is the major organ targeted by trichlorfon. However, the related molecular mechanisms remain unclear. The present study aims to obtain transcriptome information for R. chensinensis exposed to trichlorfon by using newly developed RNA-seq technology. The mRNA expression levels of four genes involved in metabolism and detoxification (eNOS, Cyt 2c, Gst kappa I,and Cyt 3a) were examined using quantitative real-time PCR (qPCR). The resulting data provides a preliminary research direction for the molecular mechanisms underlying environmental exposure, enriches the information currently in the cDNA database and identifies the genes involved in the metabolism and detoxification of trichlorfon toxicity-related mechanisms in amphibians.
2.2. RNA isolation and deep sequencing Total RNA was isolated from liver tissues from trichlorfon-exposed and unexposed frogs at 4 weeks using a TRIzol Kit (Invitrogen, CA, USA) according to the manufacturer's instructions. To obtain as many gene transcripts as possible and avoid the negative effects of individual variation, equal volumes of high-quality RNA from 12 individuals from each group were then pooled together for transcriptome sequencing. The RNA integrity was checked using an Agilent 2100 BioAnalyzer (Agilent Technologies, CA, USA), and the mRNA was purified from the total RNA using the PolyATract mRNA Isolation Systems Kit (Promega) following the manufacturer's instructions. Fragmentation was carried using divalent cations under elevated temperature in NEB Next First Strand Synthesis Reaction Buffer. First stand cDNA was synthesized using random hexamer primer and M-MuLV reverse transcriptase (RNase H-), and then the second strand cDNA synthesis was performed using DNA polymerase I and RNase H. In order to select cDNA fragments that were preferentially 150–200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). The adaptor modified fragments were selected by gel purification and amplified through PCR to create the final cDNA library. Transcriptome sequencing was conducted using an Illumina Hiseq™ 2000 according to the manufacturer's instructions. The invalid reads which contain adapters or unknown or low-quality bases, were abandoned to avoid negatively affecting the bioinformatics analysis. De novo transcriptome assembly was carried out using the short-read assembly program Trinity software (Grabherr et al., 2011). 2.3. Bioinformatic analysis
2. Materials and methods
Unigenes from the de novo assembly were further processed by sequence annotation using the NCBI (http://blast.ncbi.nlm.nih.gov/) non-redundant protein sequence (Nr) and non-redundant nucleotide (Nt) (Pruitt et al., 2007), Swiss-Prot (http://www.uniprot.org/) (Boeckmann et al., 2003), Clusters of Orthologous Groups (COG) (http://www.ncbi.nlm.nih.gov/COG/) (Koonin et al., 2004), Gene Ontology database (GO) (http://www.geneontology.org/) (Götz et al., 2008), and Kyoto Encyclopaedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/pathway.html) databases (Kanehisa and Goto, 2000). Gene names were assigned to each protein sequence based on the best Blast hit from all the Blast results.
2.1. Animals and treatments Adult Chinese brown frogs, Rana chensinensis, were obtained from the rivulet in unpolluted areas in Qinling Mountains (Shaanxi, China; N34°0′30.95′′, E109°06′48.64′′). The frogs were acclimatized in aquariums for 7 days prior to the experiment. After that, a total of 48 healthy frogs of both sexes, weighing approximately 8.967 ± 1.383 g, were randomly sampled and divided into experimental and control groups. Based on the preliminary experiment and the applied dose of trichlorfon in aquaculture, the exposure trichlorfon dose at 0.1 mg/L (90% purity, Hubei, China) was selected, and a blank control group was set. Each aquarium (40 cm × 20 cm × 20 cm) had 8 individuals with 2 L of test solution for static-renewal assay, and there were three replicates for each treatment. During the acclimation and experimental periods, the frogs were fed worms (Tenebrio molitor) three times a week. The tanks were placed on a slant to provide the option of both aqueous and dry environments. Trichlorfon-treated water was fresh prepared and renewed every 48 h to keep the trichlorfon concentration constant, and the tank was cleaned thoroughly. Throughout the experiment, the water temperature was maintained at 20 ± 2 °C, the pH values ranged from 6.9 to 7.1, the dissolved oxygen concentration remained at 7.0 ± 0.1 mg/L, and the photoperiod was set as 12:12 light:dark. GDYS-201 M multi parameter water quality analyzer (Little Swan, China) was used for the measurement of pH, dissolved oxygen and total hardness of water. The livers from 24 individuals of experimental and control groups were dissected out separately after ether anesthetized at 4 weeks after
2.4. Identification of differentially expressed genes (DEGs) The fragments per kilobase per million reads (FPKM) algorithm was used to quantify the transcript abundance in the differentially expressed genes (DEGs) analyses (Mortazavi et al., 2008). The DEGs were identified using an algorithm developed by Audic and Claverie (1997). A false discovery rate (FDR) ≤ 0.001 was used as the threshold for the pvalue in multiple tests to judge the significance of the gene expression differences. When the p-value ≤ 0.001 and a greater than two-fold change (absolute value of log2 ratio > 1) was observed in the gene expression, the transcript differences were considered significant. 2.5. Validation of RNA-seq profiles by qPCR To see the dynamic expression process of the different genes during trichlorfon stress, four differentially transcribed genes related to metabolic pathways and metabolism of xenobiotics were examined, 488
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
including induced and repressed genes. The four candidate genes were measured in the same RNA samples for validation by qPCR. Total RNA was extracted from the target hepatic tissues from the control and treated groups using the Trizol method (Invitrogen, CA, USA), as mentioned above. First-stand cDNA was synthesized using the PrimeScript RT reagent Kit with DNA eraser according to the manufacturer's protocol (Takara, Japan). The candidate gene-specific primers were designed with Primer Premier 6.0 (Supplementary file 1). The rpl8 gene was used as an internal control to normalize the targets for quantification. The qPCR analysis was conducted on an IQ5 Real Time PCR Detection system (Bio-Rad Laboratories, Richmond, CA) using SYBR Premix Ex Taq (Takara, Japan) according to the manufacturer's recommendations. The reaction mixtures were 20 μL in total volume and included 10 μL 2 × SYBR Permix Ex Taq, 0.4 μL each of the 10 μM forward and reverse primers, 2 μL cDNA, and RNase-free ddH2O for a volume of 20 μL. Three technical replicates and three biological replicates were conducted in each assay. The mRNA relative expression levels were calculated using the comparative CT method for relative quantification. Relative gene expression levels were evaluated using the 2-△△Ct method (Livak and Schmittgen, 2001). All data are expressed as the mean ± SD. The values obtained were compared based on TukeyKramer HSD test using Graph Pad Prism 5.0 (GraphPad Software Inc., USA). p < 0.05 was regarded as statistically significant.
Table 1 Summary of all unigenes annotated in the liver from R. chensinensis.
All unigenes Annotated using Nt database Annotated using Nr database Annotated using Swiss-Prot database Annotated using COG database Annotated using GO database Annotated using KEGG database All annotated unigenes
Number
Ratio
55,235 21,719 22,888 20,934 7375 15,631 16,923 27,781
100% 39.32% 41.44% 37.90% 13.35% 28.30% 30.64% 50.30%
3. Results 3.1. Sequencing and de novo assembly A total of 37,803,562 raw reads were produced from the liver in the control group, compared with 37,836,552 raw nucleotide reads from the treated group. After quality-trimming and adaptor-clipping, 37,015,892 and 37,178,282 clean reads remained for the control and treated groups, respectively. The Q20 and GC contents for the control and treated livers were 96.24%, 96.33%, 44.67%, and 45.51%, respectively. The unknown bases were 0.00%. Through assembly and elimination of redundancy, 105,047 contigs with an average length of 339 bp (N50 = 559 bp) for a total length of 35,600,355 nt were generated from the control group, whereas 103,523 contigs with an average length of 343 bp (N50 = 576) for a total length of 35,508,761 nt were created for the trichlorfon stress group. After splicing and redundancy removal, the contigs were further assembled into 55,235 unigenes with an average length of 756 bp (N50 = 1245 bp). 3.2. Functional annotation of unigenes The unigenes were annotated based on the identity of the translation frame, and 55,235 unigenes were predicted to play roles in mapping protein-coding sequences. Unigene sequences were aligned by Blast-n to the nt nucleotide database (E-value < 1e-5), and annotated by Blast-x to protein databases, such as Nr, Swiss-Prot, COG, GO and KEEG (E-value < 1e-5). A total of 21,719 (39.32%) unigenes were unambiguously aligned to the reference when Blast-n used NCBI's nt database, whereas Blast-x against NCBI's nr database or the SwissProt database returned 22,888 (41.44%) and 20,934 (37.90%) hits, respectively. Specifically, the size of the unigenes has a marked impact on the annotation efficiency (Table 1). Additionally, the threshold E-value of the annotated unigenes against the Nr database was 1e-5 (Fig. 1A). Only 15.13% of the unigenes have strong similarity (less than 1e-100) in the Nr database, whereas 84.87% of the unigenes ranged from 1e-5 to 1e-100. The similarity distributions demonstrated that approximately 40.28% of the sequences have a similarity greater than 80% and that the remaining 59.72% of the sequences have a similarity ranging from 17% to 80% (Fig. 1B). For species distributions matched against the Nr database, 53.16% of the matched unigenes showed similarities with Xenopus tropicalis, followed by Xenopus laevis (19.38%) and Rana catesbeiana
Fig. 1. Characteristics from the homology search with the Illumina sequences against the nr database. A, E-value distribution. B, similarity distribution. C, species distribution.
(2.63%) (Fig. 1C). More unigenes were similar to X. tropicalis and X. laevis because of their closer phylogenetic relationship and their abundant genomic information. Furthermore, the remaining 24.83% of matched unigenes showed similarities to other species due to the absence of genome information in amphibians. For this reason, the unigenes in the liver from R. chensinensis in this study should be further annotated with the flurry of published crustaceans genes and with more available gene background information. COG assignments were performed to predict and classify the possible unigene functions (Fig. 2). Analysis of the COG annotations showed that 7375 (13.35%) hits were assigned to the 55,235 unigenes. The hits from the COG prediction were functionally classified into 26 489
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
Fig. 2. Distribution of all unigenes annotated with COG terms.
proportion of the molecular functions (Fig. 3).
categories, in which the most enriched terms were general function (2937 unigenes), followed by replication, recombination and repair (1323 unigenes), and translation, ribosomal structure and biogenesis (1252 unigenes). As mentioned above, GO analysis was performed by Blast2GO on the putative proteins (Götz et al., 2008). After parsing the GO annotations output, 15,631 (28.3%) GO terms were finally obtained, with 52.18% for biological process, 34.83% for cellular components and 12.99% for molecular functions. The three main GO categories were classified into 58 subcategories. The distribution of the GO terms showed that cellular process, metabolic process and single-organism process were significantly enriched terms among biological processes, that the cell and cell part were the most represented terms in cellular components and that binding and catalytic activity formed the largest
3.3. KEGG pathway analysis In addition to GO and COG analyses, KEGG pathway paring was carried out. The annotations by KEGG provided various molecular pathways (Table 2). A total of 16,923 (30.64%) unigenes were mapped to 257 KEGG signalling pathways, including metabolic pathways, endocytosis, MAPK signalling pathway, and PPAR signalling pathway. The well-represented top three pathways with the greatest number of annotated sequences were metabolic pathways, pathway in cancer and regulation of actin cytoskeleton (Table 2). The number of sequences annotated in KEGG pathways ranged from 1 to 2087. These annotations provided valuable information for studying the specific biological and
Fig. 3. Gene Ontology of the unigenes annotated with the GO terms.
490
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
Table 2 The top 20 pathways with the most annotated sequences in the liver from R. chensinensis. Number
Pathways
All genes with pathway annotation (16,923)
Ratio
Pathway ID
1 2 3
Metabolic pathways Pathways in cancer Regulation of actin cytoskeleton Influenza A Focal adhesion Epstein-Barr virus infection Purine metabolism Herpes simplex infection HTLV-I infection Endocytosis RNA transport Ubiquitin mediated proteolysis Phagosome Tuberculosis Huntington's disease MAPK signalling pathway Tight junction Protein processing in endoplasmic reticulum Transcriptional misregulation in cancer Fc gamma R-mediated phagocytosis
2087 755 577
12.33% 4.46% 3.41%
ko01100 ko05200 ko04810
560 553 543
3.31% 3.27% 3.21%
ko05164 ko04510 ko05169
520 490
3.07% 2.9%
ko00230 ko05168
485 467 444 433
2.87% 2.76% 2.62% 2.56%
ko05166 ko04144 ko03013 ko04120
421 420 415 404
2.49% 2.48% 2.45% 2.39%
ko04145 ko05152 ko05016 ko04010
400 377
2.36% 2.23%
ko04530 ko04141
373
2.2%
ko05202
368
2.17%
ko04666
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Fig. 5. The relative expression of four candidate genes involved in metabolic pathways by qPCR. These data are expressed as the mean ± SD relative to the control. *p < 0.05 and ***p < 0.001.
3.5. Validation of RNA-seq profiles by qPCR To test and verify the credibility of the RNA-seq results, qPCR experiments were performed for validation. Four candidate genes involved in metabolic pathways and metabolism of xenobiotics, which included both induced and repressed genes, were selected to perform qPCR and be compared with the RNA-seq data at the gene transcript level. The results showed that the expression patterns for the four DEGs corroborated with the RNA-seq-based regulation analysis (Fig. 5). During 4 weeks of trichlorfon treatment, unigene 22540 and CL307.contig2 were significantly up-regulated, implying that these genes were involved in the response to trichlorfon. Similarly, a significant decrease in the expression of the encoding genes CL1105.contig1 and CL3031.contig2 were shown after 4 weeks of trichlorfon treatment.
metabolic process, functions and molecular mechanisms in the followup studies. 3.4. Identification of DEGs induced by trichlorfon stress To identify genes that displayed significant changes in expression during trichlorfon stress, DEGs were analysed by comparing the treated and control libraries. In all, the abundance of 3329 filtered unigenes (FDR ≤ 0.001 and log2Ratio ≥ 1) of varying overall length showed evidence of differential transcription. A higher number of genes were down-regulated than up-regulated in response to trichlorfon stress. From the data sets, a total of 1444 up-regulated and 1885 down-regulated DEGs were detected after exposure to trichlorfon for 4 weeks (Fig. 4).
4. Discussion RNA-seq technology is an ideal method for high-throughput sequence determination, which has provided efficiency and speed in gene discovery. In this study, RNA-seq was used to evaluate the transcriptional responses in the livers of R. chensinensis exposed to trichlorfon. A
Fig. 4. DEGs in the liver from R. chensinensis under trichlorfon stress. A, The number of DEGs identified in each library contrast by applying a threshold of the ratio change ≥ 2 and a q-value of < 0.05. The red and green columns represent genes up- and down-regulated by trichlorfon stress, respectively. B, Scatter plot of DEGs (FDR ≤ 0.001 and |log2Ratio| ≥ 1) illustrating the full set of genes samples. Red points are up-regulated genes, blue points are non-DEGs, and green points are down-regulated genes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
491
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
number of genes that encode P450 enzymes, which are involved in the detoxification of exogenous chemicals and the metabolism of endogenous substrates (Nebert et al., 2004). The CYP 1–3 families are involved in xenobiotic metabolism and are expressed primarily in the liver (Kawai et al., 2013). CYP gene over-expression has been attributed to pesticide resistance in several insects (Daborn et al., 2001; Yang et al., 2006; Komagata et al., 2010). These results revealed that trichlorfon increased the CL307.contig2 (CYP2C) expression levels in the liver, which could lead to enhanced trichlorfon metabolism. Meanwhile, there is strong evidence suggesting that the activation of CYP2C epoxygenase in endothelial cells is an essential step in the NO-independent vasodilatation of several vascular beds (Maier and Roman, 2001). This was identical to the increase in eNOS expression in the liver, which highlights the role of CYP2C in vascular homeostasis. In addition to the increase in CYP2C levels in the presence of trichlorfon, CL1105.contig1 (CYP3A) mRNA expression was decreased after exposure to trichlorfon for 4 weeks. The decreased levels of CYP3A may hinder the metabolism of trichlorfon, resulting in the accumulation of the compound and its toxicity. This result is consistent with the report that hepatic CYP3A1 mRNA levels were significantly decreased after renal ischemia/reperfusion (Wang et al., 2009a, 2009b). Xenobiotics should regulate/inhibit some key roles of CYP3A. For example, changes in the serotonin levels in the brain and catecholaminergic systems or in important transcriptional factors, such as the glucocorticoid receptor or the constitutive and rostane receptor, have been shown to regulate CYP genes (Rysz et al., 2015). GST constitutes an important supergene family involved in protecting against the deleterious effects of oxidative stress and xenobiotics. Down-regulation of CL3031.contig2 (glutathione S-transferase kappa 1) gene expression highlighted the fact that this enzyme probably plays divergent physiological roles during trichlorfon detoxification in R. chensinensis. In a previous study, the GST mRNA expression decreased over time with 500 μg L−1 triclosan for 168 h (Ku et al., 2014). The following are possible explanations for these phenomena: (1) different xenobiotics might affect transcription factors and then further affect the transcript of GST expression (Jacobs et al., 2005) and (2) the incomplete degradation of xenobiotics absorbed by the organisms. The results demonstrate that trichlorfon stress induces the expression of selected metabolism and antioxidant defence system-related genes in the liver, which could be involved in the response to trichlorfon stress. Furthermore, the results also validated the Illumina sequencing and assembly results. Although the specific functions of these genes remain ambiguous, our goal in the present study was to gain an overall understanding of R. chensinensis responses to trichlorfon stress and provide elementary insights into crucial signalling pathways processes. Follow-up studies will utilize the results obtained here as a foundation for more targeted studies, such as cloning the genes involved in metabolism and antioxidant defence systems, functional analysis of these genes and the regulation of associated signalling pathways.
total of 55,235 unigenes in the hepatic tissue from R. chensinensis were yielded from Illumina sequencing, 41.44% of which were aligned to the Nr database. The huge amount of sequencing data would be a stepping stone for the use of amphibian models in molecular toxicology. Furthermore, it supports the assumption that RNA-seq could generate more novel transcripts compared with traditional methods. 4.1. Functional annotations The functional annotations of the unigenes according to the COG, GO and KEGG databases provided ample numbers of candidate genes and more information about their biological features in R. chensinensis under trichlorfon stress. Therein, unigene 22540, CL307. Contig2, and CL1105. Contig1 were classified into metabolic pathways, whereas CL3031. Contig2, CL307. Contig2, and CL1105. Contig1 were classified as the metabolism of xenobiotics and drugs. Additionally, antioxidant genes including endothelial nitric oxide synthase (eNOS) and glutathione S-transferase (GST) were also identified in this transcriptome database. Previous research also confirmed that aqueous trichlorfon induced oxidative stress, affected intermediary metabolism, and subsequently increased susceptibility to pathogen infections (Venturini et al., 2015). All of this information suggests that metabolism and antioxidant defence systems should be activated to protect amphibians against pathogen infections under trichlorfon stress. We also obtained an abundance of transcripts related to detoxification. Cytochrome p450 monooxygenases (P450s) and GSTs are protein families that metabolize a wide range of xenobiotic (such as pesticides, toxic chemicals and hydrocarbons) and endogenous compounds, which may be drugs or products of oxidative stress (Ranson et al., 2002; Contreras-Vergara et al., 2004). In this study, our transcriptome data identified 295 hits coding for P450s, which were the most abundant sequences among the detoxification-related genes. There were 155 hits encoding P450s in metabolism of xenobiotics and 140 hits encoding P450s related to drug metabolism. Further, 89 hits encoded other enzymes in drug metabolism and 20 hits encoded GSTs. The results showed a consistency between our study and a transcriptome study in Nilaparvata lugens, as these proteins have been found in virtually all organisms (Contreras-Vergara et al., 2004). There were also 19 hits encoding eNOS genes, which are associated with the production of reactive oxygen species (ROS) during oxidative stress. This suggests that general physiological processes, such as enzymatic activities and oxidative metabolism, are likely disturbed in aquatic organisms exposed to organophosphorus pesticides. 4.2. Validation of RNA-seq profiles via qPCR The four unigenes subjected to qPCR analysis were selected based on the results from the annotations. The mRNA expression levels of these genes involved in metabolism and antioxidant defence systems were all significantly changed after exposure to trichlorfon for 4 weeks. Overall, the similar gene expression profiles obtained from the two distinct experiments, using different methods, showed that trichlorfon stresses induced a reproducible response in the liver of R. chensinensis. These results also further confirmed the reliability of the RNA-seq data and the accuracy of the Trinity assembly. Unigene 22540 is a highly regulated eNOS that derives NO in the liver and is essential for regulating normal hepatic blood flow (Palmer et al., 1987). The bioavailability of nitric oxide (NO) can be affected by cellular metabolism, including increases in the generation of reactive oxygen species, which will reduce NO bioavailability and in turn affect eNOS expression. Altered eNOS expression in the liver is an important factor in the progression of liver fibrogenesis to stress (Leung et al., 2008). The increased eNOS expression in this observation is in agreement with previous findings in brain tissues from rats after treatment with nebivolol (Heeba and El-Hanafy, 2012). The cytochrome p450 (CYP) gene superfamily consists of a large
5. Conclusion In this report, we have determined the transcriptomic response of R. chensinensis under trichlorfon stress with RNA-seq to examine the liver transcriptome profile and identify DEGs. Metabolic pathways-related transcripts were detected in this dataset. The results indicated that these DEGs are implicated in response to trichlorfon stress. Furthermore, the bioinformatics analysis of the transcriptome in this study generated a large number of novel genes, which will be of great value for understanding trichlorfon's mechanisms of toxicity in nonmodel species. Conflict of interests The first two authors made equivalent contributions to this work. 492
Ecotoxicology and Environmental Safety 147 (2018) 487–493
Y. Ma et al.
All of the authors read and approved the final manuscript. The authors declare that there are no conflicts of interest regarding the publication of this paper.
detoxification system in the yellow catfish (Pelteobagrus fulvidraco): expressions of CYP and GST genes and corresponding enzyme activity in phase I, II and antioxidant system. Comp. Biochem. Phys. C. 166, 105–114. Leung, T.M., Tipoe, G.L., Liong, E.C., Lau, T.Y., Fung, M.L., Nanji, A.A., 2008. Endothelial nitric oxide synthase is a critical factor in experimental liver fibrosis. Int. J. Exp. Pathol. 89 (4), 241–250. Li, B., Ma, Y., Zhang, Y.H., 2017. Oxidative stress and hepatotoxicity in the frog, Rana chensinensis, when exposed to low doses of trichlorfon. J. Environ. Sci. Health B 52 (7), 476–482. Li, X.Y., Xiao, N., Zhang, Y.H., 2014. Toxic effects of octylphenol on the expression of genes in liver identified by suppression subtractive hybridization of Rana chensinensis. Ecotoxicol 23 (1), 1–10. Li, X.Y., Liu, J., Zhang, Y.H., 2016. Octylphenol induced gene expression in testes of Frog, Rana chensinensis. Ecotoxicol. Environ. Saf. 128, 75–82. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using realtime quantitative PCR and the 2(-Delta Delta C(T)) method. Methods 25 (4), 402–408. Lundberg, R., Jenssen, B.M., Leiva-Presa, A., Rönn, M., Hernhag, C., Wejheden, C., Larsson, S., Orberg, J., Lind, P.M., 2007. Effects of short-term exposure to the DDT metabolite p,p′-DDE on bone tissue in male common frog (Rana temporaria). J. Toxicol. Environ. Heal. A. 70, 614–619. Maier, K.G., Roman, R.J., 2001. Cytochrome P450 metabolites of arachidonic acid in the control of renal function. Curr. Opin. Nephrol. Hypertens. 10 (1), 81–87. Mataqueiro, M.I., Nakaghi, L.S., Souza, J.P., Cruz, C.D., Oliveira, G.H., Urbinati, E.C., 2009. Histopathological changes in the gill, liver and kidney of pacu (Piaractus mesopotamicus, Holmberg, 1887) exposed to various concentrations of trichlorfon. J. Appl. Ichthyol. 25 (1), 124–127. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5 (7), 621–628. Nebert, D.W., Dalton, T.P., Okey, A.B., Gonzalez, F.J., 2004. Role of aryl hydrocarbon receptor-mediated induction of the CYP1 enzymes in environmental toxicity and cancer. J. Biol. Chem. 279 (23), 23847–23850. Palmer, R.M., Ferrige, A.G., Moncada, S., 1987. Nitric oxide release accounts for the biological activity of endothelium-derived relaxing factor. Nature 327 (6122), 524–526. Pruitt, K.D., Tatusova, T., Maglott, D.R., 2007. NCBI reference sequences (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucl. Acids Res. 35, D61–D65. Qiao, L., Yang, W., Fu, J., Song, Z., 2013. Transcriptome profile of the green odorous frog (Odorrana margaretae). PLoS One 8 (9), e75211. Ranson, H., Claudianos, C., Ortelli, F., Abgrall, C., Hemingway, J., Sharakhova, M.V., Unger, M.F., Collins, F.H., Feyereisen, R., 2002. Evolution of supergene families associated with insecticide resistance. Science 298 (5591), 179–181. Rysz, M., Bromek, E., Daniel, W.A., 2015. Activation of the brain serotonergic system decreases the expression and activity of liver cytochrome P450 isoforms CYP1A, 2C11 and 3A. Eur. Neuropsychopharm 25 (S2), 209–210. Soumis, N., Lucotte, M., Sampaio, D., Almeida, D.C., Giroux, D., Morais, S., Pichet, P., 2003. Presence of organophosphate insecticides in fish of the Amazon River. Acta Amazon. 33, 325–337. Thomaz, J.M., Martins, N.D., Monteiro, D.A., Rantin, F.T., Kalinin, A.L., 2009. Cardiorespiratory function and oxidative stress biomarkers in Nile tilapia exposed to the organophosphate insecticide trichlorfon (NEGUVON). Ecotoxicol. Environ. Saf. 72 (5), 1413–1424. Venturini, F.P., Moraes, F.D., Cortella, L.R., Rossi, P.A., Cruz, C., Moraes, G., 2015. Metabolic effects of trichlorfon (Masoten®) on the neotropical freshwater fish pacu (Piaractus mesopotamicus). Fish. Physiol. Biochem. 41 (1), 299–309. Wang, B.Y., Li, Q.X., Li, J., Xie, X.F., Ao, Y., Ai, Y.X., 2009a. Hepatotoxicity and gene expression down-regulation of CYP isozymes caused by renal ischemia/reperfusion in the rat. Exp. Toxicol. Pathol. 61 (2), 169–176. Wang, Z., Gerstein, M., Snyder, M., 2009b. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10 (1), 57–63. World Health Organization (WHO), 1992. Environmental Health Criteria 132: Trichlorfon, International Program on Chemical Safety. WHO, Geneva. Xu, W., Liu, W., Lu, K., Jiang, Y., Li, G., 2012. Effect of trichlorfon on oxidative stress and hepatocyte apoptosis of Carassius auratus gibelio in vivo. Fish. Physiol. Biochem. 38 (3), 769–775. Yang, Y., Chen, S., Wu, S., Yue, L., Wu, Y., 2006. Constitutive overexpression of multiple cytochrome P450 genes associated with pyrethroid resistance in Helicoverpa armigera. J. Econ. Entomol. 99 (5), 1784–1789. Zhang, M., Li, Y., Yao, B., Sun, M., Wang, Z., Zhao, Y., 2013. Transcriptome sequencing and de novo analysis for Oviductus Ranae of Rana chensinensis using illumina RNA-Seq technology. J. Genet. Genom. 40 (3), 137–140.
Acknowledgements This study was supported by the CAS “Light of West China” visiting scholar training project and the Foundation of the Shaanxi Academy of Science (2017k-05, 2017k-27). Appendix A. Supplementary material Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ecoenv.2017.09.016. References Audic, S., Claverie, J.M., 1997. The significance of digital gene expression profiles. Genome Res. 7 (10), 986–995. Boeckmann, B., Bairoch, A., Apweiler, R., Blatter, M.C., Estreicher, A., Gasteiger, E., Martin, M.J., Michoud, K., O'Donovan, C., Phan, I., Pilbout, S., Schneider, M., 2003. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucl. Acids Res. 31 (1), 365–370. Chang, C.C., Rahmawaty, A., Chang, Z.W., 2013. Molecular and immunological responses of the giant freshwater prawn, Macrobrachium rosenbergii, to the organophosphorus insecticide, trichlorfon. Aquat. Toxicol. 130–131, 18–26. Coelho, S., Oliveira, R., Pereira, S., Musso, C., Domingues, I., Bhujel, R.C., Soares, A.M., Nogueira, A.J., 2011. Assessing lethal and sub-lethal effects of trichlorfon on different trophic levels. Aquat. Toxicol. 103 (3–4), 191–198. Contreras-Vergara, C.A., Harris-Valle, C., Sotelo-Mundo, R.R., Yepiz-Plascencia, G., 2004. A mu-class glutathione S-transferase from the marine shrimp Litopenaeus vannamei: molecular cloning and active-site structural modeling. J. Biochem. Mol. Toxicol. 18 (5), 245–552. Daborn, P., Boundy, S., Yen, J., Pittendrigh, B., ffrench-Constant, R., 2001. DDT resistance in Drosophila correlates with Cyp6g1 over-expression and confers cross-resistance to the neonicotinoid imidacloprid. Mol. Genet. Genom. 266 (4), 556–563. Ekblom, R., Galindo, J., 2011. Applications of next generation sequencing in molecular ecology of non-model organisms. Heredity 107 (1), 1–15. Götz, S., García-Gómez, J.M., Terol, J., Williams, T.D., Nagaraj, S.H., Nueda, M.J., Robles, M., Talón, M., Dopazo, J., Conesa, A., 2008. High-throughput functional annotation and data mining with the Blast2GO suite. Nucl. Acids Res. 36 (10), 3420–3435. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29 (7), 644–652. Heeba, G.H., El-Hanafy, A.A., 2012. Nebivolol regulates eNOS and iNOS expressions and alleviates oxidative stress in cerebral ischemia/reperfusion injury in rats. Life Sci. 90 (11–12), 388–395. Herwig, N., 1979. Handbook of drugs and chemicals used in the treatment of fish diseases. In: Charles, C. (Ed.), A Manual of Fish Pharmacology and Material Medica. Springfield, Illinois, USA. Jacobs, M.N., Nolan, G.T., Hood, S.R., 2005. Lignans, bacteriocides and organochlorine compounds activate the human pregnane X receptor (PXR). Toxicol. Appl. Pharm. 209 (2), 123–133. Kanehisa, M., Goto, S., 2000. KEGG: kyoto encyclopedia of genes and genomes. Nucl. Acids Res. 28 (1), 27–30. Kawai, Y.K., Watanabe, K.P., Ishii, A., Ohnuma, A., Sawa, H., Ikenaka, Y., Ishizuka, M., 2013. De novo sequence analysis of cytochrome P450 1-3 genes expressed in ostrich liver with highest expression of CYP2G19. Comp. Biochem. Phys. D 8 (3), 201–208. Komagata, O., Kasai, S., Tomita, T., 2010. Overexpression of cytochrome P450 genes in pyrethroid-resistant Culex quinquefasciatus. Insect Biochem. Mol. 40 (2), 146–152. Koonin, E.V., Fedorova, N.D., Jackson, J.D., Jacobs, A.R., Krylov, D.M., Makarova, K.S., Mazumder, R., Mekhedov, S.L., Nikolskaya, A.N., Rao, B.S., Rogozin, I.B., Smirnov, S., Sorokin, A.V., Sverdlov, A.V., Vasudevan, S., Wolf, Y.I., Yin, J.J., Natale, D.A., 2004. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol. 5 (2), R7. Ku, P., Wu, X., Nie, X., Ou, R., Wang, L., Su, T., Li, Y., 2014. Effects of triclosan on the
493