The Hsp70 Gene Family in Boleophthalmus pectinirostris: Genome-Wide Identification and Expression Analysis under High Ammonia Stress

Simple Summary Heat shock proteins 70 is a family of proteins, which were expressed in response to a wide range of biotic and abiotic stressors. The development of genomic resources and transcriptome sequences makes it practical to conduct a systematic analysis of these genes. In this study, exhaustive searches of all genomic resources for Boleophthalmus pectinirostris Hsp70 genes were performed and their responses to high environmental ammonia stress were investigated. Besides, selection test was implemented on those duplicated genes, and the phylogenetic tree, gene structure, and motif analysis were also constructed to assign names of them. The result showed that there were 20 Hsp70 genes within the genome of Boleophthalmus pectinirostris, and some sites in the duplicated genes may experience positive selection, and most of Hsp70 genes were downregulated after exposure to high concentration ammonia. The present results of this study can be used as a reference for further biological studies on mudskippers. Abstract Heat shock proteins 70 have triggered a remarkable large body of research in various fishes; however, no genome-wide identification and expression analysis has been performed on the Hsp70 gene family of Boleophthalmus pectinirostris. In this study, we identified 20 Hsp70 genes within the genome of B. pectinirostris and provided insights into their response to high environmental ammonia (HEA) stress. Positive selection on stress response genes and expansion of hspa1a and hspa1a-like genes might be related to terrestrial adaptations in this species. The expression patterns of the Hsp70 gene family in the gill and liver of B. pectinirostris under HEA stress were studied by examining transcriptome data. The results showed that most Hsp70 genes were downregulated after high concentration ammonia exposure. The downregulation may be related to the hypoxic condition of the tissues.


Introduction
Heat shock proteins (Hsps) are a large group of molecular chaperones. The expressions of this gene family are induced by a wide range of biotic and abiotic stressors, some members are constitutively expressed in non-stressed cells served as housekeeping proteins [1]. The high degree of identity on the amino acid sequence and their essential physiological roles in nearly all organisms, make this group of proteins unique [2]. Based on their molecular weight, Hsps are classified into Hsp110, Hsp90, Hsp70, Hsp60, Hsp40, Hsp10, and small Hsps [3]. Among them, the Hsp70 gene family is the most extensively studied group of HSPs.

Genome-Wide Identification of Hsp70 Genes in B. pectinirostris
A total of 23 putative Hsp70 genes were initially obtained from BLAST and HMM searches. Based on the confirmation of Pfam and NCBI CDD scans, three candidate genes (without the Hsp70 domain) were discarded. Twenty residual members symbolizing the unique Hsp70 gene family in B. pectinirostris were used to create robust nomenclature following guidelines for the nomenclature of the human heat shock proteins [36]. Detailed information about members of Hsp70 gene family is shown in Table 1. To avoid confusion created by the Hsp70 names and to clarify how the names were denoted to B. pectinirostris Hsp70 genes, comparison the nomenclature among human, zebrafish, and B. pectinirostris was listed in Table 2.
Among these genes, full lengths of all coding sequences were detected in both transcriptome and genome databases. Three Hsp70 genes (XM_020937593.1; XM_020934606.1; XM_020934610.1) were identified as 'heat shock cognate 71 protein' before, and were renamed as hspa8, hspa8a.1, and hspa8a.2 after the phylogenetic tree, respectively. The hspa14 had the shortest conserved domain with 375 amino acids, whereas the longest domain was found in the hspa1a and hspa1a-like genes.
Many Hsp genes are bound by hsf or hsf-like genes and could be induced in an Hsf-dependent manner upon heat shock response. However, only hspa9 gene was found to be bound by hsf-like gene in the Hsp70 gene family of B. pectinirostris. In other words, the hspa9 gene was predicted to be Hsf-induced in B. pectinirostris.

Phylogenetic Relationships of the Hsp70 Genes among Species
An unrooted phylogenetic NJ tree was generated using the amino acid sequences of Hsp70s ( Figure 1). Names were assigned to each of them based on the clade of the NJ tree. Seven copies of B. pectinirostris hspa1a genes (hspa1a, hspa1a-like) were found to be highly orthologous to the medaka hsp70 gene, the Nile tilapia hsp70 gene and the zebrafish hsp70 and hsp70-like genes. Based on previous studies and the phylogenetic tree, the hspa1 gene in the ancestor of teleosts might be divided into two clades, and the hsp70 genes appeared earlier than hspa1b. As a result, the hsp70 genes were named hspa1a and hspa1a-like in B. pectinirostris following the guidelines for the nomenclature of the human heat shock proteins. B. pectinirostris apparently had more duplicated hspa1a-like genes than other teleosts as shown in Figure 1. The hspa8b gene was not found in this species, whereas two copies of hspa8a were present in the phylogenetic tree. All members of the B. pectinirostris Hsp70 gene family were well distributed into distinct groups, first clustered with corresponding genes of other fish species and supported by high bootstrap values ( Figure 1).
Most pairs of orthologs from B. pectinirostris and other fish species were presented, suggesting that the common ancestral genes of this gene family might have existed before the speciation of fish species. Several pairs of duplicated Hsp70 genes were found in B. pectinirostris, indicating that Hsp70 genes might undergo some duplication events after speciation.  were well distributed into distinct groups, first clustered with corresponding genes of other fish species and supported by high bootstrap values ( Figure 1). Most pairs of orthologs from B. pectinirostris and other fish species were presented, suggesting that the common ancestral genes of this gene family might have existed before the speciation of fish species. Several pairs of duplicated Hsp70 genes were found in B. pectinirostris, indicating that Hsp70 genes might undergo some duplication events after speciation.   Conserved motif analysis was performed based on the evolutionary relationships among the complete nucleotide sequences of the B. pectinirostris Hsp70s (Figure 3). Fifteen putative motifs were searched for in each Hsp70s as shown in Figure 3. After GOMo search, motifs 1, 2, 3, 9, 10, and 11 were annotated as sequence-specific DNA binding and transcription factor activity motif, and motifs 14 and 15 were annotated as K + and Mg 2+ potassium transport, respectively. As shown in Figure 3, all the identified Hsp70 genes contained motif 1, and most of them had motifs 3, 10, and 14, which might contribute to the identification of this gene family and understanding of their potential functions. Meanwhile, the Hsp70s from close evolutionary clusters shared similar motifs. The results of the motif analysis provided further support to the phylogenetic classification of Hsp70 gene family.

Expression Profiles of Hsp70 Genes in Different Tissues under HEA Stress
To study the expression regulation of Hsp70 genes in different tissues, transcriptome of gills and livers were analyzed ( Figure 4; Table 3). As shown in Table 3, 8 out of 20 Hsp70 genes were significantly involved in HEA stress responses (log 2 FC > 1.5 or < −1.5). Among these, three genes were downregulated and only one gene (hspa8) was upregulated in the liver, whereas seven genes were strongly downregulated in the gill. In addition, hspa1a.1, hspa1a.2, and hspa1al.4 genes were strongly downregulated in both tissues (log 2 FC: −4.14 to −1.75). The hspa1a.3 and hspa5 genes were downregulated only in the gill (log 2 FC: −2.21 and −1.61, respectively) after HEA treatment, whereas the hspa8 gene was upregulated only in the liver (log 2 FC: 2.92). Apparently, the duplicated hspa1a genes were more inducible than other genes (Table 3) after HEA stress. As a whole, the genes from the gill were more reactive than those from the liver (Table 3; Figure 4).

Expression Profiles of Hsp70 Genes in Different Tissues under HEA Stress
To study the expression regulation of Hsp70 genes in different tissues, transcriptome of gills and livers were analyzed (Figure 4; Table 3). As shown in Table 3, 8 out of 20 Hsp70 genes were significantly involved in HEA stress responses (log2FC > 1.5 or < −1.5). Among these, three genes were downregulated and only one gene (hspa8) was upregulated in the liver, whereas seven genes were strongly downregulated in the gill. In addition, hspa1a.1, hspa1a.2, and hspa1al.4 genes were strongly downregulated in both tissues (log2FC: −4.14 to −1.75). The hspa1a.3 and hspa5 genes were downregulated only in the gill (log2FC: −2.21 and −1.61, respectively) after HEA treatment, whereas the hspa8 gene was upregulated only in the liver (log2FC: 2.92). Apparently, the duplicated hspa1a genes were more inducible than other genes (Table 3) after HEA stress. As a whole, the genes from the gill were more reactive than those from the liver (Table 3; Figure 4).

Selection on Duplicated Hsp70 Genes
As mentioned above, B. pectinirostris had more duplicated hspa1a genes and most of them were significantly involved in HEA stress responses in both the gill and the liver. To better understand the species specific gene expansion, codon-based site models of evolution implemented in PAML Version 4.9 (Phylogenetic Analysis by Maximum Likelihood) were used. The results of comparisons of M2a vs. M1a and M8 vs. M7 showed that significance of positive selection existed among these genes ( Table 4). The Bayes empirical Bayes (BEB) can be used to identify sites under positive selection if the likelihood ratio test is significant [34]. The BEB showed that several sites were under positive selection in M2a and M8 (Table 4).

Discussion
In this study, we performed an overall analysis of the Hsp70 gene family in B. pectinirostris, including an analysis of phylogeny, gene structure, conserved motifs, expression patterns under HEA stress and selection tests. These information may be useful for genome analysis and annotation as well as for evolutionary studies in fish species.
A total of 20 Hsp70 genes were identified and annotated in this species. Compared with humans and other fish species, most of the Hsp70 gene members were found in B. pectinirostris except hsph1 and hspa8b. The absence of the hsph1 gene in B. pectinirostris was consistent with results obtained for catfish [1] and other fish species (not including zebrafish). The hspa8b gene was seemingly lost in this species, whereas two repeats of its paralogs hspa8a were found. However, it is uncertain whether the hsph1 and hspa8b genes are truly missing from the B. pectinirostris genome, although exhaustive searches of all genomic resources for this species were performed.
It was difficult to assign names to some Hsp70 genes solely by the clade of the phylogenetic tree. Nevertheless, combining the NJ tree with motif analysis, one can easily and accurately distinguish among these genes. This method has also been applied to plants [37,38] and humans [14]. Moreover, the gene structure and the type, order and number of the motif changes may reflect their specific functions that are not shared with other genes. Due to the requirements of specific capabilities, animals under certain environments might duplicate relative genes throughout their population history.
Gene duplication has been thought to play an important role in species adaptation and could provide raw genetic material for genes with new functions [3,39]. The present study suggests that local gene tandem duplication may be an important mechanism of Hsp70 gene expansion in B. pectinirostris since all of the expanded hspa1a genes exist as tandem gene clusters (Table 1). Tandem duplications of hspa1a (named hsp70 before) paralogs have previously been described in zebrafish [40], stickleback, and Tetraodon nigroviridis [41]. Nevertheless, the number of the hspa1a tandemly duplicated genes differs among species, indicating that these tandem duplications may have independent origins. Furthermore, the results of selection tests on these hspa1a and hspa1a-like genes have shown that some sites may experience positive selection, although this is rare (with a proportion of about 0.73%) ( Table 4), supporting that these genes undergo sub-functionalization or acquire new functions in the B. pectinirostris. As is known, B. pectinirostris is a typically amphibious teleost fish, and the water-to-land transition must lead to the emergence of high environmental stress. Previous studies have shown a significant increase on expression levels of hspa1a genes in other species under severe environment, such as medaka under high temperature [4], Monopterus cuchia under high ammonia stress [24], and Umbra limi under exposure to the air [42]. It is reasonable to assume that the stress response genes and expansion of hspa1a and hspa1a-like genes might have evolved terrestrial adaptations in the B. pectinirostris, which enabled them to spend a considerable part of their lives on land. Previous studies have shown three different expression patterns under thermal stress in mammalian HSP70 gene family [14,43]: (A) strictly heat-inducible HSP70; (B) cell-cycle-dependent and heat-inducible HSP70; and (C) constitutively expressed and less stress-dependent HSP70 genes. In this study, hspa1a.2, hspa1al.3, hspa1al.5, hspa1b, hspa12a, and hspa12b genes were expressed constitutively at very low levels, and hspa8a.2, hspa9, hsph2a, hsph2b, hsph4, hspa13, and hspa14 genes-all of which should be put into group C-were expressed constitutively but scarcely induced by HEA stress. Other Hsp70 genes-including hspa1a.1, hspa1al.1, hspa1al.2, hspa1al.4, hspa8, hspa8a.1, and hspa5, which were put into group B-were relatively highly expressed or showed significantly different expression levels between test and control. There was no most inducible gene response to HEA stress. Different expression patterns were also found in Hsp70 genes of channel catfish after bacterial infection [1]. Hence, the present study focused on group B genes to study the expression regulation of Hsp70 genes in different tissues of B. pectinirostris in response to HEA stress. Six Hsp70 genes of group B were significantly up-or downregulated after HEA stress, indicating their involvement in stress responses. Three genes were downregulated and only 1 was upregulated in the liver, whereas five genes were strongly downregulated and none were significantly upregulated in the gill. The genes from the gill were apparently more reactive than those from the liver. This may be because excretion of ammonia occurs mainly through the gills of fishes [27], and gills are well known as the first guard to react to unfavorable environmental conditions [44]. Moreover, some Hsp70 genes were expressed in only one tissue, showing a tissue-specific pattern.
Because the concentration of ammonia applied in this experiment only induced physiological stress but did not result in the death of individuals [45], and when we took the function of the Hsp70 genes into consideration, it was unexpected that most of them were downregulated in both tissues in this study. Previous studies have shown that B. pectinirostris can decrease the production rate of ammonia from amino acid catabolism in hepatocytes [46] to slow down the build-up of internal ammonia [47] under ammonia exposure. The protein digestion-related genes were downregulated, directly or indirectly, causing the expression of most Hsp70 genes to decline in the liver after HEA stress.
However, it was difficult to explain the expression pattern in the gills. One plausible explanation for downregulated expression was that a hypoxic condition might be generated via the exposure to elevated NH 3 levels [48], and the hypoxia of tissues might restrain the expression of Hsp70 genes. To some extent, this speculation was supported by studies examining the toxicology of ammonia nitrogen, which showed that the ammonia toxicity increased as the concentrations of dissolved oxygen (DO) decreased, whereas the resilience of aquatic organisms to environmental ammonia stress can be improved when the concentrations of DO increased [49,50]. To further verify the hypothesis, the transcriptome data of the brain in the large yellow croaker (Larimichthys crocea) under hypoxia stress was downloaded from NCBI (accession number: SRX541138, SRX541136, and SRX541132) and was re-analyzed. When the 48-h test group was compared with the 0-h group, almost all of the expression profiles of Hsp70 genes were downregulated in the brain of the large yellow croaker under hypoxia stress [51] (Figure S1), which also supported the hypothesis of the present study. The expression patterns of Hsp70 genes under HEA and hypoxia stresses showed that the functions of this gene family might require sufficient oxygen. The other possibility for the downregulated expression was that Hsp70 genes could not be triggered by high ammonia concentration in this species under the provided experimental condition. As many Hsp70 genes are induced primarily upon acute protein-damaging conditions [52], however, B. pectinirostris possesses a greater capacity to detoxify ammonia and high tolerance to environmental ammonia [28], which might make ammonia fail to induce acute protein-damaging stress.

Conclusions
In the present study, we identified 20 Hsp70 genes within the genome of B. pectinirostris and provided insights into their response to high environmental ammonia (HEA) stress. This study may contribute to illuminate the regulatory mechanism of the Hsp70 gene family in response to environmental stress. Certainly, due to the absence of data from other time periods (such as 12 h or 48 h), the deficiency of the expression profiles of Hsp70 genes under both HEA stress and high concentrations of DO, and a general lack of biological and physiological knowledge, further verifications are needed to predict the exact implications of the changes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-2615/9/2/36/s1; Figure S1: The Hsp70 expression in the brain of the large yellow croaker based on FPKM values under hypoxia stress. Brains were harvested from six fish at the 0, 12 and 48 h time points; Table S1: Species accession numbers of Hsp70 genes in the study.

Conflicts of Interest:
The authors declare no conflict of interest.