Analysis of the transcriptome data in Litopenaeus vannamei reveals the immune basis and predicts the hub regulation-genes in response to high-pH stress

Soil salinization erodes the farmlands and poses a serious threat to human life, reuse of the saline-alkali lands as cultivated resources becomes increasingly prominent. Pacific white shrimp (Litopenaeus vannamei) is an important farmed aquatic species for the development and utilization of the saline-alkali areas. However, little is known about the adaptation mechanism of this species in terms of high-pH stress. In the present study, a transcriptome analysis on the gill tissues of L. vannamei in response to high-pH stress (pH 9.3 ± 0.1) was conducted. After analyzing, the cyclic nucleotide gated channel-Ca2+ (CNGC-Ca2+) and patched 1 (Ptc1) were detected as the majority annotated components in the cAMP signaling pathway (KO04024), indicating that the CNGC-Ca2+ and Ptc1 might be the candidate components for transducing and maintaining the high-pH stress signals, respectively. The immunoglobulin superfamily (IgSF), heat shock protein (HSP), glutathione s-transferase (GST), prophenoloxidase/phenoloxidase (proPO/PO), superoxide dismutase (SOD), anti-lipopolysaccharide factor (ALF) and lipoprotein were discovered as the major transcribed immune factors in response to high-pH stress. To further detect hub regulation-genes, protein-protein interaction (PPI) networks were constructed; the genes/proteins “Polymerase (RNA) II (DNA directed) polypeptide A” (POLR2A), “Histone acetyltransferase p300” (EP300) and “Heat shock 70kDa protein 8” (HSPA8) were suggested as the top three hub regulation-genes in response to acute high-pH stress; the genes/proteins “Heat shock 70kDa protein 4” (HSPA4), “FBJ murine osteosarcoma viral oncogene homolog” (FOS) and “Nucleoporin 54kDa” (NUP54) were proposed as the top three hub regulation-genes involved in adapting endurance high-pH stress; the protein-interactions of “EP300-HSPA8” and “HSPA4-NUP54” were detected as the most important biological interactions in response to the high-pH stress; and the HSP70 family genes might play essential roles in the adaptation of the high-pH stress environment in L. vannamei. These findings provide the first insight into the molecular and immune basis of L. vannamei in terms of high-pH environments, and the construction of a PPI network might improve our understanding in revealing the hub regulation-genes in response to abiotic stress in shrimp species and might be beneficial for further studies.


Introduction
Generally, hundreds of candidate differentially expressed genes (DEGs) are unavoidably obtained after filtering and analyzing the RNA-Seq data in non-model species or in organisms that lack whole-genome sequencing [15,17,35]. The protein-protein interactions (PPIs), which are the central relationship in most biological activities [36], are available for revealing the associations among genes [37]. Important genes and related pathways can be identified by constructing PPI networks, as can be observed in the report of Zang et al. (2018), four predominant genes in uterine leiomyosarcoma (uLMS) are identified after analyzing the PPI networks of 95 differentially expressed genes [38]. The "Betweenness" of an actor i is defined as the total number of shortest paths between pairs of actors that pass through i and is an indicator of the most influential factors in the network [39][40]. The algorithm of "Betweenness" is an effective method for analyzing the candidate hub genes in PPI networks. Karbalaei et al. (2018) analyzed Alzheimeir's disease (AD) and non-alcoholic fatty liver disease (NAFLD) using PPI networks, and, based on the Betweenness values, seven hub-bottleneck node genes were discovered [41].
The balance of acid-alkali ions in water is of essential importance for aquatic crustaceans [42][43][44][45][46]. The decreasing pH value in seawaters caused by ocean acidification (OA) [47] negatively affects the mortality, growth, and reproduction of marine crustaceans [48][49][50]. Kurihara et al. (2008) evaluated that the survival rates of Palaemon pacificus were largely reduced (55%) after long-term (30 weeks) exposing to low-pH environment (pH 7.89 ± 0.05) [51]. Gao et al. (2018) reported that composition of the fatty acids (monounsaturated fatty acids, polyunsaturated fatty acids, eicosapentaenoic acid and docosahexaenoic acid), which were well known to be involved in immune modulatory effects, were significantly affected by low-pH environment (pH 7.6-7.8) in brine shrimp Artemia sinica [52]. The immunity of L. vannamei has been demonstrated to be influenced after being exposed to high-pH stress [53][54]. The report of Li and Chen (2008) showed that several immune parameters including phenoloxidase (PO) activity, superoxide dismutase (SOD) activity and total haemocyte count (THC) were significantly decreased after transferring the shrimp to high-pH stress environments, which led to decrease the resistance of L. vannamei against Vibrio alginolyticus [53]. When L. vannamei were exposed to acute stress, with pH 5.6 and pH 9.3, oxidative damage occurred by increasing the reactive oxygen species (ROS) production and the DNA comet assay value [54]. While Sookruksawong et al. (2013) proved that the PO pathway genes were also associated with resistance to Taura syndrome virus (TSV) in L. vannamei [55], Zhang et al. (2007) elucidated that the expression of SOD gene changed rapidly and dynamically in response to white spot syndrome virus (WSSV) in shrimp Fenneropenaeus chinensis [56]. These results indicated that the farming of L. vannamei in high-pH environments might weaken the immunity of shrimp to pathogens, and probably affect the final production.
Molecular signaling pathways can help to elucidate the activity and coordination of the cellexpression network in response to extracellular stresses [17], which is concerned to play key roles in revealing the regulation pattern of innate immune system in shrimp [57][58]. The important types of signaling genes, such as calcium/calmodulin-dependent protein kinases [59], mitogen-activated protein kinases (MAPKs), transcription factor (TF) and serine/threonine-protein kinases [60][61], are recognized as the key components for converting extracellular stimuli into a wide range of cellular responses in organisms, from yeast to humans [4,27,[59][60][61]. The understanding of the regulatory features of the immune factors and signaling genes involved in high-pH environmental adaptation in L. vannamei will be beneficial and necessary in further studies.
In the present study, high-pH RNA-Seq was carried out in the gill tissues of L. vannamei and the major types of signaling genes, regulation pathways and immune factors were investigated. The PPI networks of the major differentially expressed genes were constructed, and the candidate hub regulation-genes were predicted. Our results will be helpful for revealing a deeper insight into the molecular basis for the adaptation of L. vannamei to a high-pH environment.

Treatments of the shrimp samples
According to Wang et al. (2009), oxidative damage (ROS production and DNA damage increased in the haemocytes and the hepatopancreas cells) occurred when the shrimps were exposed in pH 9.3 stress environment [54]. In this study, healthy reared shrimps (supplied by Donghai Island Department, Yuehai Feed Group co., LTD, Zhanjiang, China) with an average weight of 11.09 ± 2.37 g were selected, three groups were placed in control seawater (CT group, pH 8.0 ± 0.1, 3.0% salinity, 28˚C temperature), in pH 9.3 ± 0.1 seawater for 1 hour (T_1 group, 3.0% salinity, 28˚C temperature) and in pH 9.3 ± 0.1 seawater for 48 hours (T_48 group, 3.0% salinity, 28˚C temperature) (Fig 1A), each group contained 20 shrimps and cultured in the same separate tank (50 L of seawater). The pH environments in T_1 and T_48 group were raised and maintained to 9.3 ± 0.1 by the gradual addition of 1 mol/L each of Na 2 CO 3 and NaHCO 3 in the original tank, pH was monitored over the subsequent 1 hour and 48 hours using a 228573-A01 pH electrode (Orion Research Inc., Beverly, MA, USA) attached to an Orion320 pH meter [54].
The gill tissues were separated and were immediately stored in liquid nitrogen. Six biological replicates were obtained from each group, 50-100 mg of gill tissues were extracted by 1 mL TRIzol Reagent (Invitrogen, USA) according to the manufacturer's instructions for each individual, 1 μg RNA of each individual were gained, 6 μg RNA from the six biological replicates were pooled in each group and were used for the subsequent experiments [17]. All the experiments were conducted in South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou, China.

RNA isolation and cDNA library construction
Three separate libraries (CT, T_1 and T_48) were sampled. The following steps were described by   [17]. First, the poly-A tailed mRNA was concentrated by oligo (dT)coated magnetic beads; fragmentation was performed by incubation in an NEB Next First Strand Synthesis Reaction Buffer (NEB, USA); the second strand was generated with the buffer, dNTPs, RNase H and DNA polymerase-I; and the adapters were ligated to the synthesized cDNA fragments after an end repair step. The Illumina HiSeq 4000 sequencing platform (paired-end 150 bp reads were generated) (BGI Company, Shenzhen, China) was employed to construct the libraries (CT, T_1 and T_48).
The detailed assembly methods of the transcripts were described in the reports of Li et al. (2017) [62] and   [63], briefly: the raw sequence data were filtered with SOAPdenovo software [64][65], adaptor-polluted, low-quality (A Phred quality score is used to identify the quality of the reads) [66][67] and high content of unknown base (N) reads were removed; after filtering, Trinity software (parameters:-min_contig_length150-CPU 8-min_kmer_cov 3-min_glue3-bfly_opts '-V 5-edge-thr = 0.1-stderr') [68] was employed to perform de novo assembly with the clean reads, and these were first joined to Contigs and then connected to sequences. The sequences that could not be extended were defined as Unigenes. TIGR Gene Indices clustering tool (parameters: -l40-c 10-v 25-O '-repeat_stringency 0.95-minmatch35-minscore35') [69] was used to obtain the final Unigenes for further analysis [62].

Quantification of the differentially expressed genes (DEGs)
Two comparisons of T_1/CT and T_48/T_1 were applied among the three libraries (CT, T_1 and T_48). The DEGs in comparison between T_1/CT represented the transcripts in response to an acute stress pH environment. The DEGs in the comparison between T_48/T_1 represented the transcripts in response to the endurance pH environment (endured high-pH for 48 hours). The metric fragments per kilobase of transcript per million mapped reads (FPKM) value was used to estimate transcript abundance on the basis of eliminating the influence of different gene lengths and sequencing discrepancies [70]. Therefore, the FPKM values were directly used to compare the differences in gene expressions among the samples. The algorithm developed by Audic and Claverie (1997) was used to compare the differences in various gene expression levels [71]. A value of the FDR (false discovery rate) � 0.001 and an absolute value of the Log2 ratio � 1.0 were set as the criteria values. The coexisted DEGs in the two comparisons were filtered, the Log 10 FPKM values were used to represent the expression values of transcriptional libraries (CT, T_1 and T_48), and the software Multiple Experiment Viewer (MeV 4.9.0) [72] was employed to analyze the expression profiles of the large numbers of coexisted DEGs.

Analysis of the PPI network for the major DEGs
Major DEGs from the two comparisons were obtained by a P-value < 0.05, -Log 10 FDR � 10.0 and an absolute value of Log 2 FoldChange � 3.0. The PPI networks of the major DEGs were constructed as Szklarczyk et al. (2011) [91] and Yang (2018) [92] described previously. Simply, the major DEGs were annotated by the SwissProt database, the annotated proteins were submitted to the STRING database (http://string-db.org/) for searching the experimental and interaction information, the detected proteins and interactions were exported and conserved as a file of simple tabular text output (tab separated values, TSV), and the criteria of the combined score was set to no less than 0.6 to keep the network analysis of PPI manageable [93]. The PPI network was visualized with Cytoscape version 3.4.0 software (http://cytoscape.org/) [94], the protein nodes were attributed to forming the circle layout, and the detailed user manual of the software was listed in the website http://www.gnu.org/licenses/lgpl-2.1.html. To analyze the candidate hub genes and interactions, the algorithm of "EdgeBetweenness" was employed to calculate the PPI networks [39][40], and the EdgeBetweenness values were obtained with Cytoscape version 3.4.0 software [94].

Real-time PCR validation
To verify the expression profiles of genes in our RNA-Seq results, the represented 10 DEGs in the two comparisons were selected for validation of the Illumina sequences by real-time PCR analysis. All the primers used were listed in S1 Table. The treatments that were performed were described above, and the relative transcript abundance was obtained by normalizing with the expression of the L. vannamei β-actin gene based on the 2 -ΔΔCT method [17,95].

Overview of the RNA-Seq data
After filtering, the distribution of the base composition on clean reads was monitored, all the adenine bases (A) were detected to overlap with the thymine bases (T), and all the guanine bases (G) were observed to overlap with the cytosine bases (C) (Fig 1B). The total clean reads per library ranged from 44.04 Mb to 44.58 Mb, and the clean reads ratio ranged from 82.08% to 83.98% ( Table 1). The clean reads of the three libraries (CT, T_1 and T_48) were deposited at the US National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra) with the GenBank accession No. SRP162932. The clean reads were assembled and clustered into Unigenes, the total number of Unigenes was 78,080 (Table 2), the most abundant Unigene length were clustered in the group of 0-300 bp (Fig 1C), and the coexisted Unigene number for the three libraries was 36,589 ( Fig 1D). Overall, the Unigenes returning hits in Nr, KEGG, COG, Swiss-Prot and InterPro databases were 26,850, and the compositions for each database were 26,431 (98.44%), 22,771 (84.81%), 13,614 (50.70%), 23,241 (86.56%), and 13,656 (51.38%), respectively ( Fig 1E). The coexisted gene number of the five databases was 7,944 (which occupied 29.59% of all annotated genes) ( Fig 1E). The complete nucleotide sequences of each Unigene were available in S1 File. Table 1. Summary of quality metric reads after filtering.

Quantification review of the DEGs
A total of 2,832 up-and 1,952 down-regulated genes were detected in the comparison of T_1/ CT (Fig 2A), and 2,276 up-and 3,836 down-regulated genes were discovered in the comparison of T_48/T_1 (Fig 2A). A total of 2,440 (63 + 94 +2283) coexisted DEGs (occupied 28.86% of all DEGs) were detected, 63 and 94 genes were observed to be up-regulated or down-regulated, respectively, in both T_1/CT and T_48/T_1 comparisons ( Fig 2B). To investigate the expression profiles, the 2,440 coexisted DEGs were clustered and classified into 7 gene clusters (Hierachical clustering was analyzed by Euclidean distance, the gene tree was classified by the distance threshold ranged from 2.30-2.50) (Fig 2C and 2D), of which, the top three clusters were detected (Cluster 7, Cluster 4 and Cluster 6), and the gene numbers of the top three clusters contained 82.54% of the total coexisted DEGs (Fig 2). The expression profiles of the top three cluster genes were profoundly influenced in the T_1 group and recovered to nearly the control level thereafter in the T_48 group (Fig 2D), which indicated that the expression level of the majority of the regulated genes in shrimp returned to normal conditions after long-term exposure to the pH 9.3 environment.

GO and COG analysis of the DEGs
After analyzing, the biological process was detected as the most annotated GO category in both comparisons of T_1/CT (329 DEGs) and T_48/T_1 (448 DEGs) ( Table 3). The top 15 clustered GO terms in categories of biological process, cellular component and molecular function were summarized in both the two comparison groups, respectively (Fig 3). Take the category of biological process for example, cellular process (GO:0009987), metabolic process (GO:0008152), single-organism process (GO:0008150), organic substance metabolic process (GO:0071704) and cellular metabolic process (GO:0044237) were discovered as the top 5 annotated GO terms for both libraries comparisons (Fig 3). Specifically, 241 DEGs were annotated to GO term of cellular process in the comparison group of T_1/CT and 307 genes were annotated to the same GO term in the comparison group of T_48/T_1 (Fig 3). The membrane (GO:0016020) was the most annotated GO term of the cellular component category in the T_1/CT comparison group, and the cell (GO:0005623) was the most annotated GO term in the T_48/T_1 comparison group (Fig 3). The catalytic activity (GO:0003824) was the most annotated GO term of the molecular function category in both the T_1/CT and T_48/T_1 comparison groups (Fig 3).
The COG data of the DEGs were summarized, the top three functional categories of "General function prediction only" (657 DEGs were clustered), "Translation, ribosomal structure and biogenesis" (588 DEGs were clustered) and "Carbohydrate transport and metabolism" (337 DEGs were clustered) were detected (Table 4), and the detailed DEGs information for COG annotation was available in S2 Table.

Identification of the major KEGG pathways
The signal transduction pathway was the most annotated Level-2 KEGG term in both the two comparisons (T_1/CT and T_48/T_1), 640 DEGs were annotated to signal transduction in the comparison group of T_1/CT, and 741 were annotated in comparison of T_48/T_1 (Fig 4). The cAMP signaling (KO04024), cGMP-PKG signaling (KO04022) and PI3K-Akt signaling pathways (KO04151) were the top three pathways with more KOs identified in their composition regarding the two comparison groups (S3 Table). The KEGG map and DEGs of the cAMP signaling pathway were traced. The majority DEGs in the cAMP signaling pathway were annotated to the components of cyclic nucleotide gated channel beta 1 (CNGCB1) (K04952) and patched 1 (Ptc1) (K06225) in both the T_1/CT and T_48/T_1 comparisons (S4 Table), indicated that the components of CNGCB1 and Ptc1 might be the major two KOs of the annotated components in cAMP signaling pathway (Fig 5; S4 Table).    Table).
The transcription factors (annotated gene number: 70), serine/threonine-protein kinase signaling genes (annotated gene number: 16) and calcium signaling genes (annotated gene number: 11) were selected from the total annotated genes (S6 Table), while the MAPK genes were not detected in the results, indicating that the transcription of the MAPK signaling genes  Table).

PPI network of the major DEGs
The volcano-plot figures of the DEGs were listed on the left side of Fig 8A and 8B; after filtering, 536 of the major DEGs were detected in the comparison group of T_1/CT, and 818 were detected in the comparison group of T_48/T_1 (S7 Table). Base on the STRING database, 175 protein pairs were obtained in the comparison group of T_1/CT (S8 Table), and 164 protein pairs were detected in the comparison group of T_48/T_1 (S9 Table). The main PPI network of T_1/CT was constructed; the "Heat shock 70kDa protein 8" (HSPA8), "Polymerase (RNA) II (DNA directed) polypeptide A" (POLR2A), "DEAH (Asp-Glu-Ala-His) box polypeptide 9" (DHX9), "Serine/arginine repetitive matrix 1" (SRRM1) and "Splicing factor 3a, subunit 2" (SF3A2) were detected as the top five proteins having the quantity of counted protein-interactions ( Fig 8A, the right side). The main PPI network of T_48/ T_1 was also formed; the HSPA8, SRRM1, "Serine/arginine-rich splicing factor 1" (SRSF1), "Serine/arginine repetitive matrix 2" SRRM2 and "Serine/arginine-rich splicing factor 4" (SRSF4) were discovered as the top five genes/proteins containing the number of counted protein-interaction ( Fig 8B, the right side).
With the help of Cytoscape software, the top EdgeBetweenness values of the PPI networks were summarized and listed in Table 5, and the interaction of "EP300 (interacts with) HSPA8" contained the highest EdgeBetweenness value in the comparison group of T_1/CT and showed the hub protein-interaction in the network ( Table 5). The interaction of "HSPA4 (interacts with) NUP54" contained the highest EdgeBetweenness value in the comparison group of T_48/T_1 and revealed the hub protein-interaction in the network ( Table 5). The results indicated the important protein-interactions involved in the high-pH adaptation of L. vannamei.
Based on Table 5 and on the PPI networks, the candidate hub proteins were summarized and displayed in Table 6. The POLR2A gene was discovered in six pairs of important hub protein-interactions (detailed information of the protein-interactions was listed in Table 5) in the T_1/CT comparison group, which indicated that this gene was one of the most important regulated genes in response to high-pH stress in L. vannamei (Tables 5 and 6; Fig 8A and 8B, the right side). The genes POLR2A, EP300, HSPA8, SMARCA5, DDX17, DHX9 and RTF1 were discovered as the hub regulation-genes in the comparison group of T_1/CT; the genes HSPA4, FOS, NUP54, SRSF1 and CYCS were detected as the hub regulation-genes in the comparison group of T_48/T_1 (Tables 5 and 6).

Real-time PCR validation
The treatment groups of CT, T_1 and T_48 were verified, and 10 genes were used to test the validations of the RNA-Seq results. The real-time PCR outcomes were summarized, and the Litopenaeus vannamei immune basis and the hub regulation-genes in response to high-pH stress grouped comparison results showed that all 10 candidate genes in the qPCR verification were consistent with the results of RNA-Seq technology (Fig 9).

Discussion
The efficient of the transcriptome was verified in various studies [96]. According to the studies on L. vannamei, RNA-Seq had been applied to reveal several types of stress adaptation mechanisms [15,17]. In the present study, the catalytic activity (GO:0003824) was the most annotated GO term and the signal transduction pathway was the most annotated Level-2 KEGG term in response to high-pH stress. cAMP is one of the most versatile second messengers, which mediates diverse cellular responses to extracellular signals [97]. Intracellular cAMP are induced by soluble adenylyl cyclases (sACs) and mainly effect the protein kinase A (PKA) for downstream regulation [97]. In the present study, cAMP signaling pathway was identified as the most annotated terms of signal transduction in KEGG pathway. The important role of cAMP in response to other types of stresses had been demonstrated in previous studies. The report of Zhao et al. (2016) proved that the cAMP concentrations of L. vannamei increased significantly and reached the maximum at 12 hours in low salinity (1.6%) environment [98].  [99]. The cyclic nucleotide gated channel (CNGC), a type of component in the cAMP signaling pathway, was modulated by Ca 2+ -binding proteins and directly bound to cAMP to transduce the cellular signals for downstream regulation [100][101]. Hedgehog (Hh) signaling pathway played a central role in promoting the cellular processes and activities [102], and the aberrant activation of Hh signaling led to various types of cancer and defects in the organism cells [103]. The Ptc1 protein of cAMP signaling pathway could maintain the Hh pathway in an inactive state by inhibiting translocation of the signaling effector Smoothened (Smo) [104][105]. In this paper, the CNGCB1 (K04952) and Ptc1 (K06225) were detected as the major two parts KOs of the annotated components in cAMP signaling pathway. The results indicated that the components of CNGC might be involved in transducing the cellular signal of high-pH stress and Ptc1 might be participated to maintain the Hh pathway in a proper state for cell survival in high-pH environment of L. vannamei.
Transcription factors and kinases were two types of frequently traced categories involved in stress response studies [27, [106][107]. The kinase genes, mainly comprising MAPKs, Ser/Thr protein kinases and calcium-related signaling genes, were associated with the activity and coordination of the cell-expression network in response to extracellular stresses [17]. The signaling genes of transcription factors, Ser/Thr protein kinases, Ca 2+ signaling pathway and MAPK cascades genes were summarized in the transcriptomic analysis of Chrysanthemum nankingense involved in low temperature tolerance [27]. In this study, the major signaling genes of transcription factors, Ser/Thr protein kinases and calcium-related signaling genes were detected; the MAPKs were not significantly transcribed in the adaptation of high-pH environment. Calcium impacted nearly every aspect of cellular life and fulfilled its key role in transferring extracellular stress signals into chromatin for regulation [17,59,108]. Ca 2+ concentration was determined by feedback and feed-forward regulation of the involved transport proteins [109]. The CNGCs were nonselective cation channels that did not discriminate well between alkali ions and even pass divalent cations, in particular Ca 2+ [100]. In the present Litopenaeus vannamei immune basis and the hub regulation-genes in response to high-pH stress study, the CNGC was suggested as the major annotated component in cAMP signaling pathway and the calcium-related signaling genes were detected as important type of signaling genes. Therefore, the results indicated that the molecular components of "CNGC-Ca 2+ " might play key role in transducing the high-pH stress signal in L. vannamei. . On the right side of this figure: the nodes represented the genes/proteins, the edges represented the proteininteractions, the larger and brighter the node was, the more count-numbers of proteins that interact with this node, and the thicker and brighter the edge was, the higher the value of Edgebetweenness in this edge. https://doi.org/10.1371/journal.pone.0207771.g008 Crustaceans mainly depend on a non-specific immune system that includes circulating hemocytes and various active substances [98,110]. Environmental stress might increase the vulnerability of shrimp to the pathogens normally present in the water by reducing the capacity of their immune responses [53, 80,111]. In this paper, 17 types of major immune-related categories of DEGs were queried; after analyzing, the IgSF, HSP, GST, proPO/PO, SOD, ALF and lipoprotein were detected as the top 7 categories of expressing immune factors in response to high-pH stress in L. vannamei. The IgSFs were defined as molecules containing homologous "a" indicated the experimentally determined interaction scores obtained from the STRING database; "b" represented the manageable scores for this pair of protein-interactions; "c" showed the EdgeBetweenness value obtained from the Cytoscape software. domains and thought to play significant roles in innate defense in invertebrates [79,112]. The IgG-like protein in L. vannamei was identified to contain the conserved domain with hemocyanin [112], and was involved in oxygen transporting and participating the antiviral and phenoloxidase-like activities [113], which indicated that the IgSFs might be one of the important   Litopenaeus vannamei immune basis and the hub regulation-genes in response to high-pH stress types of immune factors for high-pH adaptation in L. vannamei. The GSTs were well known to alleviate the toxicity of diverse endogenous and exogenous compounds (xenobiotics) [114] through peroxidase and isomerase activities, biosynthesis of physiologically important biomolecules, as well as through the ligand binding (non-catalytic activity) [115][116][117], and possibly were involved in response to high-pH stress in L. vannamei. Our results indicated the major transcribed immune factors of L. vannamei in term of high-pH stress and the results would be benefit for future studies referring to the farming industry of this species in saline-alkali water areas.
The protein-protein interactions (PPIs) are the available bioinformatics analysis for computing the associations between genes and for predicting the key candidate genes and pathways [36,118]. Yang et al. (2018) identified 3 candidate genes from 440 up-and 256 downregulated DEGs as the potential hub regulation-genes involved in the acute phase of myocardial infarction [92]. In this paper, the PPI networks were constructed with 536 (T_1/CT comparison) and 818 (T_48/T_1 comparison) major DEGs, and the protein-interactions of "EP300-HSPA8" and "HSPA4-NUP54" were detected as the most important interactions in response to the high-pH stress. The important biological meanings of stress-induced proteininteractions raised by the present study had partly been interpreted in previous works. The acetyltransferase, EP300, was identified to control the cellular level of activatable heat shock transcription factor 1 (HSF1) in a manner linked with the clearance of misfolded proteins for adaptation to proteotoxic stress in HeLa cells [119]. Heat shock proteins could help to escape the mRNA export inhibition via recruiting nucleoporin-related proteins, and this biological interaction during stress conditions [120][121]. The results indicated that the PPIs analyzing might be available for computing and predicting the hub genes of L. vannamei in adapting to high-pH stress.
Heat shock proteins (HSPs) are defined as molecular chaperones, which play vital roles in maintaining protein homeostasis and in facilitating protein folding and degradation [122]. According to their monomeric molecular mass, HSPs are divided into HSP100, HSP90, HSP70, HSP60, HSP40 and small HSPs (sHSPs) [81,122]. HSP70s are implicated in a wide variety of cellular processes via binding their substrate-binding domain to a short polypeptide with recognition motifs [123]. In the present paper, high-pH stress experiments were conducted on L. vannamei, the HSPs were detected as the main expressed genes among the major immune factors, and two HSP70 family genes (HSPA8 and HSPA4) were discovered as the hub regulation genes in the PPI networks. The key role of HSP70s had been proved involving in response to various stressful stimuli in previous studies [124][125]. The results of Duan et al. (2018) showed that the expression of HSP70 increased to the highest level by exposing the L. vannamei to nitrite stress for 6 hours [126]. Western blot analysis in the gill tissue of L. vannamei indicated that the expression of HSP70 was up-regulated by infecting the pathogens of white spot syndrome (WSSV) and hypodermal-hematopoietic necrosis (IHHNV) [127]. The mRNA and protein expression levels of HSP70 were observed highest in the heat-stress group of goats (Capra hircus) [128]. These results indicated that HSP70s might play essential important roles in the adaptation of high-pH stress in L. vannamei.
In conclusion, the first transcriptome analysis involved in high-pH stress was conducted on gill tissues of L. vannamei, the major high-pH associated immune factors and pathways were identified, PPI networks of major DEGs were constructed, and the hub regulation-genes were predicted. Our results provided the expression profiles and potential candidate hub genes for revealing the regulation mechanism in L. vannamei in response to high-pH stress. The construction of a PPI network was a useful method for further understanding and for revealing the hub regulation-genes in response to abiotic stress in the shrimp species.
Supporting information S1 File. The complete nucleotide sequences of each Unigene. (TXT) S1