Transcriptome and Molecular Pathway Analysis of the Hepatopancreas in the Pacific White Shrimp Litopenaeus vannamei under Chronic Low-Salinity Stress

We obtained 97.1 million reads and 98.1 billion nucleotides from the shrimp at both 3 and 30 practical salinity units (psu) (Table 1). After quality trimming and adapter clipping, 93.8 million reads accounting for 96.6% of the total reads were obtained. In addition, 26,034 genes and 38,237 unigenes with an average length of 1,610 bp were obtained by de novo assembly using Trinity software after splicing and removing redundancy. Among 26034 genes, 855 genes were significantly up or down regulated, and all the differential expression genes were show in S1 Table. A total of 2,341 unigenes were down-regulated, and 2,363 unigenes were up-regulated in low salinity with the absolute fold change >2. Among the unigenes, the largest and smallest unigenes were 24,554 bp and 351 bp, respectively (Table 2); the length distribution of unigenes is shown in S1 Fig. To assess the quality and coverage of the transcriptome data, we mapped the assembly unigenes by using Bowtie software, and the mapping data accounted for over 90% of the data, suggesting that the transcriptome dataset provided good gene coverage and enriched the transcriptome information of L. vannamei in the present study.

The predicted sequences (predicted open reading frame nucleotide sequences) and unpredictable sequences (unpredictable nucleotide sequences) were annotated using BLASTp and BLASTx, respectively, and then blasted to protein databases, including the nonredundancy (NR), STRING, COG and KOG databases (BLAST 2.2.28+, E-value < 1e—5) (Table 3). Among the annotated predicted sequences, a total of 17,232 (76.83%), 6298 (28.08%), 3720 (16.58%) and 302 (1.35%) sequences were unambiguous alignments relative to the reference when BLASTx was used against the NR, STRING, KOG, and COG databases, respectively. However, among the unpredictable sequences, only 2,235 (14.14%), 746 (4.72%), 509 (3.22%), 313 (1.98%), 243 (1.54%), 128 (0.81%), and 56 (0.35%) of the 15,806 sequences were matched against the NR, gene ontology (GO), NT, STRING, KOG and COG databases, respectively.

After parsing the GO annotation output, a total of 8,779 unigenes were finally annotated with GO terms, with 50.02% annotated to biological processes, 28.93% annotated to cellular components and 21.05% annotated to molecular functions. The distribution of GO terms showed that cellular processes and metabolic processes were the well-represented terms among the biological processes. Cells and cell parts were significantly enriched in cellular components, and analytic activity and binding consisted of a large proportion of molecular functions (). S2 Fig

Clusters of orthologous groups (COGs) of proteins were determined to predict and classify the possible functions of unigenes. The COG annotation analysis showed that three types of function were obtained, including information storage and processing, cellular processes and signaling, and metabolic pathways. The hits from the COG analysis were functionally classified into 25 categories, and the most enriched terms were related to general functions and then transcription and signal transduction mechanisms (). S3 Fig

Various molecular pathways were obtained by Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation. A total of 9,621 unigenes were mapped onto 317 pathways, and the most enriched sequences were metabolic pathways, which was followed by the biosynthesis of secondary metabolites, spliceosome and RNA transport. The top 20 pathways with the greatest number of annotated sequences are shown in Table 4. A total of 47 significantly changed GO terms were obtained, and the most significant change was in molecular functions, which was followed by catalytic activity, histone H4 acetylation, structural constituents of cuticle and chitin binding. All of the significantly changed (P < 0.05) GO pathways are listed in Table 5. A KEGG pathway enrichment analysis was performed for the gene expression between salinity treatments to identify the number of significantly changed samples along the pathway that were relevant to the background number. The most significantly changed KEGG pathways were the glycosphingolipid biosynthesis, lysine degradation, glycosaminoglycan biosynthesis and malaria pathways. The gene recorded as B3GNT1,2 was significantly up-regulated both in the glycosphingolipid biosynthesis pathway and glycosaminoglycan biosynthesis pathway when shrimp were exposed to salinity at 3 psu (Figs 1 and 2). In addition, fatty-acid biosynthesis (Fig 3) was significantly enhanced, especially in short-carbon-chain fatty acids (C8-C18). In addition, low-salinity conditions enhanced polyunsaturated fatty-acid (PUFA) biosynthesis, especially that of highly unsaturated fatty acids such as ARA, EPA and DHA (Figs 4 and 5). All of the significantly changed (P < 0.05) KEGG pathways are listed in Table 6. These annotations provide valuable information for studying the specific biological and metabolic processes and functions and molecular mechanisms under salinity stress in L. vannamei.

Twenty randomly selected genes were measured in the same hepatopancreas RNA samples by qPCR, and the expression levels of these genes were significantly associated with the RNA-seq results (R = 0.91, Fig 6). These results confirm the reliability of RNA-seq and accuracy of the Trinity assembly.

Shrimp gill is one of the main tissues for osmoregulation and icon exchange, and the polyunsaturated fatty acids synthesized in hepatopancreas are the main components of shrimp gill. Previous study has shown that under salinity stress, shrimp gill needs extra energy (most from saturated fatty acids) and polyunsaturated fatty acids to ensure the osmoregulation and icon exchange in L. vannamei [38]. In this study, a number of lipid metabolism pathways in the hepatopancreas of L. vannamei were influenced by salinity stress, including the fatty-acid biosynthesis pathway, arachidonic acid metabolism pathway (S4 Fig), adipocytokine signaling pathway (S5 Fig), glycerophospholipid metabolism pathway (S6 Fig), ether lipid metabolism pathway (S7 Fig) and ketone body synthesis and degradation pathway (S8 Fig). These pathways are found in the Chinese mitten crab Eriocheir sinensis under osmotic stress [31]. Lipids, including the fatty-acid structure, composition and metabolism, reduce osmotic shock in aquatic animals (Sui et al., 2007) by providing sufficient energy to maintain the ion balance and regulate the structure of biological membranes [6,39]. This process can directly or indirectly function in “compensatory processes” and “limiting processes," which are the two major osmoregulation strategies in crustaceans [4, 40].

In this study, most of the significant pathways related to salinity adaptations were associated with the above two major strategies. The limiting process is a strategy to adjust the permeability of the boundary structures to maintain hemolymph osmolality/ions in gill membranes, which is effective at reducing ion diffusion and water influx because the ion transport mechanism requires additional energy [4]. The ability to regulate gill permeability is an adaptive response that is crucial for decapod crustaceans because it facilitates long-term habitation in environments with variable salinity [4, 40]. Therefore, the cell membranes of the gill should play an important role in osmoregulation [41]. The “compensatory process” strategy is accomplished via the active exchange of solutes in hemolymph to counterbalance passive diffusion and maintain osmolality [4]. Because this process is energetically costly, a number of energy metabolism pathways and ion regulation pathways must be involved.

Fatty-acid biosynthesis was significantly enhanced, especially for short-carbon-chain fatty acids (C8-C18), and arachidonic acid was also used to convert fatty acids or other metabolic products in the present study. When shrimp suffer from ambient salinity stress, additional energy is required through nutrient intake to maintain homeostasis by osmoregulation via the “compensatory process,” in which lipids play significant roles [42–45].

However, ketone bodies are indispensable for energy that is metabolized from fatty acids, which was revealed by the KEGG analysis in this study. Ketone bodies are produced by the liver from fatty acids during periods of low food intake or carbohydrate restriction. When carbohydrates are scarce, energy must be obtained from the breakdown of fatty acids from body tissue instead of glucose [46, 47]. Interestingly, in the adipocytokine signaling pathway, we found that the gene FACS functioned in long-chain fatty-acid biosynthesis and was significantly down-regulated under low-salinity stress. CPT-1 is another relevant gene for long-chain fatty-acid β-oxidation, and it was significantly down-regulated. Therefore, it is reasonable to assume that L. vannamei prefers to use shorter-chain fatty acids for energy supplementation and selectively retains longer-chain unsaturated fatty acids [39, 45]. The enhancement of saturated fatty acid biosynthesis in hepatopancreas (the main site for lipogenesis) would provide sufficient energy for osmoregulation in shrimp gill at low salinity, and this result is consistent with our previous findings in L. vannamei cultured at low salinity [38]. Penaeus monodon also prefers to use shorter-chain fatty acids in energy metabolism and selectively retains longer-chain unsaturated fatty acids [48].

In this study, fatty acid elongation (n > 18) was significantly enhanced, and the biosynthesis of α-linolenic acid, EPA, DHA, arachidic acid, arachidonic acid and other PUFAs was up-regulated at low salinity. Saturated fatty-acid biosynthesis was enhanced, and PUFA biosynthesis was also strengthened at low salinity. PUFAs can improve the resistance to osmotic shock in aquatic animals because PUFAs are mainly incorporated in cell membranes and can increase membrane permeability and fluidity [49, 50]. Free fatty acids, especially LC-PUFA, have the potential to modulate fatty-acid composition in the gill membrane and increase enzymatic efficiency [45, 51]. Modifications to the fatty-acid composition in the gills that increase the level of LC-PUFA can increase the gill area to enhance the osmoregulatory capacity in shrimp at low salinity, thereby increasing survival [39]. In addition, arachidonic acids in fish can enhance the branchial Na⁺/K⁺-ATPase activity and influence the ion balance [52]. Furthermore, arachidonic acid metabolites can influence the regulation of ion balances across the gill membrane [53].

In this study, the most influential pathways were associated with glycosphingolipid biosynthesis, glycosaminoglycan biosynthesis and glycosaminoglycan degradation, and they are all related to gill permeability. Glycosphingolipids (including lacto-/neolactoglycolipids and sphingolipids) function to protect the cell surface by maintaining the stability of the membrane or plasma membrane to protect against harmful environmental factors by forming a mechanically stable and chemically resistant outer leaflet for the plasma membrane lipid bilayer [54, 55]. Glycosaminoglycans are classified into four groups: keratan sulfate, heparin/heparan sulfate, chondroitin/dermatan sulfate and hyaluronic acid. One of the main functions of keratan sulfates is the maintenance of tissue hydration and implantation and migration of endothelial cell [56].

When shrimp were exposed to salinity of 3 psu, the B3GNT1,2 gene, which participates in both the glycosphingolipid biosynthesis pathway and glycosaminoglycan biosynthesis pathway, was significantly up-regulated. The FUT1_2 gene was also significantly up-regulated in the glycosphingolipid biosynthesis pathway. Thus, we believe that the genes B3GNT1,2 and FUT1_2 can promote glycolipid biosynthesis (Fig 1). Glycolipids constitute the lipid bilayer of the plasma membrane and play an important role in maintaining the physical state of the membrane [57]. Thus, glycolipid biosynthesis might promote the “limiting process” and provide resistance to low-salinity stress at 3 psu salinity. However, the effect of keratan sulfate and glycosaminoglycan degradation in low-salinity ambient osmoregulation is still not clear and requires further study.

Glycerophospholipids are glycerol-based phospholipids that are the main component of biological membranes [58], and their biosynthesis was significantly enhanced at low salinity (Fig 1). Although direct interactions were not observed between glycerophospholipid and osmoregulation in this study, glycerophospholipids could indirectly improve osmoregulation by affecting membrane permeability. Ether lipids are ubiquitous and constitute a major portion of the cell membranes in mammals [59]. In addition, these lipids play an important role in the generation of lipid second messenger systems, such as prostaglandins and arachidonic acid, which are important in signal transduction [60]. Ether lipids also act directly in cell signaling [61] and are involved in osmotic stress signal transduction. Another possible function of ether lipids is as an antioxidant against oxidative stress, which has been demonstrated in cell culture under salinity stress. Therefore, these lipids might play a role in serum lipoprotein metabolism [62] and lipid metabolism may play an important and indispensable role in osmoregulation [38].

Because of the complexity of the physiological response to low-salinity stress in L. vannamei, several pathways (in addition to lipid metabolism) show potential importance in the shrimp's ability to cope with salinity stress. However, clear evidence of the direct involvement of these pathways during salinity challenges has not been observed; thus, the putative functions of these pathways including lysine degradation, cholinergic synapse, drug metabolism pathway, steroid hormones metabolism pathway, phosphonate and phosphinate metabolism, are only briefly discussed.

Lysine is metabolized in eukaryotes to yield acetyl-CoA via lysine acetylation [63, 64]. Acetyl-CoA participates in osmoregulation as an intermediate metabolite commonly produced from energy metabolism, lipid metabolism and carbohydrate metabolism. In our results, the gene expression of 3-hydroxyacyl-CoA dehydrogenase, which can produce acetyl-CoA, was significantly up-regulated, resulting in more acetyl-CoA entering the citrate cycle for energy production. The above pathway could indirectly influence ion transfer or energy metabolism and promote “compensatory processes," which is consistent with reports on the Chinese mitten crab under salinity stress [31]. On the other hand, the choline in animal tissues is a primary component of neurotransmitter acetylcholine and functions with inositol as a basic constituent of lecithin [65]. Choline prevents the formation of fat deposits in the liver and facilitates the movement of fat into cells [66]. Thus, considering the result of our previous study [38], it seems that more fat may be used in supplying energy through β-oxidation.

Various physiological and pathological factors can affect drug metabolism, including age, individual variations, enterohepatic circulation, nutrition, intestinal flora and sex differences. Cytochrome P450 influences arachidonic acid metabolism [67, 68] and fatty-acid metabolism [69] and may have an indirect influence on osmoregulation by influencing the regulation of fatty-acid metabolism and other physiological and biochemical processes. However, the effect of drug metabolism on osmoregulation is still unknown and requires further study. Among the malaria pathways, we found that low-density lipoprotein receptor-related protein 1 (LRP1) was significantly changed when shrimp suffered from salinity stress. LRP1 might have an important role in low-density lipoproteins and thereby influence osmoregulation.

Steroid hormones help control metabolism, inflammation, immune functions and salt and water balance, influence sex characteristics and promote illness and injury prevention [70–72]. Moreover, Birukawa found that steroid hormones are involved in the osmoregulation of cetaceans [73]. Phosphonates are effective chelating agents that remain stable under harsh conditions. Phosphonates are also regularly used in reverse-osmosis systems [74]. However, the interaction or correlation between osmoregulation and phosphonates is still not clear and requires further study in aquatic animals.

Juvenile white shrimp (1.98 ± 0.28 g) were obtained from the Shenzhen base of the South China Sea Fisheries Research Institute, Shenzhen, China, and they were stocked in six tanks at a density of 40 shrimp per tank (500 L) at 17 psu salinity for one week. The shrimp were acclimated to 3 psu and 30 psu in three tanks through daily increments of 2 psu prior to the start of the 8-week experiment. During the acclimation and experimental periods, the shrimp were fed a commercial diet (10% moisture, 40% crude protein, 8% crude lipid, 12% ash, 30% carbohydrates, 16.7 kJ g^-1 digestible energy) three times daily at 08:00, 16:00 and 22:00 h. Based on the amount of feed left over from the previous day, daily rations were adjusted to a feeding level slightly over satiation. The unused feed was removed daily with a siphon tube. The photoperiod was 12 h light and 12 h dark. Seawater was pumped from the Dayawan Coast in Shenzhen and filtered through an activated carbon cartridge for at least 3 d before entering the culture system. The tap water was aerated before it was added to the tank to adjust the salinity level. During the experiment, water equal to 1/3 of the tank volume was exchanged once daily. The water quality parameters were monitored 2–3 times a week throughout the feeding trial and maintained at pH 7.5–7.9, temperature 26–28°C, dissolved oxygen 4.8–6.4 mg L^-1, and total ammonia nitrogen < 0.02 mg L^-1 during the trial.

At the end of the experiment, shrimp were fasted for 24 h before sampling. Five shrimp at the intermolt stage C in each tank were dissected to obtain the hepatopancreas tissue for the transcriptome analysis.

Total RNA was extracted from the hepatopancreas using the TRIzol reagent (Invitrogen) according to the manufacturer’s instructions, and genomic DNA was removed using DNase I (TaKaRa). RNA quality was determined by a 2100 Bioanalyzer (Agilent), and the RNA was quantified using the NanoDrop 2000 spectrophotometer (ND-2000, Gene Company Ltd.). Only high-quality RNA samples (OD_260/280 = 1.8–2.2, OD_260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, > 10 μg) were used to construct the sequencing library.

The RNA-seq transcriptome library was prepared following the TruSeq RNA (Illumina, San Diego, CA) sample preparation instructions using 5 μg total RNA. Briefly, messenger RNA was isolated according to the polyA selection method by oligo beads [75] and then segmented (100 bp to 400 bp) by a fragmentation buffer. Next, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen) with random hexamer primers (Illumina). Then, the synthesized cDNA was subjected to end repair, phosphorylation and ‘A’ base addition according to Illumina’s library construction protocol. The libraries were size-selected for cDNA target fragments of 200–300 bp on 2% low-range ultra-agarose followed by PCR amplification using Phusion DNA polymerase (New England Biolabs) for 15 PCR cycles. After quantification by a TBS-380 fluorometer, the paired-end RNA-seq library was sequenced with the Illumina HiSeq 2000 system (2 × 100 bp read length). Raw reads were archived at the National Center for Biotechnology Information (NCBI) Sequence Read Archive under the accession No. SRP048814.

The raw paired-end reads were trimmed and quality controlled by SeqPrep (https://github.com/jstjohn/SeqPrep↗) and Sickle (https://github.com/najoshi/sickle↗) with default parameters. Clean data from the samples were then used to perform RNA de novo assembly with the program Trinity (http://trinityrnaseq.sourceforge.net/↗) [76]. All of the assembled transcripts were searched against the NCBI protein NR, STRING and KEGG databases using BLASTx to identify the proteins that had the highest sequence similarity with the given transcripts to retrieve their function annotations, and a typical E-value cut-off was set at < 1.0×10⁻⁵. BLAST2GO (http://www.blast2 go.com/b2 ghome↗) [77] was used to obtain the GO annotations of uniquely assembled transcripts for describing biological processes, molecular functions and cellular components. A metabolic pathway analysis was performed using KEGG databases (http://www.genome.jp/kegg/↗) [78, 79].

To identify differentially expressed genes (DEGs) between two different samples, the expression level of each transcript was calculated according to the fragments per kilobase of exon per million mapped reads (FRKM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/↗) [80] was used to quantify gene and isoform abundances, and the R statistical package software EdgeR (empirical analysis of digital gene expression in R, http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html↗) [81] was used for the differential expression analysis. A functional-enrichment analysis that KEGG was performed to identify the DEGs that were significantly enriched in GO terms and metabolic pathways (at a Bonferroni-corrected P-value 0.05) relative to the whole-transcriptome background. The GO functional enrichment and KEGG pathway analysis were conducted by Goatools (https://github.com/tanghaibao/Goatools↗) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do↗) [82].

Twenty randomly selected genes with significant expression in KEGG pathways were validated with quantitative real-time PCR (qPCR) and gene-specific primers designed using Primer Premier 6 (Table 7). Total RNA was extracted from the target hepatopancreas tissues using a TRIpure Reagent kit (Aidlab, RN01). Samples of polyadenylated RNA were reverse-transcribed using a TaKaRa kit (Cat. No. RR036A). The reactions were conducted in a total volume of 20 μl with the following reaction components: 2 μl 5X PrimeScript RT Master Mix (Perfect Real Time), 1 μg total RNA, and RNase-free dH₂O up to 20 μl. The protocol for reverse transcription was 37°C for 15 min and 85°C for 5 s. The qPCR analysis was conducted in the CFX96 Real-Time PCR system (Bio-Rad Laboratories, Richmond, CA) using the Ultra SYBR Mixture (WCBIO, CW0957). The amplifications were performed in a 96-well plate in a 20 μl reaction volume containing 10 μl Ultra SYBR Mixture (WCBIO, CW0957), 0.4 μl each gene-specific forward and reverse primer, 8.4 μl RNase-free water and 0.8 μl cDNA. The thermal profile for the Ultra SYBR Mixture PCR analysis was 95°C for 10 min followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. The β-actin gene was used as the reference, and each gene had three replicated wells. Relative fold changes were calculated in the Relative Expression Software Tool version 2009 based on the cycle threshold values generated by qPCR [83].

All relevant data are within the paper and its Supporting Information files.