Abstract Cypermethrin (CYP), a persistent synthetic pyrethroid, poses a significant threat to aquatic life. The present work aims to identify the effect of CYP on the brain and liver transcriptome profile of fish Labeo catla. CYP (0.7 µg/L) was long-term exposed to L catla, and differentially expressed genes (DEGs) were investigated using RNA-sequencing (Illumina HiSeq2500). A total of 2665 and 18,020 unigenes and 333 and 454 DEGs were identified in the brain and liver transcriptome, respectively. DEGs associated with MAPK signaling and apoptosis pathways were co-expressed in both brain and liver transcriptome, according to pathway enrichment analysis. Concurrently, steroid and terpenoid backbone biosynthesis were overexpressed in liver, suggesting induction of apoptosis and steroid production, respectively. Among the identified DEGs, fosab, nr4a1, dhcr24 and tm7sf2 were significantly enriched and molecular docking studies further supported our findings. Long-term CYP exposure also altered the levels of antioxidant enzymes superoxide dismutase (SOD), catalase (CAT), peroxidase (POD), and malondialdehyde (MDA) content in L. catla. The present study provides significant insights into the consequences of CYP-induced toxicity in fish and identified functional genes that could be screened for estimating the toxic load of CYP in future ecotoxicological risk assessment studies. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-86513-x. Keywords: Cypermethrin; L. catla; Brain, liver; Apoptosis; Steroid biosynthesis Subject terms: Next-generation sequencing, Bioinformatics, Gene expression analysis __________________________________________________________________ A sustainable increase in food production in parallel to the growing world’s population has been possible due to the use of pesticides in the agricultural sector, which has become an indispensable part contributing to a remarkable increase in agricultural production in the 20th century^[38]1. In the 1950s, the annual production of pesticides was 0.2 million tonnes, which has increased to more than five million tonnes by 2000; the production of pesticides increased globally at a pace of around 11% each year^[39]2. Every year, 3 billion kilogrammes of pesticides are used globally, but only 1% are employed to control insect pests on targeted plants effectively^[40]3,[41]4. The vast majority of pesticides still present are absorbed by or reach unintended plants and animals. Pesticide contamination has subsequently produced environmental pollution and adverse effects on human well-being^[42]5. Cypermethrin (CYP), a synthetic pyrethroid (Type II), is used against diverse anthropod classes and is utilised for controlling insect pests of various agricultural crops and also for mosquito, cockroaches and lice eradication^[43]6,[44]7. Additionally, CYP is used to check pests in the processing of wool, cotton, sheep dipping, soybeans and moths, the management of salmonids by controlling sea lice, in forestry, etc^[45]8–[46]11. Because of this, farmers became interested in CYP’s wide range of activities, leading to increased global usage^[47]12. CYP is a significant agricultural product due to its decreased toxicity to mammals and birds and minimal soil accumulation. However, rampant use of this pyrethroid has led to non-target creatures, including humans and animals exposed to it. As a result, there is a heightened risk of accumulation in the non-target animals through contaminated water^[48]13. CYP concentration in surface water has been observed to range within 0.01 to 9.8 ug/L, with farming runoff (riparian drainage canals of the river Ardas and Erythropotamos) reporting the highest value at about 194 µg/L^[49]13–[50]15. Among the non-targeted creatures, fish are more susceptible to pyrethroids because they lack carboxylesterase, an enzyme that breaks down pyrethroids; thus, can retain pyrethroids for longer periods than other animals^[51]16. It is crucial to conduct research on how different fish species react to pesticides for several reasons. Primarily, it might consider fish as a food commodity as well as economic activity. In addition, given that the majority of the pesticides can bioaccumulate as well as biomagnify along food chains^[52]17and fish enjoy higher trophic positions, this fact may not only showcase the negative effects of pesticides and the reactions of non-target fish species with economic relevance but can provide conclusive idea about the condition of the entire trophic system^[53]18. Labeo catla(Catla), an Indian major carp, is one of the most favoured fish because of its rapid growth and delicious meat. It makes up a significant portion of the composite fish culture production in India and the polycultures of Nepal, Bangladesh, Pakistan, and Myanmar^[54]19. Additionally, it helps many fish farmers in the southern part of Asia, and the Indian subcontinent support their families. It is one of the top twenty fish species produced around the globe, accounting for 84.2% of all finfish produced worldwide^[55]20. Recent reports showed that CYP could inflict neurotoxicity^[56]21, hepatotoxicity and apoptosis in L. catla fingerlings even at sub-lethal doses i.e. 0.12 and 0.41 µg/L for a duration of 45 days^[57]22. Our recent investigation revealed that the 96 h LC[50] concentration of CYP to be 7 µg/L in L. catla which is relatively a higher concentration than previous reports^[58]23. However, the molecular mechanism behind the induced neurotoxicity and hepatotoxicity is still unexplored. Due to its lower cost, high throughput capabilities, and comparability between species and pollutant impacts, transcriptome profiling has become more widely used in the last ten years better to comprehend the underlying effect of toxins in aquatic system. Using downstream bioinformatic functional annotation software in conjunction with RNA sequencing, one may relate changes in gene expression to diseases, projected ontological changes, and related modes of action responsible for these reactions^[59]24. The toxic effect of CYP using a transcriptomic platform has been reported earlier^[60]7,[61]25–[62]27. However, the effect of chronic exposure to CYP on the two most crucial organs, i.e. brain and liver used for toxicological assessment of any toxicant as well as the mechanism behind the induced toxicity, remained unaddressed. Recent studies showed that CYP exposed female mice at human acceptable daily intake or chronic reference dose (60 µg/kg/day) for 44 weeks induced primary ovarian insufficiency by up-regulating pro-apoptotic protein Bcl-2 (B cell lymphoma/leukemia type 2) modifying factor as well as cell cycle inhibitor p27 as revealed from the ovarian transcriptome^[63]28. Similarly, after 15 days of postnatal exposure to CYP (5 and 20 mg/kg body weight), a brain transcriptomic study of mice exhibited developmental neurotoxicity by affecting pathways related to mitochondrial dysfunction, biogenesis, protein folding, trafficking and degradation^[64]29. However, in mice liver transcriptomic study, post-liver injury induced by beta-CYP, the significantly enriched pathways were retinol and linoleic acid metabolism and the JaK-STAT signaling pathway^[65]26. Exposure to an environmentally relevant dose (0.3 ng/bee) of CYP to Apis mellifera (honeybees) showed overexpression of pathways regulating muscular and cellular processes and metabolic enzymes^[66]25. The present study intends to fill the knowledge gap and comprehend the mechanism of CYP-induced toxicity in teleost, reorientation of critical genes and related pathways using transcriptomic analysis post long-term CYP exposure to L. catla followed by validation of selective genes by qPCR. Given that the brain and liver are the two vital organs and indicators of their physiological status, these organs were selected for comparative transcriptomic study. A clear understanding of the toxicological mechanism of CYP to L. catla may also identify novel genes and/or pathways that have implications in higher vertebrates under similar toxic conditions and also for environmental risk assessment studies. Moreover, this approach could also help identify adverse outcome pathways (AOPs) for future studies^[67]30 to understand the potential interaction of CYP with the brain and liver of fish, the two main organs for toxicological studies. Materials and methods Ethics statement The present study was done per approval from the Institute Animal Ethics Committee (IAEC) (ICAR-CIFRI/IAEC-21–22/01). The experimental design, sample collection and fish’s sacrifice were done following the ethical standards and were in accordance with the relevant guidelines and regulation of the study country and also all methods reported in the present study were in accordance with ARRIVE guidelines. Animals and chemical L. catla (length: 14 ± 2.56 cm; weight: 16 ± 2.65 g) were procured from a fish hatchery at Chamta, West Bengal, India and maintained in the hatchery of ICAR-CIFRI following preventive measures. The fish were fed twice on a daily basis (3% of body weight) with commercial feed (moisture: 10%, crude fat: 3%, crude protein: 28%, and crude fiber: 4%) and acclimatised for two weeks at 25.0 ± 1 °C in dechlorinated water (pH 7.5 ± 0.2), ensuring proper aeration and 12 h alternate day and night cycle before the exposure study. The leftover feed and faecal waste were removed daily. Cypermethrin (CAS No. 52315-07-8, > 97% purity) was purchased from M/S Sigma-Aldrich. A stock solution (100 mg/L) was prepared in hexane and further diluted according to the required dosage from the stock. All the solvents were first purified and redistilled before preparing stock solutions. Water quality characteristics Water quality characteristics (pH, conductance, temperature, salinity, total dissolved solids (TDS), turbidity) were analysed using a Multiparameter water quality probe (YSI ProDSS), and hardness was estimated by titration^[68]31 on a weekly basis. Experimental design For the transcriptomics study, 0.7 µg/L (1/10th of LC[50]) of CYP was selected according to the 96 h LC[50]value obtained from our previous study^[69]23. Previously, 96 h LC[50] value was determined using CYP concentration that lies between the safe and lethal dose i.e. 1, 3, 5, 7, 10 µg/L) for L. catla^[70]23. L. catla fingerlings (n = 60) were randomly allocated (ten fish per tank) in rectangular tanks, one experimental and one control group, each group having three replicates. The control group consisted of only dechlorinated water with no toxicant, and the experimental group was exposed to 0.7 µg/L CYP in a semi-static renewal mode where the given toxicant at the specified concentration is replenished everyday during the exposure period at a fixed time (12 noon) after addition of half fresh solution. Fish were fed twice daily at a particular time ( 1.00 pm) with commercial feed, and feeding was stopped two days before the 30th day of exposure when the fish were sacrificed. During the period of exposure, no mortality was observed. The fishes were anaesthetized using Tricaine (Sigma-Aldrich) after completion of the exposure period and tissues were collected accordingly. For antioxidant enzyme assay, brain and liver sections were kept in 2% sucrose solution until further processing in −40 ºC. Two fishes were taken from each tank of control and experimental group, respectively and a total of six individual fish samples (n = 6) were analysed. For transcriptomic study, whole brain and liver tissues were carefully taken out from individual fishes from control and experimental group. Three replicate were prepared for each group consisting of pooled tissues from three individual fishes. Therefore, each replicate consist of pooled samples from three individual fishes. The replicates were then kept in RNALater (Sigma), and stored at −40 ºC until further processing. Antioxidant enzyme assay The antioxidant enzymes i.e. superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD) activities as well as malondialdehyde (MDA) content was evaluated in brain and liver tisues. The total protein concentrations were measured in the brain and liver tissues following the method of Lowry et al.^[71]32. The SOD, CAT and POD activity were assessed following methodologies stated earlier^[72]33–[73]35 and MDA content were measured using commercially available Fish MDA ELISA Kit (MyBioSource, San Diego, USA) using the kit instructions. The MDA assay detection limit was 5.39 pg/mL.The enzyme activities were expressed in terms of units/milligrams of protein. RNA isolation, cDNA library preparation and sequencing Extraction of total RNA was done from brain and liver tissues using RNeasy Micro Kit (Qiagen, Germany) per the stated guidelines. The extracted RNA’s concentration as well as purity were evaluated using a Nanodrop Spectrophotometer (Thermo Scientific, 2000), and the RNA integrity was calculated with the help of Agilent Bioanalyzer 2100 (Agilent, CA, USA) using the manufacturer’s instructions. RNA integrity number (RIN) value greater then 7.5 were acceptable for further processing^[74]36. Illumina compatible NEBNext Ultra Directional RNA Library Preparation Kit (New England BioLabs, MA, USA) was utilised for RNA sequencing libraries preparation. For removal of any DNA contamination, DNase I treatment were given, and to further purify it, magnetic beads along with oligo(dT) (Invitrogen, USA) were utilised. To break the mRNA into short fragments of approximately 300 bp, a fragmentation buffer was added. cDNA synthesis was accomplished with short mRNA fragments and random hexamers (first strand), and the cDNA second strand was eventually synthesized by using DNA pol I and RNase H. The amplification of adaptor-ligated, size-selected cDNA was done using PCR, and to construct the final cDNA library, it was isolated on 2% agarose gels. Bioanalyzer (Agilent 2100) and qPCR (ABI StepOnePlus System) were used for qualitative and quantitative evaluation of the sample library. The library products were then sequenced using Illumina HiSeq2500 in a 2 × 150 bp paired-end format to generate 40–50 M reads per sample. Image data output was then transformed to raw reads by base calling and eventually stored in fastq format. The quality of the generated raw reads from sequencing was determined by fastp 0.23.2 program. After removing the adapters and low-quality reads, only the high-quality reads were utilised for further analysis^[75]37. For alignment of clean data with the reference genome, the BWA mem version 0.7.17 program was used^[76]38. The genome of Danio rerio (zebrafish, GRCz11) was utilised for mapping. The union mode in HTSeq version 0.9.1 was utilized to evaluate the gene expression level in this study^[77]39. The sequence data for the brain and liver generated in this study is submitted to NCBI, Sequence Read Archive (SRA), PRJNA981973 and PRJNA992386, respectively. Differential expression analysis and functional annotation False discovery rate (FDR) was minimised for the identification of differentially expressed genes (DEGs)^[78]40. DEGs were statistically analysed using the DESeq2 version 1.34.0^[79]41 with an adjusted P value 0.05 and an absolute Log2fold change value 1 or −1. For enrichment-based pathway analysis, the DEGs were mapped according to the terms of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using DAVID (Database for Annotation, Visualization and Integrated Discovery). Validation of selected DEGs by qPCR To validate the transcriptomic data, expression of the selected genes from the identified DEGs was analyzed by qPCR and gene-specific primers were designed employing Primer 3 software using Danio rerio as the reference genome. Six selected genes from DEGs associated with MAPK signaling, apoptosis and steroid biosynthesis pathway were used to validate brain and liver transcriptomics data and subjected to qPCR analysis. The genes that were considered for qPCR validation in the present study were proto-oncogene protein c-fos (fosab), nuclear receptor subfamily 4 group A member 1 (nr4a1), growth arrest and DNA-damage-inducible protein (gadd45ba), metastasis suppressor KiSS-1 (kiss1), delta 24-sterol reductase (dhcr24), and delta 14-sterol reductase (tm7sf2) after CYP exposure (Table [80]S1). The total RNA from brain and liver tissue was extracted using TRI Reagent (Sigma-Aldrich) as per the manufacturer. The quantification of the extracted RNA and its purity were determined by spectrophotometer (Multiskan SkyHigh Microplate Spectrophotometer, Thermo Scientific) and Bioanalyzer (2100, Agilent, Santa Clara, CA, USA), respectively. RNA samples with absorbance between 1.5 and 2 at 260/280 nm were selected for qPCR analysis. The quality of the RNA was rechecked by the presence of two clear bands, i.e. 28 S and 18 S, on agarose gel electrophoresis (1.5%). Then, the cDNA synthesis was carried out by using a cDNA synthesis kit by Qiagen (QuantiNova™ Reverse Transcription kit). Before the qPCR analysis, the primer efficiency was checked by running serial dilutions of the primers designed, and then constructing standard curves from the generated data and the primer efficiency ranged within 95–100%. CFX Connect (Biorad, UK) real-time PCR was employed for qPCR analysis, and each sample’s reaction volume was 20 µL. A particular reaction was composed of iTaq Universal SYBR Green Supermix (BioRad) (10 µL), 10 pM of each primer for the selected gene (1 µL), 10 times diluted cDNA (5 µL), and RNase/DNase free water (3 µL). The conditions for amplification used for qPCR were: 3 min of denaturation (95 °C), again 15 s denaturation (95 °C) for 40 cycles, followed by annealing and extension at 60 ºC for 30 s. The relative-fold change was calculated with the help of 2^−ΔΔCtmethod^[81]42. β-actinwas chosen as the housekeeping gene in this investigation and its stability was previously assessed using the comparative delta Ct method^[82]23. Protein and ligand structure prediction and preparation The genes identified as potential biomarkers in the present study for assessing CYP-induced toxicity were further selected for molecular docking studies. The complete target protein sequences of fosab, nr4a1, dhcr24, tm7sf2 genes of Danio rerio i.e. v-fos (NP991132), Nr4a1([83]AAI64696), Dhcr24 ([84]AAI65211), TM7SF2 (NP001008597) were retrieved from the NCBI database in FASTA format ([85]https://www.ncbi.nlm.nih.gov/) as Danio rerio was used as a reference genome in the present study. The templates suitable for homology modelling were searched for similar sequences against protein data bank (PDB) using BlastP and the PDB ID: 2wt7, 3ej8, 4quv, 7wnh with a resolution of 2.30 Å, 2.55 Å, 2.74 Å and 3.10 Å were selected for v-fos, Nr4a1, Dhcr24 and TM7SF2 respectively as the best templates for 3D modelling having a high percentage of sequence identity and query coverage ([86]http://www.rcsb.org/). Using CLUSTAL X, alignment was done for these query protein sequences with 3D modelling templates (PDB: 2wt7, 3ej8, 4quv, 7wnh), and crystal structures for the query sequence were identified. Homology modelling was needed to construct the structure of the target proteins, and MODELER10.3 package with a minimum model score, the lowest discrete optimised protein energy score ([87]https://salilab.org/modeller/download_installation.html) was used for predicting the 3-D structures of the templates (PDB ID: 2wt7, 3ej8, 4quv, 7wnh). Galaxy Loop carried out the structural refinement process and side chain optimisation of the modelled structure, and the loop-refined 3D model was refined subsequently using Galaxy Refine (https://galaxy.seoklab.org/cgi-bin/submit.cgi? type=REFINE. The optimised model was subjected to energy minimisation using Swiss-Pdbviewer_4.10 ([88]https://spdbv.unil.ch/). Three minimised models were validated using PROCHECK ([89]https://saves.mbi.ucla.edu/) considering amino acids both in the allowed and disallowed region as well as overall G-factor. The structure of CYP available in the PubChem database ([90]https://pubchem.ncbi.nlm.nih.gov/) in SDF format was taken and then changed to PDB format by SMILES Translator ([91]https://www2.chemie.uni-erlangen.de/services/translate/). Then, it was subjected to the Avogadro platform ([92]https://avogadro.cc/) for energy minimization. On these compounds, an 500 steps steepest descent algorithm along with MMFF94 force field was applied. Insilico molecular docking studies Autodock 4.2.6 software was used for ligand-protein structure-based docking for v-fos, Nr4a1, Dhcr24 and TM7SF2 proteins with ligand i.e. CYP^[93]43. Before the in silico analysis process, protein structure and ligand were prepared after removing water molecules, assigning Kollman charges, adding partial charges, and hydrogen atoms with polar contacts. The target protein’s active sites were identified through CASTp 3.0 software^[94]44. Electrostatic grid maps were estimated using Auto Grid version 4.0 to determine the protein’s binding site. To potentiate the binding of CYP (ligand) to the target proteins, grid maps with 0.375 Å, 0.375 Å, 0.375 Å, 0.564 Å spacing were generated for v-fos, Nr4a1, Dhcr24 and TM7SF2, respectively. The X, Y, and Z dimensions and the X, Y, and Z centre were generated that cover the entire macromolecule. For v-fos, the grid size of 114 × 124 × 126 xyz points, with center coordinates of 65.516, −80.958, 9.227 (X, Y, Z), Nr4a1 (120 × 122 × 124 xyz points, grid center of −10.721, −0.522, −25.862 (X, Y, Z), Dhcr24 (126 × 124 × 126 xyz points, grid center of −34.113, −6.756, 14.070 (X, Y, Z) and TM7SF2 126 × 126 × 126 xyz points, grid center of −20.736. 26.358, 23.379 (X, Y, Z) were selected for molecular docking. For best conformers, the Lamarckian Genetic Algorithm (LGA) was used. In total, 100 LGA runs were carried out to identify the target proteins’ binding site, and the ligand’s docking conformation with the lowest binding energy was chosen for further analysis. The docked complexes were visualised, and the receptor-ligand interactions like hydrogen bonds, hydrophobic interactions, etc., were analysed with the help of Discovery Studio Visualize software ([95]https://discover.3ds.com/discovery-studio-visualizer-download). Further, for a better understanding the stability of the target proteins (v-fos, Nr4a1, Dhcr24, TM7SF2) with CYP, a 100ns MD simulation was performed for each of the CYP-protein complexes and their details is provided as supplementary material (S1). Statistical analysis The results of antioxidant enzyme assay and qPCR analysis are shown as mean ± standard error. The data were analysed using Student’s t-test using SPSS (version 16.0). Statistical significance was measured at P < 0.05. DAVID 2021 ([96]https://david.ncifcrf.gov/), gene enrichment tool was used for GO and KEGG pathway analysis and the statistical parameters applied were Benjamini multiple test correction^[97]45,[98]46. The identified significant DEGs between the control and CYP-treated group were also represented via a Venn diagram using Venny 2.1 ([99]https://bioinfogp.cnb.csic.es/tools/venny/). Results Water quality parameters The water quality characteristics in the present study were found to be in the following range: pH 7.85–8.10, temperature 26.1–27.9 ºC, conductance 340–460 µS cm^−1, salinity 185–220 mg L^−1, TDS 245–325 mg L^−1, turbidity 0.35–2.10 NTU, hardness 145–155 mg L^−1. No significant difference in water quality characteristics was observed in CYP-treated as well as control samples. Antioxidant response on CYP exposure The antioxidant enzymes SOD, CAT and POD activity (Fig. [100]1A) and MDA content (Fig. [101]1B) was found to increase significantly (P < 0.05) in the exposed group (0.7 µg/L) in respect to control in both brain and liver tissues of L. catla on chronic exposure to CYP. Fig. 1. [102]Fig. 1 [103]Open in a new tab Comparative antioxidant profiles of SOD, CAT, and POD activity (A) and MDA content (B) in brain and liver tissues of L. catla on exposure to CYP in respect to control. Values are depicted as mean ± SE (n = 6). Different letters are indicative of significant differences at P < 0.05. Transcriptome analysis of brain and liver tissues and identification of DEGs Transcriptome sequencing of libraries from L. catla ‘s brain and liver tissues on experimental exposure to CYP with respect to individual control in triplicate was done via Illumina HiSeq2500. On removing short reads, low-quality reads and reads related to mitochondria, from the total reads, clean reads were obtained for individual samples (Table S2). DESeq2 program identified a total of 26665 genes, which showed differential expression in brain samples of L. catla on experimental exposure to CYP in respect to control (Fig. [104]2A). A total of 132 genes among the total DEGs were up-regulated significantly, and 201 genes were down-regulated in brain samples of L. catla in experimental group exposed to CYP in respect to control (Fig. [105]3). Similarly, 18020 unique genes showed differential expression in liver samples of L. catla in experimental group exposed to CYP in respect to control (Fig. [106]2B) among which 186 genes were significantly up-regulated and 268 genes were down-regulated (Fig. [107]3). Fig. 2. [108]Fig. 2 [109]Open in a new tab Volcano plot of brain (A) and liver (B) samples of L. catla on exposure to CYP in respect to control. The volcano plot is generated based on the p-values. The p-values < 0.5 and Log2Fold Change > 1 were considered up-regulated, and p-values < 0.5 and Log2Fold Change <−1 as down-regulated. Other genes were considered as neutral. R package ‘ggplot2’ has been used for generating the plots. Fig. 3. [110]Fig. 3 [111]Open in a new tab Venn-diagram of up-regulated and down-regulated genes in brain and liver transcriptome of L. catla on exposure to CYP. Functional annotation of DEGs Functional annotation of significant DEGs in brain and liver transcriptome of L. catla on exposure to CYP were done in terms of GO, including cellular component (CC), molecular function (MF) and biological processes (BP) and are depicted in Fig. [112]4. GO annotation of the significant DEGs identified in the brain transcriptome revealed the involvement of maximum number of genes in the intracellular anatomical structure followed by organelle and intracellular organelle under the cellular component. In terms of biological process, the highest number of genes were involved in biological regulation, followed by regulation of biological process and regulation of cellular process. Under the molecular function category, most of the genes were associated with binding, followed by organic cyclic compound binding and heterocyclic compound binding (Fig. [113]4A). Fig. 4. [114]Fig. 4 [115]Open in a new tab GO annotation of significant DEGs in terms of cellular component, molecular function and biological process of L. catla ‘s brain (A.) and liver (B.) tissues exposed to CYP. GO annotation of the significant DEGs identified in the liver transcriptome revealed the involvement of a maximum number of genes in mitochondrion followed by membrane-enclosed lumen and organelle lumen under the cellular component. In terms of molecular function, the maximum number of genes were involved in oxidoreductase activity, followed by transition metal ion binding and hydrolase activity, acting on ester bonds. In the biological process category, the highest number of genes were involved in the small molecule metabolic process, followed by the catabolic process and lipid metabolic process (Fig. [116]4B). Pathway enrichment analysis Pathway enrichment analysis of the identified DEGs in L. catla ‘s brain and liver transcriptome on exposure to CYP was done in terms of GO and KEGG to identify biological pathways and gene functions that were substantially more expressed in each tissue comparison. According to GO terms of identified DEGs in brain transcriptome, enrichment analysis revealed myofibril, contractile fiber and U1snRNP as the top three subcategories under the cellular component category. In the molecular function category, glucocorticoid receptor binding, steroid hormone receptor binding and RNA polymerase II regulatory region sequence-specific DNA binding were dominant. Under the biological process category, pyrimidine dimer repair, copper ion transport, and excitatory postsynaptic potential were most significant (P < 0.05) (Supplementary Fig. 1A). Similarly, enrichment analysis identified DEGs in liver transcriptome revealed rRNA methyltransferase activity, enoyl-CoA hydratase activity, and 3-hydroxy acyl-CoA dehydrogenase activity were most dominant under the molecular function category. In the category belonging to biological processes were cholesterol biosynthetic, sterol biosynthetic, and sterol metabolic processes were most significant (P < 0.05). Under the cellular component, cytosolic proteasome complex, mitochondrial respiratory chain complex III and small nucleolar ribonucleoprotein complex were highly expressed (Supplementary Fig. 1A and 1B). KEGG pathway enrichment analysis of the identified DEGs in the brain transcriptome showed that the DEGs that were up-regulated were mainly involved in the MAPK signaling pathway (ko04010) and p53 signaling pathway (ko04115), and the down-regulated genes were associated with neuroactive ligand-receptor interaction (ko04080) and p53 signaling pathway (ko04115). Similarly, the up-regulated DEGs identified in the liver transcriptome were mainly related to metabolic pathways (ko01100), steroid biosynthesis (ko00100) and MAPK signaling pathway (ko04010) and the down-regulated DEGs in the liver transcriptome of L. catla were mainly involved with metabolic pathways (ko01100), peroxisome proliferator-activated receptors (PPAR) signaling pathway (ko03320) and fatty acid metabolism (ko01212) (Table [117]1)^[118]47. Table 1. KEGG pathways (top fifteen) of DEGs identified in the brain and liver transcriptome of L. catlaon exposure to CYP^[119]47. Sl no. Pathways Pathway ID Gene number Genes Brain 1. Neuroactive ligand-receptor interaction ko04080 8 chrnb1l, npy1r, gh1, p2rx4a, adra2db, rxfp3.3a1, gabrz, tacr1b 2. MAPK signaling pathway ko04010 7 nr4a1, gadd45ba, mras, fosab, dusp1, cdc4212, tradd 3. p53 signaling pathway ko04115 7 gadd45ba, zgc:92658, ddb2, apaf1, ccnb1, steap3, chek1 Liver 4. Metabolic pathways ko01100 83 gda, mvk, dhcr24, g6pca.2, rdh12, idi1, lpin1, csad, lipg, tm7sf2, ppcdc, gatm, mat2ab, st3gal5, tsta3, fdps, cyp51, sc5d, cyp27b1, ebp, pomgnt1, aldob, msmo1, aldoca, lss, gck, alox12, chac1, pmvk, aldh5a1, aldh1l1, pnp5a, sqlea, hsd17b7, mvda, glula, fdft1, hmgcs1, dhcr7, ogdhb, miox, impa1, dhfr, elovl6, impdh2, hprt1, cahz, adssl, comta, pla2g12b, acadvl, fads2, cyp2u1, hmox1a, impdh1b, cmpk2, scd, nt5c2l1, ndufa6, tmem86b, alas2, ndst3, ca4a, suclg2, ca12, cyp8b1, lipca, adob, hadhaa, si: dkey-91i10.3, arg1, hadhab, galnt11, cyp7a1, ckba, ldhbb, kl, zgc:77938, cox8b, cyp8b2, pnpla3, cyp24a1, pcyt1bb 5. Steroid biosynthesis ko00100 13 dhcr24, tm7sf2, cyp51, sc5d, cyp27b1, ebp, msmo1, lss, hsd17b7, fdft1, dhcr7, lipf, soat2 6. PPAR signaling pathway ko03320 13 apoa1a, fads2, scd, angptl4, cpt2, plin2, cyp8b1, si: dkey-91i10.3, fabp1b.1, cyp7a1, cyp8b2, lpl, slc27a1b 7. MAPK signaling pathway ko04010 8 nr4a1, gadd45bb, igf2a, gadd45ba, fosab, igf2b, jun, gadd45aa 8. Fatty acid metabolism ko01212 7 elovl6, acadvl, fads2, scd, cpt2, hadhaa, hadhab 9. Fatty acid degradation ko00071 6 acadvl, cyp2u1, cpt2, hadhaa, zgc:77938, eci2 10. Insulin signaling pathway ko04910 6 g6pca.2, socs3a, pik3r1, ppp1r3b, socs2, gck 11. Apoptosis ko04210 6 gadd45bb, gadd45ba, fosab, pik3r1, jun, gadd45aa 12. Lysosome ko04142 6 lipf, ctsk, tpp1, ctsd, nots, pla2g15 13. Toll-like receptor signaling pathway ko04620 6 fosab, pik3r1, jun, stat1a, ctsk, stat1b 14. Terpenoid backbone biosynthesis ko00900 5 mvk, idi1, fdps, mvda, hmgcs1 15. FoxO signaling pathway ko04068 5 gadd45bb, g6pca.2, gadd45ba, pik3r1, gadd45aa [120]Open in a new tab Identification of DEGs altered with CYP exposure The identified DEGs were further sorted according to fold change, p-value, a amalgamation of pathway enrichment analysis, annotation and manual literature searches to differentiate the most significantly affected pathways and the corresponding genes. The top five up-regulated and down-regulated genes in the brain and liver transcriptome of L. catla according to log2fold change value are depicted in Fig. [121]5. The DEGs fosab (K04379), gadd45ba (K04402) and nr4a1 (K04465) involved in MAPK signaling pathway (ko04010) and apoptosis pathway, ko04210 (gadd45ba and fosab) were significantly up-regulated in both brain and liver transcriptome. Similarly, the genes btg2 (protein BTG2, K23140) and sik1 (serine/threonine-protein kinase SIK1, K19008) were significantly up-regulated in both brain and liver transcriptome (Table S3) and the gene kiss1 (Kisspeptin 1(metastasis-suppressor KiSS-1, K23140)), up-regulated in the brain transcriptome plays a significant role in GnRH secretion (map04929). Among the down-regulated genes, npy1r (neuropeptide Y receptor type 1, K04465) and gh1 (growth hormone, K05438) were significantly altered and are associated with neuroactive ligand-receptor interaction (ko04080), and apaf1(apoptotic protease activating factor-1) is involved in p53 signaling pathway (ko04115). Fig. 5. [122]Fig. 5 [123]Open in a new tab Top five up-regulated and down-regulated genes identified in L. catla ‘s brain and liver transcriptome on CYP exposure. The green and ash color indicates the top five upregulated genes in the brain and liver transcriptome, respectively and the red and orange color indicates the top five down-regulated genes in the brain and liver transcriptome, respectively. Abbreviations: proto-oncogene protein c-fos (fosab); growth arrest and DNA-damage-inducible protein (gadd45ba); metastasis-suppressor KiSS-1 (kiss1); protein BTG2 (btg2); nuclear receptor subfamily 4 (nr4a1); calcium homeostasis (calhm3); FMS-like tyrosine kinase 4 (flt4); deoxycytidylate deaminase (dctd); neuropeptide Y receptor type 1 (npy1r); methylenetetrahydrofolate reductase (NADPH) (mthfr); cytochrome P450 family 1 subfamily A (cyp1a); serine/threonine-protein kinase SIK1 (sik1); Delta14-sterol reductase (tm7sf2); diphosphomevalonate decarboxylase (mvda); isopentenyl-diphosphate (idi1); arginase (arg1); apolipoprotein A-I (apoa1a); cholesterol 7alpha-monooxygenase (cyp7a1); perilipin-2 (plin2); elongation of very long chain fatty acids protein 6 (elovl6). The DEGs showing significant up-regulation in the liver transcriptome were mainly associated with steroid biosynthesis pathway (ko00100), including dhcr24 (K09828), tm7sf2 ([124]K00222) and cyp51 (sterol 14 alpha-demethylase, K05917) genes and mvk (mevalonate kinase, [125]K00869), mvda (diphosphomevalonate decarboxylase, [126]K01597) and idi1 (isopentenyl-diphosphate delta-isomerase, K01823) were involved in terpenoid backbone biosynthesis pathway (ko00900). Furthermore, the down-regulated DEGs were mainly involved in PPAR signaling pathway (ko03220) including apoa1a (apolipoprotein A-I, K08757), scd (stearoyl-CoA desaturase, [127]K00507), plin2 (perilipin-2, K17284), cyp8b1(cytochrome P450 8B3 (sterol 12-alpha-hydroxylase, K07431)), cyp7a1(cholesterol 7 alpha-monooxygenase, [128]K00489), cyp8b2 (cytochrome P450 8B2 (sterol 12-alpha-hydroxylase, K07431)) and fatty acid metabolism pathway (ko1212) including elovl6 (elongation of very long chain fatty acids protein 6, K10203), scd and hadhaa (enoyl-CoA hydratase / long-chain 3-hydroxyacyl-CoA dehydrogenase, K07515). Verification of RNASeq data using qPCR The qPCR results, i.e., fold change as well as direction of change (up or down-regulated), concur with the brain and liver transcriptomic data for the candidate genes selected for the present study (Supplementary Fig. 2). Molecular interactions of CYP with selected significantly enriched DEGs Pathway enrichment analysis showed that few DEGs were associated with multiple pathways. For example, fosab and nr4a1 are among the top up-regulated genes related to MAPK signaling and apoptosis pathways identified in both the brain and liver transcriptome. Similarly, in the liver transcriptome, the metabolic pathways followed by the steroid biosynthesis pathway was the top most enriched pathways in which dhcr24 and tm7sf2 genes are among the top up-regulated genes. Therefore, these genes were identified as potential biomarkers and were tested by molecular docking studies to look for potential binding and further validation. To confirm the interactions, if any, between the significantly enriched DEGs fosab, nr4a1, dhcr24 and tm7sf2 corresponding to proteins v-fos, Nr4a1, Dhcr24 and TM7SF2, respectively, with CYP as ligand, molecular docking studies were performed. Docking studies showed that the binding affinity of v-fos with CYP was − 6.14 Kcal/mol and formed one hydrogen bonding (LYS-257) accompanied by a hydrophobic contact (LEU-248) and a Pi-Lone Pair (Thr-247). The binding affinity of Nr4a1 and CYP was − 8.50 Kcal/mol. The amino acid residues of target protein Nr4a1 (PHE-211, LEU-412, PRO-154, ILE-186, ILE-232, ARG-230, TYR-180, HIS-185, LEU-416, PRO-421) interacted through hydrophobic bonds and one electrostatic interaction (GLU-204). Dhcr24 and CYP exhibited the highest affinity with a binding energy of −10.10 kcal/mol. LYS-314, ARG-318 and ASN-353 formed hydrogen bonds and TYR-354, TRP-346, MET-251, PHE-307, LYS-314, LEU-340 were associated with hydrophobic interactions. The binding energy of Tm7SF2 with CYP was − 7.96 Kcal/mol and formed three hydrogen bonds (THR-173, GLN-182 and PRO-192) and five hydrophobic interactions (ALA-174, VAL-172, PRO-192 and PRO-188) (Fig. [129]6). Fig. 6. [130]Fig. 6 [131]Open in a new tab Molecular docking results of CYP and target proteins (A) v-fos, (B) Nr4a1, (C) Dhcr24, (D) TM7SF2. The active pocket residue in the target proteins is in orange, and the ligand is in yellow. The green and pink dotted lines represent hydrogen bonding, hydrophobic interaction, and pipi conjugations. Discussion Transcriptome analysis identifies DEGs under challenging conditions, resulting in a clear comprehension of the genes and/or pathways connected to it and also gives a dynamic depiction of the biological processes occurring inside a living system^[132]36. In the present study, comparative brain and liver transcriptomics of L. catla after long-term exposure for 30 days to CYP, equivalent to 1/10th 96 h LC[50]dose (0.7 µg/L)^[133]23 was generated using Illumina HiSeq2500 which is the first study investigating the effect of chronic exposure of CYP on the commercially important food fish, L. catla brain and liver transcriptomic alterations. Mostly pesticides target pests primarily by acting on their nervous system; however, long-term exposure to CYP has been linked to a variety of adverse effects, including chronic and persistent neurotoxic outcomes, immunosuppression, teratogenic effects, increase in oxidative stress through greater production of reactive oxygen species (ROS), and alterations in the antioxidant enzyme system, which results in oxidative damage^[134]48, also observed in the present study with an increase in SOD, CAT, POD as well as MDA content on exposure to 0.7 µg/L concentration, similar to our previous study^[135]23. Since the liver and brain are primarily in charge of controlling the aforementioned effects of toxic exposure, transcriptome studies were performed to obtain a comprehensive picture of the impact of long-term exposure to CYP on L. catla^[136]24,[137]27. Pathway enrichment analysis revealed that the identified DEGs in the brain and liver transcriptome of L. catla following 30 days of CYP exposure were able to significantly affect pathways such as neuroactive ligand-receptor pathway, signaling pathways like MAPK, p53 and PPAR followed by steroid biosynthesis, fatty acid metabolism and degradation, terpenoid backbone biosynthesis as well as apoptosis pathway. Molecular docking studies revealed that all the selected proteins exhibited strong affinity towards CYP, further supporting our findings that these genes/proteins could serve as early markers to assess CYP-induced toxicity in fish even at chronic exposure to low concentrations for a long time. Together, these findings may elucidate the molecular mechanism of toxicity induced by CYP in L. catla; however, further validation is necessary to be more conclusive. Previously, Zhang et al.^[138]49. exposed zebrafish embryos to CYP for 96 h and reported extracellular matrix proteins (ECM)-receptor interaction and DNA replication as the most affected pathway on CYP exposure to zebrafish embryos. Similarly, Ranjani et al.^[139]7. reported CYP-induced transcriptomic changes involving disruption of the phototransduction pathway in zebrafish embryos at 48 hpf. Liver transcriptome analysis following CYP exposure at different concentrations (0.625 and 1.25 ng/L) to Cherax quadricarinatus, red claw crayfish revealed that the antioxidant, biotransformation modulation, mitochondrial dysfunction, inflammatory response, detoxification and autophagy impairment, and reproduction were mainly affected with long-term exposure^[140]27. Brain transcriptomics studies in Oncorhynchus mykiss (rainbow trout) following exposure to another environmental predominant pyrethroid, bifenthrin, for two weeks at low concentrations, affected neuroendocrine functioning, impaired extracellular matrix’s structural integrity, disrupted the brain’s cell signaling pathways and induced neuronal death^[141]24. The previous reports of the effect of CYP exposure on transcriptomic alterations in fish model are mostly on short-term exposure; this is perhaps the first study reporting the consequences of CYP exposure for long-term exposure on L. catla. MAPK signaling pathway (ko04010) was the only pathway that was up-regulated in both brain and liver transcriptome of L. catla, and the co-expressed DEGs associated with this pathway were fosab and junbb (part of activator protein-1 (AP-1) family of transcription factors), gadd45ba (growth arrest and DNA-damage-inducible protein), and nr4a1(a pro-apoptotic transcription factor associated with cellular apoptosis). MAPK signaling pathway mainly controls cellular reactions and signal transduction processes, and the three primary MAPKs are p38 MAPK, c-Jun N-terminal kinase (JNK), and extracellular signal-regulated kinase 1 and 2 (ERK1/2). They also control mitochondrial division and the mitochondrial apoptotic pathway and are critical players in cell proliferation and death^[142]50. Considering the DEGs associated with this pathway, significantly up-regulated in the brain transcriptome were fosab, mras, dusp1, gadd45ba, nr4a1, cdc4212 and tradd. dusp1 inactivates extracellular signal-regulated kinase (ERK) by dephosphorylation and is a pro-apoptotic gene promoting apoptosis in addition to its phosphatase activity and is reported to play a significant role in controlling inflammation and human monocytic cell line (THP-1) macrophage-mediated apoptosis in response to BCG infection^[143]51. fosab encoding Fosab protein, the zebrafish orthologue of c-Fos, a proto-oncogene, was found to accumulate in brain neurons following organophosphorus diisopropyl fluorophosphate exposure to zebrafish larvae due to neuronal hyperexcitation and over-expression of fosab indicates neuronal apoptosis^[144]52. nr4a1 is one of the three genes that encode the nuclear transcription factor NR4A, which induces apoptosis on translocation from the nucleus to mitochondria and is reported to be over-expressed on exposure studies (in vitro) on mononuclear cells from the bone marrow (BMMCs) when insecticide malathion was exposed at a concentration of 0.1 µM^[145]53. Growth arrest and DNA-damage inducible protein (gadd45aa, gadd45ba, gadd45bb), co-expressed in both brain and liver transcriptome along with fosab and jun, are associated with both MAPK signaling and apoptosis pathway, which seems to be induced simultaneously in L. catla post chronic CYP exposure (Supplementary Fig. 3) indicating the toxic effect of CYP even at a low concentration (1/10th of LC[50]). gadd45a was also found to be up-regulated on alpha CYP exposure to study the induction of apoptosis and oxidative stress in human neuroblastoma cell line, SH-SY5Y^[146]54. The DEGs btg2, sik1, fosab, gadd45ba and junbb co-expressed in both brain and liver transcriptome of L. catla were also among the identified co-expressed DEGs in the brain transcriptome of zebrafish after three and twenty-one days of post mild traumatic brain injury^[147]55. In the liver transcriptome of L. catla following chronic exposure to CYP, in addition to metabolic pathways, steroid biosynthesis and terpenoid backbone biosynthesis were up-regulated, which are interrelated as terpenoids, the most prevalent lipids and starting materials for the production of secondary metabolites, including sterols and steroids and the majority of the genes up-regulated in this pathways are also common to the metabolic pathway (ko01100). Most of the crucial genes essential for terpenoid biosynthesis were up-regulated following CYP exposure (mvk, idi1, fdps, mvda, hmgcs1) encoding important enzymes of the mevalonate pathway which is a crucial pathway for the production of the precursors isopentenyl-pyrophosphate (IPP) and dimethylallyl-pyrophosphate (DMAPP), thought to be the only pathway for terpenoid backbone biosynthesis. DMAPP and IPP were then converted into farnesyl-pyrophosphate (FPP), which was then further catalysed by farnesyl diphosphate synthase (FDPS), encoded by the fdps gene to generate terpenoids^[148]56. Similar to the present findings, most of these genes were up-regulated following salinity stress in the liver transcriptome of Cynoglossus semilaevis (half-smooth tongue sole)^[149]57. The steroid biosynthesis pathway was one of the significantly enriched pathways on chronic exposure to CYP in L. catla liver transcriptome, and some major genes associated with this biosynthesis pathway were up-regulated (dhcr24, tm7sf2, cyp51, sc5d, cyp27b1, ebp, msmo1, lss, hsd17b7, fdft1, dhcr7). fdft1 encoding the enzyme farnesyl-diphosphate farnesyltransferase catalyses the first committed step in sterol biosynthesis followed by lss (lanosterol synthase) responsible for the synthesis of lanosterol, cyp51 (sterol 14 alpha-demethylase) for the production of a 14-demethylation product of lanosterol, tm7sf2 (delta14-sterol reductase) associated with sterol reduction and hsd17b7 (3beta-hydroxysteroid 3-dehydrogenase) regulating the biosynthesis of zemosterol and cholesterol^[150]57. Interestingly, it was previously reported that in TM7SF2 mutant mouse, the normal functioning of cyp51 and msmo1 is compromised due to positional and expressional regulation by tm7sf2^[151]58. The Bloch pathway’s final phase, or the initial step of the Kandutsch-Russell (K-R) pathway’s cholesterol synthesis, is catalysed by the enzyme DHCR24 (encoded by dhcr24). In addition, theoretically, DHCR24 can operate on any intermediate, including desmosterol and lanosterol, to move intermediates from the Bloch pathway to the K-R pathway. In addition, DHCR24 can work in concert with DHCR7 (dhcr7), the last essential enzyme in the K-R pathway, to regulate its activity, ensuring coordinated control of cholesterol synthesis^[152]59. Studies carried out in vitro also showed that DHCR24 knockdown resulted in the suppression of cholesterol biosynthesis, the reduction of plasma membrane cholesterol, and the amount of intracellular cholesterol and also plays a central role in the pathogenesis of Alzheimer’s disease^[153]60. The critical genes associated with the steroid biosynthesis pathway identified in this study were also found in Oreochromis niloticus and Cynoglossus semilaevis on transcriptomics study post salinity stress^[154]57,[155]61. Pathway enrichment analysis of DEGs identified in the liver transcriptome revealed that the PPAR signaling pathway, along with fatty acid metabolism be significantly down-regulated on long-term CYP exposure to L. catla, and the associated down-regulated genes were mainly associated with fatty acid biosynthesis, i.e. fads2 (K10226) and scd ([156]K00507) (co-expressed in both the pathways) and elovl6 (K10203). These genes are mainly responsible for long-chain polyunsaturated fatty acids (LC-PUFAs) biosynthesis, such as arachidonic acid (C20:4, ARA), eicosapentaenoic acid (C20:5, EPA) and docosahexaenoic acid (C22:6, DHA) produced internally from C18 precursors through a chain of desaturation and elongation events that are catalysed by the enzymes fatty acyl elongase (elovl6) and desaturase (fads2 and scd), respectively^[157]62. A downregulation of these enzymes post-CYP exposure in L. catla indicates a significant decrease in unsaturated fatty acid levels, which may be due to the usage of fatty acids to meet the energy demands caused by the stress the pesticide generated^[158]18. Long-term exposure to CYP, even at low doses, can also elucidate toxic effects on the aquatic environment as well as its inhabitants as observed in L. catla brain and liver transcriptomic study, including (1) induction of apoptotic pathway in both brain and liver tissues (2) over-expression of steroid biosynthesis and terpenoid backbone pathways in response to CYP induced stress. Therefore, exposure to CYP, even at low concentrations, could alter the aquaculture environment, make it more difficult to produce aquatic goods safely, or even be detrimental to humans because it may enter the food chain. Conclusions The present study highlighted the brain and liver transcriptome’s response to long-term exposure of CYP to L. catla using next-generation sequencing. Pathway enrichment analysis of DEGs showed MAPK signaling, apoptosis, steroid and terpenoid backbone biosynthesis to be significantly enriched. It further suggested that CYP is toxic even at low concentrations as it induced cellular apoptotic pathways and stimulated steroid biosynthesis in response to the toxic shock. Therefore, the present investigation elucidated the mechanism underlying CYP-induced toxicity in L. catla and also identified functional biomarkers for transcriptional studies (fosab, nr4a1, tm7sf2, dhcr24) for early risk assessment, which have implications in studying other toxic pesticides prevalent in the aquatic ecosystem and also in other teleosts. Electronic supplementary material Below is the link to the electronic supplementary material. [159]Supplementary Material 1^ (2MB, docx) Acknowledgements