Transcriptome analysis of the brain provides insights into the regulatory mechanism for Coilia nasus migration

Background Coilia nasus (C. nasus) is an important anadromous fish species that resides in the Yangtze River of China, and has high ecological and economical value. However, wild resources have suffered from a serious reduction in population, attributed to the over-construction of water conservancy projects, overfishing, and environmental pollution. The Ministry of Agriculture and Rural Affairs of the People’s Republic of China has issued a notice banning the commercial fishing of wild C. nasus in the Yangtze River. Wild C. nasus populations urgently need to recover. A better understanding of C. nasus migration patterns is necessary to maximize the efficiency of conservation efforts. Juvenile C. nasus experience a simultaneous effect of increasing salinity and cold stress during seaward migration, and the brain plays a comprehensive regulatory role during this process. Therefore, to explore the early seaward migration regulation mechanism of juvenile C. nasus, we performed a comparative transcriptome analysis on the brain of juvenile C. nasus under salinity and cold stress simultaneously. Results Relevant neurotransmitters, receptors, and regulatory proteins from three categories of regulatory pathway play synergistic regulatory roles during the migration process: neuronal signaling, the sensory system, and environmental adaptation. The significant differential expression of growth-related hormones, thyroid receptors, haptoglobin, and prolactin receptors was similar to the results of relevant research on salmonids and steelhead trout. Conclusions This study revealed a regulatory network that the brain of juvenile C. nasus constructs during migration, thereby providing basic knowledge on further studies could build on. This study also revealed key regulatory genes similar to salmonids and steelhead trout, thus, this study will lay a theoretical foundation for further study on migration regulation mechanism of anadromous fish species.


Background
The Coilia fish belongs to the family of Engraulidae and the order of Clupeiforme, and is distributed in the midwest Pacific and Indian oceans. As a popular Coilia fish species for consumers in China, Coilia nasus (C. nasus) is a precious fish species in the Yangtze River. It is one of the "Three Delicious Species in the Yangze River", with Reeve's shad (Tenualosa reevesii) and obscure pufferfish (Takifugu fasciatus) being the other two species [1,2]. However, it has suffered from a serious population reduction in recent years as a result of the overconstruction of water conservancy projects, overfishing, and environmental pollution [3][4][5]. Consequently, the catch yield has reduced by 60% and continues to drop yearly [6]. It has been included on the "National Key Protective Species List" of China. The Ministry of Agriculture and Rural Affairs of the People's Republic of China has issued a notice banning the fishing of wild C. nasus in the Yangtze River for production. The restoration of wild C. nasus is urgently needed.
C. nasus is an important anadromous fish species. In February, mature adults return to their native Yangtze River and its tributaries to spawn. Their offspring move to the estuaries, where they will remain until autumn, and then migrate to the ocean for growth and fattening [7,8]. Therefore, during this process, juvenile C. nasus is simultaneously exposed to increased salinity and cold stress. There has been very few research on regulation mechanism of C. nasus during migration, which were mainly on regulatory pathways and function of key regulatory genes that function during spawning migration, such as the comparative transcriptome analysis on brain and liver of wild adult C. nasus during spawning migration [9] and function analysis on FoxL2 and Cyp19a1of C. nasus during anadromous migration [10]. The results indicated that many neurotransmitter signaling pathways in brain and relevant receptors, transporters, and regulatory proteins were significantly upregulated. Meanwhile, most pathways in liver were downregulated and indicated its function in energy conservation during spawning migration. The brain serves as the center of the nervous system in vertebrates and exerts a more comprehensive regulatory function than other tissues of perception system regulation, learning, and memory muscle activity, through which the organism responds to the changing environment [11,12]. Therefore, research on the influence of environmental factor variation on the brain transcriptome will be beneficial for revealing the comprehensive regulatory network that is formed during C. nasus migration.
Traditionally, research on the effects of temperature and salinity as environmental stressors in fish has been carried out in the liver and gills due to the pivotal roles of these organs in energy supply and osmoregulation.
Recent studies that investigated the strengthening of the brain regulatory function in response to salinity and cold stress have indicated that the expression of hormones, neurotransmitters, receptors, and key regulatory proteins was upregulated [13][14][15][16][17][18]. Xu et al. [19] investigated the effect of cold exposure on the brain transcriptome of the Yellow rum (Nibea albiflora). The results indicated that the most significantly enriched pathway was involved in signal transduction. Salmonids, such as Atlantic salmon (Salmo salar), coho salmon (Oncorhynchus kisutch), and steelhead trout (Oncorhynchus mykiss gairdneri), in addition to C. nasus, are also economically important anadromous fish species. In order to explore their regulatory mechanisms during smoltification, some research has been carried on trout, and resident and migratory salmonids, including comparative transcriptome analyses of the brain, liver, gill, kidney, and olfactory rosettes [20][21][22][23][24]. The results of these analyses indicated that differentially expressed genes (DEGs) were mainly involved in development and metabolism [20,21]. Relevant research on Atlantic salmon indicated that DEGs were involved in electron transport, oxygen transport and endocrinology, there was no change in the expression of thyroid-stimulating hormone (TSH), which is different from the results of similar research on steelhead trout and coho salmon [20,[22][23][24]. Additionally, a comparative transcriptome analysis on coho salmon in freshwater and early marine environments showed that differential regulatory pathways in the brain were mainly involved in protein synthesis and MHC1-mediated antigen presentation [24]. These studies indicated that anadromous fish species have differential regulatory mechanisms during seaward migration. Therefore, it is essential to explore the regulatory patterns in different anadromous fish species to reveal the potential universal regulatory mechanisms.
Research on the regulatory mechanism of C. nasus during migration is still in its infancy. Juvenile C. nasus seaward migration is an important part of the species' life history, but relevant research has not been carried out. Given the simultaneous effects of salinity and cold stress that juvenile C. nasus experiences during seaward migration, we performed a comparative transcriptome analysis of the brain under saline and cold stress, to investigate the regulatory role that the brain of juvenile C. nasus plays during migration. We aimed to reveal key regulatory pathways and genes, in order to construct a regulatory network; lay the theoretical foundations for further research on regulatory mechanisms during C. nasus migration and for the optimization of artificial breeding of C. nasus, which is beneficial for providing high-quality fry fish for proliferation and release; and contribute to efforts towards the restoration of wild C. nasus. This study will also lay a theoretical foundation for research on the regulation patterns of global Coilia fish during migration. Combined with existing reports on anadromous fish, this study will collect basic information on the regulation mechanism of anadromous fish species during migration.

Results
To comprehensively explore regulation mechanism of juvenile C. nasus during seaward migration, we performed comparative transcriptome analysis on juvenile C. nasus under saline and cold stress simultaneously. Top 10 GO terms, top 10 KEGG pathways and key DEGs were obtained after library construction, sequencing, data filtering, assembly, annotation and differential expression analysis. Correlation analysis on intraclass difference in the control and stressed group was made, validation of RNA-Seq data was carried out with quantitative real-time polymerase chain reaction (qPCR).

Transcriptome assembly and statistics of unigenes
The average RIN (RNA Integrity Number) for six brain samples was 9.5. After quality filtering, the RNA-Seq of six brain samples yielded around 46.36 million highquality sequence data. The Q value (Q30) was used as the cutoff for quality control. The Q30 values of the samples reached up to 93.03%, and the GC-content of each sample reached around 48.5% (Table 1). The clean reads obtained from the six transcriptome libraries were assembled to full-length transcripts, and a total of 436, 325 unigenes were obtained after the elimination of redundant transcripts. The average transcript length was 795 bp, and N50 was 1001 bp. The average clean ratio reached 99.8%.
Correlation analysis on intraclass differences in the control and stressed group CORREL function was used to analyze difference of FPKM (Fragments per kilobase of transcript per million mapped reads) of DEGs in the three replicated groups of control group (C1-C3), as well as in the stressed group (E1-E3) (Additional file 5: Table S4). The correlation analysis results of C1-C2, C2-C3, C1-C3, E1-E2, E2-E3 and E1-E3 were as follows, y = 0.835x + 0.9861and R 2 = 0.8554 (correlation coefficient r = 0.924863193), y = 1.1849x -1.2712 and R 2 = 0.9373 (r = 0.968150821), y = 1.0599x -0.7987 and R 2 = 0.92 (r = 0.959142331), y = 0.8144x + 1.0789 and R 2 = 0.7603 (r=, 0.969855973), y = 0.9511x + 0.9828 and R 2 = 0.9081 (r = 0.935047937), y = 0.9119x + 0.204 and R 2 = 0.889 (r = 0.937179862). The results indicated that replicated groups in the control group had strong correlation, as well as in the stressed group, intraclass difference were both small in these two groups. These differences were mainly caused by the individual differences of experimental animals and operation difference during experiment, which are normal and acceptable difference. Therefore, Intra-group differences did not affect the further analysis on differences between the control and stressed group.

Top 10 gene ontology (GO) enrichment analysis on DEGs
Based on the GO enrichment analysis, 38,579 unigenes were categorized into 62 functional groups from three categories: biological processes (BP), molecular functions (MF), and cellular components (CC) (Additional file 1: Figure S1). Then, we conducted a further GO enrichment analysis on DEGs and obtained the top 10 GO terms from each of the three categories ( Fig. 1). Most BP terms, with the exception of some involved in the general function (Nos. 1, 2, 5, and 6), were related to neuronal signal transduction (Nos. 3 and 7) or the sensory system (Nos. 4 and 10). Most MF and CC terms were relevant to synaptic transmission or the sensory system and relevant components, such as neuropeptide binding, glutamate receptor activity, the synaptic vesicle membrane, the cell junction, retinol binding, photosystem I, the interphotoreceptor matrix, etc. The DEGs of each term are shown in Additional file 2: Table S1.

Top 10 Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis
In total, we obtained 4721 DEGs (Additional file 2: Table S1); 2020 DEGs were downregulated and 2701 DEGs were upregulated. As shown in Fig. 2, five pathways were involved in neuronal signaling, including neuroactive ligand-receptor interaction, the calcium signaling pathway, glutamatergic synapse, Clean ratio equals clean reads/raw reads retrograde endocannabinoid signaling, and the serotonergic synapse. Two pathways were related to the sensory system-olfactory transduction and phototransduction-and two were relevant to environmental adaptation-circadian entrainment and ECM-receptor interaction. The DEGs involved in these pathways are shown in Additional file 2: Table S1.

Functional analysis on DEGs
According to a pathway hierarchy (Additional file 3: Table S2), the top 10 GO and top 10 KEGG terms indicated that the DEGs were mainly involved in three categories: neuronal signaling, the sensory system, and environmental adaptation. The relevant DEGs are listed in Table 2.

Validation of RNA-Seq data by qPCR
Ten DEGs were randomly selected from the RNA-Seq data of upregulated and downregulated genes. As shown in Fig. 3, expression of the genes were normalized to beta-actin. The genes and primers used for qRT-PCR were shown in Additional file 4:  (Fig. 4), These ten DEGs exhibited a concordant direction in both the RNA-Seq and qPCR analyses. The results indicated that key pathways and DEGs were mainly involved in the neuronal signal transduction,

Discussion
In this study, we performed a comparative brain transcriptome analysis of C. nasus exposed to salinity and cold stress simultaneously. A comprehensive analysis of the top 10 GO and top 10 KEGG terms indicated that the most significantly enriched pathways were involved in three categories: neuronal signaling (neuroactive ligand-receptor interaction, the calcium signaling pathway, the glutamatergic synapse, retrograde endocannabinoid signaling, and the serotonergic synapse), the sensory system (olfactory transduction and phototransduction), and environmental adaptation (ECM-receptor interaction and circadian entrainment). Neurotransmitters, receptors, and key regulatory proteins were shown to be upregulated, constructing a regulatory network in the C. nasus brain that plays a role in C. nasus migration.

Quality of the de-novo assembled transcriptome
As shown in Table 1, the ratio of clean reads reached 99.8%, and the reported clean read ratio of high-quality de-novo assembly was around 96% [25,26]. To some extent, this indicated that we obtained a high-quality assembled dataset in this study. Due to the fact that the reads to be assembled in this study were from different individuals, this added complications, considering that Trinity has a self-correction ability [27]. We also took some actions to refine our dataset, such as filtering the longest open reading frames for each transcript, clustering similar transcripts, and making annotations on filtered transcripts via comparing them with closely related species. Therefore, the assembly result in this study is credible.

Neuronal signal transduction in C. nasus brains under salinity and cold effect
Glutamate is a major excitatory neurotransmitter of the central nervous system. It mediates rapid interneuronal communication and comparatively slow synaptic plasticity, whilst also regulating the ionotropic glutamate receptor (iGluR) and the metabotropic glutamate receptor (mGluR) [28]. In this study, salinity triggered the upregulated expression of glutamate dehydrogenase in the brain, which may be the reason for the upregulation of the glutamate receptor [13]. On the contrary, serotonin exerted a regulatory role in inhibitory neuronal signal transduction and maintained the excitability balance, together with glutamate [29]. The retrograde endocannabinoid signaling pathway was also significantly upregulated. Retrograde endocannabinoid signaling, which plays a homeostatic role in the balance of neuronal activity between the glutamate and inhibitory serotonergic neurotransmitters, is mediated by retrograde messengers, including anandamide and 2-arachidonoyl glycerol (2-AG). These signaling pathways have pivotal regulatory functions in synaptic plasticity, learning, and memory in response to environmental stress [30][31][32].
In addition to the above-mentioned pathways, many other neurotransmitters, receptors and regulatory proteins were upregulated, such as, galanin, neuropeptide K (NPK), neurotensin (NT), nociceptin and its corresponding receptors, cerebellin-1 (CBLN1), leucine-rich repeat transmembrane neuronal protein 4 (LRRTM4), regulating synaptic membrane exocytosis protein 1 (RIMS), signalinduced proliferation-associated 1-like protein 1 (SIPA1L1), Synaptosomal-associated protein 25 (SNAP25), Teneurin-1 (TENM1), etc. (Table 2). Of these compounds, there is evidence that galanin is involved in the regulation of coping with stressful events [33]. Neurotensin, a key neurotransmitter, plays a regulatory role in easing pain in the brain, for example, by regulating the aversive memory of temperature nociception. Nociceptin and its receptor also play roles in numerous brain activities, such as pain sensation and fear learning [34]. CBLN1, essential for synapse integrity and synaptic plasticity, plays a regulatory role in maintenance of excitatory synapses [35]. LRRTM4, it plays a pivotal role in the development and maintenance of nervous system of vertebrate [36]. RIMS, it plays a regulatory role in neurotransmitter release at the active zone during short-term synaptic plasticity and plays an important role in maintaining the normal probability of neurotransmitter release [37]. SIPA1L1, it functions in promoting reorganization of the actin cytoskeleton and plays a regulatory role in dendritic spine morphogenesis [38]. SNAP25, it plays a key regulatory role in synaptic function of specific neuronal systems, and associates with proteins involved in vesicle membrane fusion [39]. TENM1, as a signal transducer, it plays a regulatory role in the establishment of proper connectivity in the nervous system [40].

Regulation of the sensory system of C. nasus under salinity and cold effect
Salinity is a significant environmental factor during fish migration that can affect the expression of the olfactory receptors, which could be involved in olfactory imprinting. Particularly, for anadromous fish, the sensory system plays a key role in migration [17]. Several studies on the migration mechanism of Salmonids (sockeye salmon, Atlantic salmon, etc.) have concluded that the olfactory epithelium is essential for imprinting and subsequent homing migration as fish return to their natal streams and undergo rapid salinity changes [41,42]. Under salinity and cold stress condition, as a result, olfactory transduction pathway of juvenile C. nasus was significantly upregulated and many relevant regulatory genes were significantly enriched including cyclic nucleotide-gated olfactory channel, claudin, anoctamin, guanylyl cyclase-activating protein, guanylyl cyclase GC-E, etc. (Table 2). Further, the olfactory epithelium associated genes of juvenile C. nasus were also shown to be significantly elevated during salinity and cold stress. In contrast, in the Xu et al. [19] study on the Yellow drum (Nibea albiflora), the pathways involved in the sensory system were not upregulated, possibly because Yellow drum (Nibea albiflora) is a non-anadromous fish species. In C. nasus, as an anadromous fish, the sensory system plays an important role in guiding its returning behavior during migration, so its regulatory mechanisms are different from other non-anadromous fish.
Meanwhile, many visual system-relevant genes of juvenile C. nasus were shown to be significantly down-regulated during stress exposure, such as rod cGMP-specific 3',5'-cyclic phosphodiesterase subunit beta, retinol dehydrogenase 7, ventral anterior homeobox, and rhodopsin. Rod cGMPspecific 3',5'-cyclic phosphodiesterase, a protein involved in the transmission and amplification of visual signals, plays a vital role in the formation of a functional phosphodiesterase holoenzyme. The ventral anterior homeobox is a transcription factor involved in retinal ganglion cell axon guidance. Rhodopsin, a light sensitive G-protein-coupled receptor protein, enables vision in dark conditions [43]. Hale et al. [20] also carried out comparative transcriptome on brain in juvenile resident and migratory smolt rainbow trout (Oncorhynchus mykiss) also found that phototransduction signaling pathway and key regulatory genes were expressed significantly. The light-brain-pituitary (LBP) axis transmits information about environmental photoperiod to endocrine systems, the subsequent hormone actions regulate physiological behaviour such as parr-smolt tansformation, etc. [44].

Response regulation under salinity and cold effect
Somatostatin is a peptide hormone that inhibits growth via inhibiting the release of secondary hormones [45]. Cortistatin is a neuropeptide with a high structural similarity and similar function to somatostatin. Furthermore, cortistatin can bind with all somatostatin receptors to mediate its inhibitory effect on downstream growth-associated signaling pathways [46]. In this study, the somatostatin receptor was upregulated in the experimental group to perform an inhibitory role in growth, which is consistent with reports showing that high salinity stress inhibits the growth of tilapia fish [47]. This may be because fish obtain a better energy distribution in response to stress via growth suppression. Heat shock protein 70 (HSP70), a main molecular chaperone, combines with unfolded or misfolded proteins and also promotes the degradation and removal of denatured proteins [48].
In this study, HSP70 was significantly upregulated, and many studies have indicated an elevation of HSP70 in response to environmental stress [49,50]. HSP70 is also regarded as having the function of inhibiting growth, as well as playing a synergistic role with somatostatin and cortistatin. Research on the regulation mechanism of the brain during salmonid and steelhead trout smoltification showed the differential expression of growth hormone and thyroid hormone receptors [20,22]. This study also showed the significant differential expression of hormones and receptors related to growth regulation, including cortistatin and the somatostatin receptor, which indicated the importance of growth regulation during anadromous fish migration. Haptoglobin, a stress-activated acid glycoprotein, can bind to free plasma hemoglobin and present itself to the mononuclear macrophage system for processing, in order to maintain homeostasis of the blood micro-environment under stressful conditions [51]. In this study, the haptoglobin concentration increased significantly after stress. A comparative transcriptome analysis of the brain of Atlantic salmon during smoltification showed a significant upregulation of hemoglobin [21]. In this study, hemoglobin was not upregulated, but haptoglobin, which has a regulatory function, was upregulated. Given its binding and presenting function on hemoglobin, this may indicate that Atlantic salmon and C. nasus have a similar regulatory pattern in maintaining homeostasis of the blood micro-environment during seaward migration.
In this study, many osmoregulatory-related genes were significantly upregulated, such as claudins, the sodium-and chloride-dependent glycine transporter 2, the solute carrier family (SLC), and the aquaporin and prolactin receptors. Sodium-and chloride-dependent glycine transporter 2 is necessary for Na + /Cl --dependent neurotransmitter transportation. SLC is a group of membrane transport proteins, predominantly located in the cell membrane, which play roles in osmoregulation. Aquaporins are integral membrane proteins that control the movement of water in and out of cells and function in osmoregulation [52]. Wang et al. [53] reported that Aquaporin I performs an important regulatory function in Coilia nasus under conditions of salinity stress. Prolactin has also been reported to have a regulatory function when freshwater fish swim into seawater [54]. In this study, the prolactin receptors were upregulated to assist with the anti-stress response of C. nasus. Prolactin is considered to be an important hormone for salmonid and steelhead trout during the smoltification adaptation process. In this study, the significant up-regulation of the prolactin receptor in C. nasus also reflected the important regulatory role of this gene during anadromous fish seaward migration.
Cold stress causes damage to fish through the accumulation of free radicals. Methionine is an effective scavenger for free radicals. It is also one of the methyl donors and participates in the synthesis of adrenaline, creatine, choline, carnitine, and DNA methylation. Inadequate methionine will affect the metabolism of fish and thus weaken its cold resistance [55]. In this study, many genes in the methionine and cysteine receptor pathways were significantly elevated, which also indicated the important role that the methionine metabolism of juvenile C. nasus plays in response to cold stress.
According to the theory of the phase transition of membrane lipids, it is believed that the content of unsaturated fatty acids in membrane lipids can affect membrane fluidity and, thus, the cold resistance of fish [56]. Beyond that, cholesterol also plays an important role in maintaining the fluidity of biological membranes [57]. Sterol, a cholesterol-like substance, is widely regarded as a key functional substance in the response to environmental stresses such as extreme temperature and high salinity stress [58,59]. Sterol 14-demethylase is a key synthetase in the synthesis of cholesterol and phytosterol [60]. In the present study, its expression level was significantly elevated after stress, indicating that the cholesterol level was elevated to enhance the anti-cold resistance of juvenile C. nasus.

Conclusions
In this study, we discovered that neuronal signal transduction played an important regulatory role in maintaining excitability homeostasis. Additionally, the function of the olfactory system was strengthened and the visual signaling pathways were weakened during juvenile C. nasus seaward migration. The identified significant differential expression of the growth hormone and receptors; the upregulation of the genes which function in maintaining homeostasis of the blood environment, such as haptoglobin and the prolactin receptors, are similar results to those obtained in research on salmonids and steelhead trout smoltification adaptation. This study revealed an early regulation pattern during juvenile C. nasus seaward migration as well as some potential universal regulatory genes which are important during anadromous fish migration. However, this study is more beneficial in terms of revealing the regulation mechanism during juvenile C. nasus early seaward migration, and provides a glimpse of the regulation mechanism during C. nasus seaward migration. At present, there has not been much research on the regulation mechanism of juvenile C. nasus seaward migration, so the pathways and genes identified in this study represent a starting point for further research, and do not represent an exhaustive list of all important candidate genes. Therefore, the regulation mechanism of juvenile C. nasus during seaward migration requires more in-depth studies.

Experimental fish and synergistic stress of salinity and low temperature
Experimental fish were five-month-old juvenile C. nasus (sex haven't been distinguished) whose average body length was 133.5 ± 0.65 mm and body weight was 9.66 ± 0.93 g. The fish were collected from Yixing, an experimental base of the Freshwater Fisheries Research Center of the Chinese Academy of Fishery Sciences, in November. Wild juvenile C. nasus always migrate to the ocean during this time. One month before the experiment, juvenile C. nasus were netted and randomly transferred to six aquariums with circulating water systems whose volume was 1096 × 470 × 670 mm 3 . Each aquarium contained six fish. The aquariums were continuously aerated and the water quality was monitored daily: the water temperature was 20 ± 0.95°C, the pH was 7.2, and the concentration of dissolved oxygen was 7.9 ± 0.63 mg/ L. The experimental fish were fed with compound feed three times each day, at 7:00 a.m.,12:00 p.m., and 6:00 p.m. Feeding stopped 1 day before the stress experiment.
The control group (C) and the stressed group (S) were established. The stressed group conditions included a salinity level of 15‰, achieved by adding sodium chloride, and a temperature of 15°C, achieved by adding crushed ice. Salinity was monitored by using a SALT6 salimeter (Eutech/Oakton, Illinois, USA). The water temperature was monitored by using a RHXL3SD thermometer (OMEGA, Elgin, USA). As shown in Fig. 5, each group had three replicates. It was observed that juvenile C. nasus swam gently after 3 h stress treatment., showing an obvious difference in behavior from the initial stress period. Therefore, we sampled at 3 h to discuss internal regulation and adaptation mechanism of juvenile C. nasus after stress. Two fish were sampled from each aquarium at 3 h after the stress conditions were induced. A barrel containing 40 mg/L MS-222 (Kuer Bioengineering, Beijing, China) was prepared in advance. Fish were netted quickly and placed in it for 10 s for immersion anaesthesia at the time of collection at 10 am. We cut an incision from fish's neck for bloodletting to sacrifice fish, then we cut the skull above brain from the incision with scissors, took out the brain with forceps, put it in liquid nitrogen for snap-frozen, and then stored at -80°C.
Total RNA extraction, cDNA library construction, and Illumina sequencing Total RNA was extracted from each brain sample using RNAiso reagent (Takara, Kusatsu-Shiga, Japan), in accordance with the manufacturer's instructions. Then, equal amounts of total RNA of individuals from each aquarium were pooled. Finally, six RNA brain samples were obtained (Fig. 5), and the amount of RNA of each sample was 1 μg. The concentration and quality of extracted RNA were checked on an Agilent Bioanalyzer 2100 with RNA 6000 Nano Labchips (Agilent technologies, Santa Clara, USA), and contaminant genomic DNA was removed with Recombinant DNaseI (Takara, Kusatsu-Shiga, Japan). The mRNA was isolated with magnetic beads and then broken into fragments and reverse transcribed into cDNA with added adapters. Finally, six cDNA libraries were constructed and sequenced on the Illumina HiSeq™ 2500 platform (Illumina, San Diego, USA). Paired-end data were used for further analysis. The obtained raw data were submitted to NCBI (NCBI, Bethesda, USA) with the accession number SRP078492.

Data filtering and assembly
The raw data were tested for quality control using FASTQC (Babraham Institute, Cambridge, UK). Some low-quality vectors (including adapters and/or primers), contaminated reads, low-quality bases at the 3' end, empty reads, and ambiguous 'N' nucleotides were removed, and the cutoff value for length control was set as 35 bp. NGS QC TOOLKIT v2.3.3 (Roche, Pleasanton, USA) (http://59.163.192.90:8080/ngsqctoolkit/) [61] was used to filter the above-mentioned data. Transcriptome assembly was performed using the Trinity software (Broad institute, Cambridge, UK) [27].
Function annotation, gene quantification, and differential expression analysis The similarity alignment was based on the BLAST algorithm (NCBI, Bethesda, USA). For homology annotation, the obtained non-redundant sequences were aligned in the following priority order: non-redundant protein (Nr), non-redundant nucleotides (Nt), Swiss-prot (http://www. uniprot.org/downloads), clusters of orthologous groups for eukaryotic complete genomes (KOG, ftp://ftp.ncbi. nih.gov/pub/COG/KOG/kyva), gene ontology (GO, Fig. 5 Diagram of the sampling method. C1-C3 indicated three replicated groups in the control group, S1-S3 indicated three replicated groups in the stressed group http://www.geneontology.org/), and the Kyoto Encyclopedia of Genes and Genomes.
Differential expression quantification was calculated using the DESeq software package (http://bioconductor. org/packages/release/bioc/html/DESeq.html) [68]. Fold change was calculated as the ratio of the expression level of genes in the stressed group and control group samples, and |log 2 foldchange| ≥ 1 and FDR < 0.05 (Benjaminie-Hochberg false discovery rate) [69] were set as the cutoff thresholds to determine the significantly differentially expressed genes. GO and KEGG pathway enrichment analyses were performed on differentially expressed genes (FDR < 0.05). The GO and KEGG terms were sequences with -log10 (p-value), and terms with a list hit of less than two were excluded (list hit refers to number of DEGs enriched in one GO or KEGG term). Finally, we obtained the top 10 GO terms and top 10 KEGG pathways in the brain. According to the pathway hierarchy that the KEGG website published (Additional file 3: Table S2), the top 10 KEGG pathways were further categorized. Considering the top 10 GO terms, the top 10 KEGG pathways, key pathways, and DEGs relevant to the anti-stress response, a regulatory network was constructed.

Quantitative real-time PCR validation
A quantitative real-time polymerase chain reaction (qPCR) was performed with the extracted RNA to validate the accuracy of the transcriptome sequencing analysis by the ABI 7500 Real-Time PCR System (ABI, New York, USA) using the Power SYBR™ Green PCR Master Mix. Ten differentially expressed genes were randomly selected with the Randbetween Function and the primers were designed with Primer Premier 5 software (Biosoft, California, USA). β-actin was used as the internal reference. The amplifications were conducted with the following procedure: 95°C for 30 s, 40 cycles of 95°C for 5 s, 60°C for 34 s, and 72°C for 50 s. Every sample was analyzed in triplicate, and the permitted variation of the cycle threshold in replicates was 8%. The 2 -ΔΔCT method was used to calculate the gene expression level [70].

Statistical analysis
The statistical analysis for qPCR experimental data was completed with SPSS 21.1 software (SPSS, Chicago, USA), and values are shown as means ± SD. All data were subjected to one-way ANOVA and Tukey's multiple range tests. The α-level was set at 0.05.