Transcriptomic response to three osmotic stresses in gills of hybrid tilapia (Oreochromis mossambicus female × O. urolepis hornorum male)

Background Osmotic stress is a widespread phenomenon in aquatic animal. The ability to cope with salinity stress and alkaline stress is quite important for the survival of aquatic species under natural conditions. Tilapia is an important commercial euryhaline fish species. What’s more tilapia is a good experimental material for osmotic stress regulation research, but the molecular regulation mechanism underlying different osmotic pressure of tilapia is still unexplored. Results To elucidate the osmoregulation strategy behind its hyper salinity, alkalinity and salinity-alkalinity stress of tilapia, the transcriptomes of gills in hybrid tilapia (Oreochromis mossambicus ♀ × O. urolepis hornorum ♂) under salinity stress (S: 25‰), alkalinity stress(A: 4‰) and salinity-alkalinity stress (SA: S: 15‰, A: 4‰) were sequenced using deep-sequencing platform Illumina/HiSeq-2000 and differential expression genes (DEGs) were identified. A total of 1958, 1472 and 1315 upregulated and 1824, 1940 and 1735 downregulated genes (P-value < 0.05) were identified in the salt stress, alkali stress and saline-alkali stress groups, respectively, compared with those in the control group. Furthermore, Kyoto Encyclopedia of Genes and Genomes pathway analyses were conducted in the significant different expression genes. In all significant DEGs, some of the typical genes involved in osmoregulation, including carbonic anhydrase (CA), calcium/calmodulin-dependent protein kinase (CaM kinase) II (CAMK2), aquaporin-1(AQP1), sodium bicarbonate cotransporter (SLC4A4/NBC1), chloride channel 2(CLCN2), sodium/potassium/chloride transporter (SLC12A2 / NKCC1) and other osmoregulation genes were also identified. RNA-seq results were validated with quantitative real-time PCR (qPCR), the 17 random selected genes showed a consistent direction in both RNA-Seq and qPCR analysis, demonstrated that the results of RNA-seq were reliable. Conclusions The present results would be helpful to elucidate the osmoregulation mechanism of aquatic animals adapting to saline-alkali challenge. This study provides a global overview of gene expression patterns and pathways that related to osmoregulation in hybrid tilapia, and could contribute to a better understanding of the molecular regulation mechanism in different osmotic stresses.


Background
Saline-alkali water accounts for a considerable proportion of water resources in the world. Salinity has been long recognized as one of the fundamental factors affecting aquatic species distribution and influencing physiological processes of marine and estuarine organisms, such as survival, hemolymph osmolarity, and tissue water content [1][2][3]. Furthermore, salinity has significant effects on physiology of aquatic organisms, and salinity adaptation is a complicated process that involves a series of physiological responses to the environment with different osmotic regulation requirements [4,5]. Carbonate alkalinity stress was considered as a major risk factor for fishes surviving in saline-alkaline water [6,7]. Although many studies have examined osmoregulatory physiology in freshwater and marine populations of fish [8][9][10][11][12], it still remains unknown about the molecular mechanism of osmotic stress tolerance [13]. Therefore, the identification and characterization of the genes and the regulatory factors involved in hyperosmoregulation are now essential for increasing the production and the efficiency of selective breeding programs for some important fish species.
Tilapia is an important commercial fish in China, which is euryhaline fish species and it is the second most important fish after carps in aquaculture [14]. Owing to ease of aquaculture, marketability and stable market prices, the production of tilapia has quadrupled over the past decade [15], and which is promoted and farmed in more than 100 countries and regions [16]. Owing to its importance in tropical and subtropical aquaculture and its extreme euryhaline ability, tilapia has received a high level of scientific interest and is one of the most popular model species for research on fish osmoregulation [17][18][19]. With the increasing scarcity of freshwater available for aquaculture, tilapia that can live in brackish or hyperosmotic water would enable expanding the range of culture and increasing the production of global tilapia [20]. In China one new hybrid tilapia strain, named Mo-Ho tilapia "Guangfu No.1", is a hybrid from Mozambique tilapia Oreochromis mossambicus female × Wami tilapia O. urolepis hornorum male and is tolerant of hyper-osmotic water and grows well in saline pools [21]. It has been widely farmed on the southern coast of China.
With regard to ion-regulation, a set of tissues and organs such as the gills, kidney and gut plays a vital role in the teleosts [22][23][24]. Gill contains complex transport epithelia functions as aquatic gas exchange, acid-base regulation and excretion of nitrogenous wastes [25]. Several studies have shown that the epithelium of gill is reorganized, the distribution, number and size of chloride cells both have changed, considerably to meet the requirements of ion transport and permeability [26,27]. When subjected to various abiotic stresses from ambient water, such as salt and alkali, gills are the first site for sensing stress signals and initiating signaling cascades at the molecular level in response to adverse environments. Through these responses, fish can regulate gene expression of regulatory, functional proteins, and morphology of gill epithelium [28][29][30] to enhance stress tolerance. It is important for current fish biology researches to study fish response to osmotic stress.
Recent years, with the prosperous development of the next-generation sequencing, high-throughput sequencing has been widely used to identify conserved and novel functional genes and signaling pathways in aquatic animals, such as Litopenaeus vannamei [13,31], Eriocheir sinensis [32], Palaemon caridean shrimps [33] and Mytilus trossulus [34]. Guo et al identified three possible osmoregulation-related signaling pathways as lipid metabolism related pathways, tight junction pathway and thyroid hormone signaling pathway in Siberian sturgeon Acipenser baeri in response to salinity stress [35]. Lv et al showed that 552 differentially expressed genes including amino acid transport proteins and ion transport enzymes were obtained in Portunus trituberculatus under low salinity stress [36].
In this study, the transcriptome of the hybrid tilapia (O. mossambicus♀ × O. urolepis hornorum♂) gills was sequenced by using Illumina/Hiseq-2000 RNA-seq technology for the first time. The transcriptomic data between control group and salinity challenge group, alkalinity challenge group and salinity-alkalinity challenge group were compared and analyzed to identify osmoregulation-related genes and pathways. The discoveries of this study can be conductive to illustrate the mechanism of osmotic regulation for aquatic animals adapting to salinity challenges, alkalinity challenge and salinity-alkalinity challenge.

Sequencing and the analysis of reads
To identify mRNAs in hybrid tilapia in response to osmotic stress challenges, mRNA libraries derived from 12 groups were constructed and sequenced using the Illumina deep-sequencing platform. Using high-throughput sequencing, 54,418,529, 61,482,829, 48,353,150 and 44, 867,651 average raw reads were obtained from gill tissue of control group(C), saline challenge group (S), alkali challenge group (A) and saline-alkali challenge group (SA) in hybrid tilapia, respectively (Table 1)

Differential gene expression (DGE) analysis
Based on the expression patterns of DEGs in this study, the values of log 10 FPKM (FPKM: fragments per kilo bases per million reads) were used to cluster the 12 groups, respectively, to measure and identify the similarities among the groups (Fig. 1). Based on the criteria that |logFC| > 1, P-value < 0.05 and FDR ≤ 0.05, the edegR and DESeq2 method were used to conducted differential expression analysis, the difference between these two methods were not significant. According to Sahraeian et al study [37], DESeq2 provided the most accurate differential analysis, so the analysis results of DESeq2 method were used to conduct the following analysis ( Fig. 2). 3772 significant DEGs were detected in group S compared with control group, of which there were 1958 up-regulated genes and 1814 down-regulated genes in group S (Additional file 1: Fig. S1A). We also identified 3412 significant DEGs in group A compared with control group including 1472 up-regulated genes and 1940 down-regulated genes in group A (Additional file 1: Fig.  S1B). Simultaneously, 3050 significant DEGs were detected in group SA compared with group C, of which there were 1315 up-regulated genes and 1735 down-regulated genes in group SA (Additional file 1: Fig. S1C). In addition, 2504 significant DEGs were detected in group SA compared with group S, of which there were 940 up-regulated genes and 1564 down-regulated genes in group SA (Additional file 1: Fig. S1D). Meanwhile, 609 differentially expressed genes in group SA compared with group A, which include 312 upregulated genes and 297 down-regulated genes in group SA (Additional file 1: Fig. S1E). Further analysis indicated that a total of 1250 unigenes showed significantly different expression levels in salinity stress (S), alkalinity stress (A) and salinity-alkalinity stress (SA) (Fig. 3a). And a total of 46

Gene ontology (GO) analysis of significant DEGs
In this study, GO functional analysis indicated 2300 (C vs. S), 2110 (C vs. A) and 1960 (C vs. SA) significant genes in total were classified into 40, 35 and 37 sub-categories of three major categories: biological processes, cellular components and molecular function, respectively (Fig. 4). Through the comparative analysis of GO enrichment results between different groups, it was found that some significant GO terms, such as protein binding transcription factor activity (GO:0000988), behavior (GO:0007610), signaling (GO:0023052) and molecular function regulator (GO:0098772) were enriched uniquely in salt tolerant. While in salt-alkalinity group, reproductive process (GO:0022414) and cell aggregation (GO:0098743) were uniquely enriched. We also found some significant GO terms are most likely to play an essential role in regulating osmotic stress tolerance in hybrids, such as immune system process (GO:0002376), response to stimulus (GO: 0050896), transporter activity (GO:0005215) and electron carrier activity (GO:0009055). And these results of GO enrichments will provide the research clues for our following research.

Kyoto encyclopedia of genes and genomes analysis
The top 30 KEGG pathways with the most number of annotated sequences were listed (Additional file 1: Table  S1, S2, S3). These pathways consisted of cytokine-cytokine receptor interaction, protein processing in endoplasmic reticulum, autoimmune thyroid disease, pyrimidine metabolism, amino sugar and nucleotide sugar metabolism, p53 signaling pathway, biosynthesis of amino acids, glycolysis / gluconeogenesis, DNA replication, proximal tubule bicarbonate reclamation, proteasome, plycine, serine and threonine metabolism, alanine, aspartate and glutamate metabolism, TNF signaling pathway, IL-17 signaling pathway, and so on.

Gene set enrichment analysis (GSEA) analysis of genes
To date, GSEA is being used to identify specific biological processes involved in disease outcomes and bioinformatics analysis in the medical literature and the biological  research. GSEA enrichment tool was used to analysis genes between different groups, and then GO analysis was performed to find which have biological regulate function gene set. After GSEA analysis in this study, multiple GO terms were enriched in salinity group, alkalinity group and salt-alkalinity group (Additional file 1: Table S4, S5, S6) and a vast number genes were annotated to take participated in osmotic stress tolerance.

Identification of significant DEGs and KEGG pathways related to osmoregulation
In addition, in all significant DEGs, some functional genes that related to osmoregulation, including ATP1A (sodium / potassium -transporting ATPase subunit alpha), PRLR (prolactin receptor), NKCC1 (sodium / potassium / chloride transporter), SLC9A3, NHE3 (sodium / hydrogen exchanger), TRPV4 (transient receptor potential cation channel subfamily V member 4), SLC4A2, AE2 (anion exchanger) and CLCN2 (chloride channel 2), were also identified separately. Some genes that may play important roles in osmoregulation in salinity stress, alkalinity stress and salinity-alkalinity stress were listed in Table 2.
A comparative analysis of the KEGG pathways enriched in the three groups revealed that 25 of the pathways were significantly enriched in all the three groups, including steroid biosynthesis (ko00100), IL-17 signaling pathway (ko04657), glycine, serine and threonine metabolism (ko00260), Cell cycle (ko04110), and so on. However, only two pathways, proteasome (ko03050) and the p53 signaling pathway (ko04115), were enriched in the salinity group and the alkalinity group. Meanwhile, three pathways, fructose and mannose metabolism (ko00051), glycolysis / gluconeogenesis (ko00010) and alanine, aspartate and glutamate metabolism (ko00250), were enriched in the salinity and saline-alkali group. While in alkali stress group and salinity-alkalinity group, there were 18 pathways were significantly enriched, include natural killer cell mediated cytotoxicity (ko04650), synthesis and degradation of ketone bodies (ko00072), glycosaminoglycan degradation (ko00531), and so on ( Table 3).

Significant DEGs analysis in different groups
Hierarchical cluster, GO function analysis and KEGG pathway analysis were carried out for significant DEGs in each interval in Venn chart (Fig. 3a). Hierarchical clusters were analyzed of significant DEGs that were expression uniquely in salinity stress (S), alkalinity stress (A), salinity-alkalinity stress (SA), salinity stress and alkalinity stress (S & A), salinity stress and salinity-alkalinity stress (S & SA), alkalinity stress and salinity-alkalinity stress (A & SA), the heat map of significant DEGs were showed in Additional file 1: Fig. S2.
Furthermore, GO analyses were performed to identify similarities and differences among differentially expressed genes only in salinity stress (S), alkalinity stress (A), salinityalkalinity stress (SA), alkalinity and salinity-alkalinity stress (A, SA), salinity and salinity-alkalinity stress (S, SA) and three osmotic stress shared (S, A, SA). By comparing the analysis results of GO enrichments, we found that transcription factor activity, protein binding (GO:0000988) and behavior (GO:0007610) were uniquely enriched in salinity stress (A). We also found that cell aggregation (GO: 0098743), multi-organism process (GO:0051704) and membrane-enclosed lumen (GO:0031974) were uniquely enriched in salinity-alkalinity stress (SA), salinity-alkalinity stress (S, SA) and three osmotic stress (S, A, SA), respectively. All these analysis results will provide research data and clues for the future research directions, and are essential for the research of the molecular mechanism of osmotic regulation in teleost. The results of GO enrichment were listed as Fig. 5.
Meanwhile, KEGG analyses were also performed to identify similarities and differences among differentially expressed genes uniquely in salinity stress (S), alkalinity stress (A), salinity-alkalinity stress (SA), alkalinity and salinity-alkalinity stress ( . We found that cytokine-cytokine receptor interaction, ECM-receptor interaction, glycolysis / gluconeogenesis, ether lipid metabolism, linoleic acid metabolism, carbon fixation in photosynthetic organisms, beta-alanine metabolism, alpha-linolenic acid metabolism, quorum sensing and styrene degradation pathway were uniquely annotated in salinity stress. However, the pathways enriched only in alkalinity stress are mostly related to immunity and disease. Under the salinity-alkalinity stress, mostly pathways were enriched related to metabolism and immunity. The results suggested that compared with salinity stress, alkalinity stress makes tilapia more susceptible to pathogen infection, which more likely leads to diseases.

The validation of differently expressional genes by qRT-PCR
To validate the veracity and reliability of differentially expressed genes identified by mRNA-Seq, we randomly selected 17 genes for qRT-PCR validation from those with different expression patterns based on functional enrichment and pathway results. Melting-curve analysis revealed a single product for all tested genes. Log 2 FCs from qRT-PCR were compared with the mRNA-Seq expression analysis results (Fig. 6). The results of qRT-PCR were significantly correlated with the mRNA-Seq

Discussion
Osmoregulation in tilapia is a complex process because of the diverse range of salinities and alkalinities that they are exposed to in their natural habitat. To be euryhaline, tilapia must be able to cope with salt gain and osmotic water loss when it is in hyper-osmotic water stress. Edwards et al discussed the principles of ion and water transport in euryhaline fishes, and discussed the role of including gills, kidney, urinary bladder, gastrointestinal tract, rectal gland of elasmobranchs, skin, opercular membrane, and yolk sac organs in regulated the osmotic stress [38]. N'Golo et al reported the changes of ionocyte morphology and function after transfer from fresh to hypersaline waters in the tilapia Sarotherodon melanotheron [39]. High-throughput RNA-seq is a good method in illuminating the underlying molecular mechanisms of osmoregulation and immunization [40,41], and the transcriptomic analysis with osmotic stress has been realized in some aquatic vertebrate species and invertebrate species, such as Oreochromis mossambicus [42], Acipenser baeri [35], Oryzias latipes [43], Litopenaeus vannamei [13,31], Eriocheir sinensis [32], Crassostrea gigas [44] and Mytilus trossulus [34]. However, most of these studies were directed to only one osmotic stress treatment, but the natural water environment is complex and changing, so it is important to perform various osmotic stress treatments on fish. In this study, we conducted high-throughput sequencing of three osmotic stresses, including saline stress, alkaline stress and mixed of saline-alkaline stress on hybrid tilapia, in order to understand the osmoregulation mechanism of tilapia under different osmotic stress.
In the present study, 3772, 3412 and 3050 significantly different expression genes were identified in salinity stress, alkalinity stress and saline-alkaline stress, respectively. GO and KEGG analysis were conducted to identified genes and pathways that related to osmoregulation. GO analysis results indicated that cellular process, cell part and binding were the dominant group in the category of biological processes, cellular components and molecular functioning, respectively. These results were consistent with the researches in aquatic animals such as Pacific White Shrimp Litopenaeus vannamei [31], Siberian sturgeon Acipenser baeri [35], Eriocheir sinensis [32] and Ruditapes philippinarum [39]. Numerous ion transport enzymes and ion transporters associated with osmoregulation were identified in the three osmotic stress groups. Osmoregulation related genes, such as NKA, PRL, PRLR, TRPV4, AQP1, NHE3, and solute carrier family have been reported in many studies. Zhu et al found that the NKAα1 mRNA expression levels in the gills of tilapia reached peak level at 24 h after transfer from freshwater to seawater [45]. Yoko et al research showed that PRLR2 facilitates acclimation responses to increased extracellular osmolality [46]. Yuan et al observed the synchronous expression trend of the renal PRLR with pituitary PRL (5d) and the asynchronous expression peaks between branchial (8d) and renal PRLR (5d) in olive flounder Paralichthys   [48]. Muhammad et al suggested that cytoplasmic carbonic anhydrase (ChqCAc) gene is involved in response to changes in pH and in systemic acid-base balance in freshwater crayfish [49]. Furthermore, gene NHE3 was detected and analyzed in Sciaenops ocellatus [50] and Macrobrachium species [51]. Some other genes slc12a1, NKCC2 and AQP1 were also detected and analyzed in Oryzias latipes after transfer from fresh water to seawater [52]. Tatsuya et al also found Na + / K + − ATPase activity increased in specific gills of Octopus ocellatus after transfer from 30-ppt normal seawater to 20-ppt salinity [53]. By matching to KEGG pathway database, possible functions of the significant DEGs were analyzed in order to learn more about the physiological responses of hybrid tilapia to simultaneous salinity stress, alkalinity stress and salinity-alkalinity stress. 39, 51 and 58 KEGG pathways were found in salinity stress, alkalinity stress and salinity-alkalinity stress, respectively. 25 pathways were enriched in the three osmotic stress group. Among these pathways, seven pathways were related to metabolism, eight pathways were related to disease and immunity, and the remaining was related to biosynthesis, cell cycle and proximal tubule bicarbonate reclamation. In deal with hyper-osmotic stress and the energetically demanding salt extrusion, signaling related to DNA replication, mis-match repair, protein biosynthesis and degradation, metabolism and cell cycle regulation were increased [42]. Jiang et al research indicated that amino acid metabolism played crucial roles in osmoregulation [54].
In the proximal tubule bicarbonate reclamation pathway of salinity stress, the gene expression of carbonic anhydrase 4 (CA4), aquaporin-1 (AQP1), sodium-coupled neutral amino acid transporter (SNAT3, solute carrier family 38 member 3), malate dehydrogenase (MDH1) and phosphoenolpyruvate carboxykinase (GTP) (PEPCK) were significantly upregulated. However, the expression of sodium bicarbonate cotransporter (NBC1, solute carrier family 4 member 4) was significantly decreased. While in the alkalinity stress, the expression of NBC1 was significantly up-regulated, but the expression of NHE3 was similar with the control group. However, in salinity-alkalinity stress, the expression of NHE3 was significantly downregulated, NBC1 gene was up-regulated, but the expression of AQP1 was similar to the control group. Those genes plays pivotal role in maintaining the balance of proximal tubular cell and the regulation of blood pH. The  renal proximal tubule is responsible for most of the renal sodium, chloride, and bicarbonate reabsorption [55].
Alexander et al showed that NHE3 participates in the reabsorption of Na + , bicarbonate, and water from the proximal tubule, contributing to the maintenance of intravascular volume, blood pressure, and acid-base homeostasis [56]. Masashi et al found that a majority of renal proximal bicarbonate absorption is accomplished by coordinated operation of the apical NHE3 and the basolateral electrogenic Na + / HCO 3 − cotransporter NBC1 [57]. Defects in NHE3 or NBC1 can theoretically reduce the reabsorption of filtered HCO 3 − and cause proximal renal tubular acidosis (pRTA) [58][59][60][61]. Gordon et al reported that AQP1 enhances the reabsorption of HCO 3 − by the renal proximal tubule by increasing the CO 2 permeability of the apical membrane and the macroscopic cotransport of Na + plus two HCO 3 − occurs as NBC transports Na + plus CO 3 2− and AQP1 transports CO 2 and H 2 O [62]. Our results indicated that proximal tubule bicarbonate reclamation pathway was evidently crucial for maintaining systemic acid-base homeostasis in osmotic stress of tilapia.  Nine pathways were uniquely enriched in salinity stress, such as ECM-receptor interaction, protein processing in endoplasmic reticulum, nitrogen metabolism and arachidonic acid metabolism. Among these pathways, arachidonic acid metabolism playing important role in osmoregulation has been reported in Eriocheir sinensis [32], Litopenaeus vannamei [31], Ruditapes philippinarum [39] and Litopenaeus vannamei [63]. Six KEGG pathways were characteristic enriched in alkalinity stress group. Meanwhile, 12 KEGG pathways, including pantothenate and CoA biosynthesis, mineral absorption, valine, leucine and isoleucine biosynthesis, cell adhesion molecules (CAMs), glyoxylate, dicarboxylate metabolism and propanoate metabolism, etc., were uniquely annotated in salinity-alkalinity stress. In the mineral absorption pathway, TRPV6, chloride channel 2 (clcn2), zinc transporter (ZnT1), Na + / K + -exchanging ATPase and femitin (FTH1) were upregulated significantly. Meanwhile, NHE3 and heme oxygenase 1 (HMOX1) were downregulated significantly. Through the network formed of all genes, this pathway played vital role in ion regulation and maintain the balance of osmotic stress in internal and external of the tilapia.
From KEGG pathway analysis of significant DEGs uniquely in salinity stress group, three enriched pathways related to lipid metabolism including alpha-linolenic acid metabolism, linoleic acid metabolism and ether lipid metabolism were examined. The long-chain polyunsaturated fatty acid (LC-PUFA) level raised in euryhaline and migratory teleosts during salinity acclimation has been reported by many researches, such as Oncorhynchus mykiss [64], Dicentrarchus labrax [65], and Oncorhynchus masou [66]. LC-PUFAs, as structural components of plasma membrane,  enhanced the fluidity of cell membranes [67], results in the changes of other physical membrane properties and affects the membrane-protein binding [68]. These results could be used to explain the increase of LC-PUFAs level in euryhaline and migratory teleosts during salinity acclimation. While in alkalinity stress, the enriched pathways were mostly related to disease, pathogen defense and innate immune response, suggesting that tilapia was more vulnerable to alkaline stress than saline stress. While in salinity-alkalinity stress group, most pathways were related to biosynthesis and metabolism, such as pantothenate and CoA biosynthesis, neomycin, kanamycin and gentamicin biosynthesis, nicotinate and nicotinamide metabolism and 2-Oxocarboxylic acid metabolism. In saline & alkaline, the enriched pathways of significant DEGs mostly related to antigen processing and antibody production. In saline & saline-alkaline, most pathways were related to secretion and metabolism, such as gastric acid secretion, pancreatic secretion, nitrogen metabolism and propanoate metabolism. While in alkaline & salinealkaline, most pathways were related to immunity, and mismatch repair, nucleotide excision repair, DNA replication, glycine, serine and threonine metabolism, cell cycle, sphingolipid metabolism, glycosaminoglycan degradation, and so on, which were related to osmoregulation. Inokuchi et al showed that the apoptosis of freshwater-type mitochondrial rich cells would be enhanced and stimulated the recruitment of mitochondrial rich-cell complexes to develop hypo-osmoregulatory ability for seawater adaptation after transfer from hypo-osmotic to hyper-osmotic water in Oreochromis mossambicus [69]. Moreover, osmotic stress has been reported to cause oxidative stress, protein damage, DNA double-strand breaks and oxidative base modification [70][71][72]. These results are in consistent with the enrichment of canonical pathways such as mismatch repair, nucleotide excision repair, DNA replication, glycine, serine and threonine metabolism and cell cycle in response to osmotic tolerance of the hybrid tilapia. In saline & alkaline & salinity-alkalinity, 28 pathways were found, of which 10 pathways were related to metabolism, such as glutathione metabolism, glycine, serine and threonine metabolism, purine metabolism and pyrimidine metabolism.

Conclusions
Here, the transcriptome of the hybrid tilapia response to saline stress, alkaline stress and saline-alkaline stress was examined following the application of RNA-seq. Through more than 40 million sequence reads and the expression data of the four libraries, the study provided some useful insights into signal transduction pathways in hybrid tilapia and offered a number of candidate genes as potential markers of tolerance to osmotic stress for tilapia. The present research described for the first time the effects of salinity, alkalinity and salinity-alkalinity on gilled classic osmoregulation-related genes and signaling pathways during acclimation to saline and alkaline environment in tilapia. Our results provide a deep insight into the gill osmotic stress response of the hybrid tilapia to saline stress, alkaline stress and the mix of saline-alkaline stress that will be helpful in identifying reliable osmoregulation target genes and pathways for fish in osmotic stress.

Sample collection
The gill tissue samples were obtained from the hybrid tilapia (O. mossambicus♀ × O. urolepis hornorum♂) specimens at the Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Guangzhou, China. All fishes were acclimated in the laboratory environment (28°C) for two weeks before the experiment. The average weight of the experimental fishes was 51.27 ± 3.78 g. Three treatments with three replicates each were set up. The Following anesthesia, the fishes were decapitated and the gills were dissected immediately and immersed into liquid nitrogen for high-throughput sequencing. After 24 h of osmotic stress, the remaining fishes were released.

RNA isolation, library construction and illumina sequencing
The TRIzol® reagent (Invitrogen, USA) was used to extract the total RNA from four samples according to the manufacturer's protocol and then treated by DNase I to implement the DNA digestion and obtained pure RNA products. Finally, the integrity of the total RNA was confirmed using 1% agarose gel electrophoresis and an Agilent 2100 Bio analyzer (Agilent Technologies, USA). Equal amounts of RNA from different individuals in the same group were pooled for library construction. Next generation sequencing library preparations were constructed according to the manufacturer's protocol (NEBNext® Ultra™ RNA Library Prep Kit for Illumina®). Subsequently, the deep sequencing of the four mRNA libraries was performed using Illumina HiSeq 2000 according to the manufacturer's instructions at the GENEWIZ (Genewiz, Suzhou, China).

Basic analysis of sequencing data
In order to remove technical sequences, including adapters, polymerase chain reaction (PCR) primers, or fragments thereof, and the quality of bases lower than 20, the pass filter data in the fastq format were processed by Trimmomatic (v0.30) to provide high-quality, clean data. Firstly, the reference genome sequences and the gene model annotation files of relative species were downloaded from genome website, such as NCBI. Secondly, Hisat2 (v2.0.1) was used to indexed with reference genome sequence. Finally, we mapped the clean reads to reference sequences of Oreochromis niloticus (National Center for Biotechnology Information, NCBI; accession numbers: GCA_001858045.3, https://www.ncbi.nlm.nih.gov/assembly/GCA_001858045.3 ) via software Hisat2 (v2.0.1). In the beginning transcripts in fasta format are converted from known gff annotation file and indexed properly. Then, with the file as a reference gene file, HTSeq (v0.6.1) estimated gene and isoform expression levels from the pair-end clean data.

Differential expression analysis and bioinformatics analysis of genes
DESeq2 and edegR method were used to conducted differential expression analysis. FPKM (Fragments Per Kilobases per Millionreads) method can eliminate the influence of gene length and sequencing amount on the calculated gene expression, and the calculated gene expression level can be directly used to compare the expression differences between different genes, so it was used to normalized the data. After adjusting by BH (fdr correction with Benjamini/Hochberg) for controlling the false discovery rate, the P-values < 0.05 and |log2FC| > = 1 were set to detect significant differentially expressed ones [73,74]. The gene ontology (GO) enrichment analysis indicated that the significant differentially expressed genes with a significant P-value less than 0.05 were enriched for Fig. 6 Validation of RNA-seq data by qPCR. To validate the RNA-seq data, the relative mRNA levels of 17 randomly selected DEGs in the gills of hybrid tilapia were examined by qPCR. The mRNA levels by qPCR are presented as the fold change compared with the mock-treated control after normalization against β -actin. The relative expression levels from the RNA-seq analysis were calculated as log 2 FC values. The gene nkai1 was calculated in two different groups. (Note: nkain1-1 was in a group C vs. A, nkain1-2 was in group C vs. SA; trpv4-1 was in group S vs. SA, trpv4-2 was in group C vs. S) biological processes, molecular functioning, and cellular components through the software GO-TermFinder. KEGG (Kyoto Encyclopedia of Genes and Genomes) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances (https://www.genome.jp/ kegg/kegg3.html) [75]. We used scripts in house to enrich significant differential expression gene in KEGG pathways. The Venny online software (http:// bioinfogp.cnb.csic.es/tools/venny/) was used to combine the analyses of differentially expressed genes. Gene Set Enrichment Analysis (GSEA) was also used to detect the expression changes of the whole gene set. This enrichment tool can fully detect genes that are not significantly different in expression levels but are plays important biological role at the overall level.

RNA-Seq data validation by qPCR and statistical analysis
To validate and measure the differential expression of mRNAs identified by high-throughput sequencing, 16 differentially expressed mRNAs were randomly selected for real-time quantitative PCR analyses. Real-time quantitative PCR was performed using the SYBR Green PCR Master Mix (Life Technologies, USA) according to the manufacture's protocol in an ABI 7300 Real-Time PCR System (Applied Biosystems, Foster City, CA) using the following conditions: 50°C for 2 min, 95°C for 2 min, 40 cycles of 95°C for 15 s, 60°C for 15 s and 72°C 30s. At the end of qPCR assay, the melting curve was analyzed to evaluate the amplification specificity. All primer sequences are listed in Table 11. The relative expression was calculated using the 2 -△△Ct method and β-actin was used as the internal control to normalize mRNAs expression. The results are presented as the mean ± standard error. The analysis of all data was implemented using the statistical software SPSS 22.0 (SPSS, Chicago, IL, USA) and the expression levels of mRNAs were analyzed using one-way ANOVA method.
Additional file 1. Summary results of significant DEGs analysis and the results of GO, KEGG enrichment.

Acknowledgements
We would like to thank the reviewers for their helpful comments on the original manuscript.
Author's contributions HS and HZ designed the experiment and drafted the manuscript, HS performed the experiments and analyzed the data. DM, ZL and FG participated in experimental design and coordination and contributed to the manuscript preparation. All authors reviewed and approved the manuscript.  Availability of data and materials All the data supporting our findings are contained within the manuscript.
The clean sequence reads dataset supporting the conclusions of this article are available on NCBI BioProject with accession number PRJNA558948 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA558948) and will be made publicly available upon acceptance of this manuscript for publishing.

Ethics approval and consent to participate
The hybrid tilapia used in this study was produced by at the Pearl River Fisheries Research Institute, Chinese Academy of Fishery Sciences, Guangzhou, China. All experiments were performed according to the Guide for the Care and Use of Laboratory Animals of China. This study was also approved by the Animal Care and Use committee of the Centre for Applied Aquatic Genomics at the Chinese Academy of Fishery Sciences.

Consent for publication
Not applicable.