Comparative transcriptome analysis of Triplophysa yarkandensis in response to salinity and alkalinity stress

Journal Pre-proof Comparative transcriptome analysis of Triplophysa yarkandensis in response to salinity and alkalinity stress

Sheng-Ao Chen, Jilun Hou, Na Yao, Congxin Xie, Dapeng Li PII:

S1744-117X(19)30158-3

DOI:

https://doi.org/10.1016/j.cbd.2019.100629

Reference:

CBD 100629

To appear in:

Comparative Biochemistry and Physiology - Part D: Genomics and Proteomics

Received date:

25 March 2019

Revised date:

11 September 2019

Accepted date:

12 September 2019

Please cite this article as: S.-A. Chen, J. Hou, N. Yao, et al., Comparative transcriptome analysis of Triplophysa yarkandensis in response to salinity and alkalinity stress, Comparative Biochemistry and Physiology - Part D: Genomics and Proteomics(2019), https://doi.org/10.1016/j.cbd.2019.100629

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier.

Journal Pre-proof Title: Comparative transcriptome analysis of Triplophysa yarkandensis in response to salinity and alkalinity stress Sheng-Ao Chen a, b, †, Jilun Houc, d, †, Na Yaob, Congxin Xie a, Dapeng Lia, * a

College of Fisheries, Huazhong Agricultural University, Wuhan 430070, China

b

College of Animal Science, Tarim University, Alar, 843300, China

c

Key Laboratory of Aquatic Genomics, Ministry of Agriculture, Beijing, 100141, China

d

Beidaihe Central Experimental Station, Chinese Academy of Fishery Sciences, Qinhuangdao,

066100, China *

ur

na

lP

re

-p

These authors contributed equally to this work as first authors.

Jo

†

ro of

Corresponding author. E-mail: [email protected](D. Li)

Journal Pre-proof Abstract: Triplophysa yarkandensis, a fish belonging to the family Nemacheilidae, is distributed in the Tarim River, China, immediately north of the Qinghai-Tibet Plateau. Due to increasing salinity and alkalinity in the Tarim River, the habitats of T. yarkandensis have been seriously altered. To identify the genes and pathways that are important for responding to salinity and alkalinity stress, the gill transcriptomes of fish living under different salinity and alkalinity conditions were obtained using RNA sequencing. A total of 1123448964 clean reads were obtained and assembled into 177271 unigenes, with an average length of 1703 bp. Around 13526 unigenes showed differential expression when comparing different salinity concentrations with the

ro of

controls, 6967 of which were upregulated and 6559 were downregulated. When comparing different alkalinity concentrations with the controls, there were 17475 unigenes that showed

-p

differential expression, of which 10457 were upregulated and 7018 were downregulated. Only 146 unigenes were both differentially expressed in salinity and alkalinity groups compared to the

re

control. The results of KEGG enrichment showed that there were five upregulated and 12

lP

downregulated pathways in fish subject to salinity treatment. For fish exposed to alkalinity treatment, 15 pathways were upregulated and 13 downregulated. There were four upregulated and

na

four downregulated pathways that were shared by fish subject to salinity and alkalinity treatments. To our knowledge, this is the first study on the T. yarkandensis transcriptome; the information

ur

presented here will provide further understanding of the fish’s response to salinity and alkalinity stress, as well as further insight into the T. yarkandensis genome.

Jo

Keywords: salinity, alkalinity, RNA-Seq, Tarim River

1. Introduction

The Tarim River is the longest inland river in China, with a total length of 1321 km. It is a typical seasonal river in an arid area and does not produce runoff. The biodiversity in the Tarim River ecosystem is relatively concentrated and fragile. The main fish species caught in this river are Aspiorhynchus laticeps, Schizothorax biddulphi, and Triplophysa yarkandensis. Triplophysa yarkandensis, locally named “dog-headed fish,” belongs to the order Cypriniformes in the Nemacheilidae family, within the Triplophysa genus and Hedinichthys subgenus(Nie et al., 2013). It is an endemic and dominant species, widely distributed in the Tarim River. T. yarkandensis is

Journal Pre-proof also one of the main components of the Qinghai-Tibet Plateau fish population and one of the highest-altitude fish species(Chen et al., 2010). Due to the high nutritional and economic value of T. yarkandensis, this species is subjected to overfishing, resulting in the reduction of this species. In addition, serious damage to the Tarim River ecosystem and a decrease in the amount of water in Tarim River are causing an increase in salinity and alkalinity, which is seriously impacting the habitats of indigenous fish. In 2010, the ionic concentrations in the Tarim River basin were eight-fold higher than in 1965; the Cl- concentration in particular, was 21.7-fold higher(Wang et al., 2013). In addition, the dissolved oxygen content, salinity, pH, and conductivity were found to be

ro of

high in the Alar region of the Tarim River(Chen et al., 2014).

Salinity and alkalinity are two important environmental factors that can have direct effects on

-p

the physiological status of aquatic animals. Salinity changes can cause oxidative damage due to an increase in reactive oxygen species (ROS)(Choi et al., 2008; Li et al., 2014; Lushchak, 2011;

re

McNamara et al., 2004). The effect of alkalinity stress on fish involves synergistic toxicity of total

lP

alkalinity, pH, and the resulting changes in metabolism of ammonia nitrogen. which HCO3concentration in fish blood increases, which in turn increases the pH, disrupting the acid-base balance in vivo. The increase in blood pH leads to a decrease in H+ concentration, such that the

na

metabolic NH3 is unable to combine with H+ to be excreted as NH4+, resulting in “ammonia

2012).

ur

poisoning”(Gilmour et al., 2007; Wilkie et al., 1994; Wilkie and Wood, 1996, 1994; Yao et al.,

Jo

Recently, RNA sequencing (RNA-seq) has been widely used to study gene expression profiles and to conduct pathway analysis on aquatic animals’ subject to salinity and alkalinity stress. RNA-Seq has been proven to be an adequate tool for capturing the gene and pathway responses(Chen et al., 2015; Li et al., 2014; Wang et al., 2015; Xu et al., 2017, 2013a, 2013b; Yao et al., 2012; Zhang et al., 2019; Zhao et al., 2015). To our knowledge, no transcriptome study has been performed on T. yarkandensis to date. In the present study, we performed RNA-seq to investigate the gene expression profile of T. yarkandensis that was challenged with different concentrations of salinity and alkalinity. Our results will further enhance the understanding of T. yarkandensis response mechanisms to salinity and alkalinity stress, and provide basic data for the protection and culture of this species.

Journal Pre-proof 2. Materials and methods 2.1 Ethics Statement This study was approved by the Animal Care and Use Committee of Tarim University. 2.2 Fish and fertilization In this study, fish from the second filial (F2) generation of wild captured T. yarkandensis were used as male and female parents. The fish were reared in the aquarium room of Tarim University. Spawning was induced using a single injection of hCG, with a dosage of 20 IU/g body weight.

ro of

After 12–14 h of rearing at 21 ± 0.5 °C, the gametes were collected. Eggs and sperm were mixed, and ambient tap water was added to obtain fertilized eggs, which were incubated for hatching. 2.3 Salinity and alkalinity exposure

-p

Fish with body lengths of 69.13 ± 11.52 mm and body weights of 5.02 ± 1.36 g were used for

re

salinity and alkalinity exposure. Prior to exposure, the fish were reared for 7 days in the experiment environment. During the experiment, the water temperature was 20 ± 0.5 °C.

lP

Salinity concentrations were prepared using NaCl, and alkalinity concentrations were prepared using NaHCO3. According to the salinity and alkalinity of Tarim River (4～5‰ for

na

salinity, 2.7～14.4 mmol/L for alkalinity), also to the safety concentrations of salinity (1.5‰) and alkalinity (6.9 mmol/L) obtained from previous study (unpublished data), the salinity

ur

concentrations were set as 2.0, 4.0, and 6.0‰, and the alkalinity concentration gradients were set

Jo

to 6.9, 9.9, and 12.9 mmol/L (pH: 7.2~7.4). The exposure test was performed in 30 L tanks. Each concentration was repeated in triple with 20 individuals in each tank. During exposure, each tank was supplied with artificial diets three times a day (10% of body weight). One hour after feeding, the residual feed and dirt was removed, fresh water was added, and the salinity and alkalinity concentrations were readjusted. 2.4 Measuring of growth trait After 30 days of exposure, all fish in each concentration were measured for total length. The absolute growth rate (AGR), relative growth rate (RGR), instantaneous growth rate (IGR) was calculated according to Yin (1995). 2.5 Tissue sampling and RNA extraction After 30 days of the experiment, gills from three individuals of each salinity and alkalinity

Journal Pre-proof concentration together with those from the control group were collected and stored in RNAlater (AM7021, Ambion, CA, USA), and stored at -20 °C. In total, 21 samples (nine samples from the salinity conditions, nine from alkalinity, and three from control) were collected. Total RNA extraction was carried out using TRIZOL Kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer instructions. The extracted RNA was treated with RNase-free DNase I to remove potential genomic DNA. Finally, the quality of each RNA was determined by 1% agarose gel electrophoresis and a NanoPhotometer spectrophotometer (IMPLEN, CA, USA). 2.6 Library preparation and sequencing

ro of

A total sample of 1.5 μg RNA was used as the input material for the RNA-seq library preparations. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina

-p

(NEB, USA). Following manufacturer’s recommendations, index codes were added to attribute sequences to each sample. The quality of the prepared library was assessed using the Agilent

re

Bioanalyzer 2100 system.

lP

After preparation of the library, the clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego,

na

CA, USA) according to the manufacturer’s instructions. After cluster generation, the libraries were sequenced using an Illumina Hiseq 2000 platform and paired-end reads were generated.

ur

2.7 De novo assembly and unigene annotation In order to obtain clean reads, the adaptor sequences in the raw reads were trimmed, and reads

Jo

containing poly-N (percentage of N > 10%) together with reads that contained more than 10% of bases with a PHRED quality ≤ 20 were removed. All downstream analyses were based on clean reads with high quality. De novo assembly of the transcriptome was performed using Trinity software (v2.4.0)(Grabherr et al., 2011) with min_kmer_cov set to 3 and all other parameters set to default. Corset software (v1.05) was then used to cluster the transcripts into unigenes using default parameters(Davidson and Oshlack, 2014). Subsequently, all the assembled unigenes were compared for genetic similarity to the NCBI non-redundant protein sequences database (NR) (parameters: e-value = 1e-5, --more-sensitive), the euKaryotic Ortholog Groups (KOG) (parameters: e-value = 1e-3, --more-sensitive), and the Swiss-Prot protein database (parameters: e-value = 1e-5, --more-sensitive) using DIAMOND software (v0.8.22)(Buchfink et al., 2015). Unigenes were also compared to the NCBI nucleotide

Journal Pre-proof sequences database (NT) using the NCBI blast program (v2.2.28+) (parameter: e-value = 1e-5), to the Protein family database (Pfam) using the hmmscan package (HMMER 3) (parameter: e-value = 0.01), to the Kyoto Encyclopedia of Genes and Genomes database (KEGG) using the KEGG Automatic Annotation Server (KAAS) (r140224) (parameter: e-value = 1e-10)(Moriya et al., 2007), and to the Gene Ontology database (GO) using blast2go (b2g4pipe_v2.5) (parameter: e-value = 1e-6)(Götz et al., 2008). 2.8 Differential expression analysis Gene expression levels for each sample were estimated by mapping back clean reads onto the

ro of

assembled transcriptome using Bowtie2(Langmead and Salzberg, 2012). The read count for each gene was then obtained from the mapping results using RSEM(Li and Dewey, 2011). The

-p

identification of differentially expressed genes (DEGs) was performed using DESeq2 (v1.6.3)(Love et al., 2014). Genes with an adjusted P value < 0.05 and |log2FoldChange| > 1 found

re

by DESeq2 were assigned as differentially expressed. For each salinity or alkalinity concentration,

lP

the DEGs were obtained by comparing to control or other concentration. All salinity concentration groups were taken together, and regarded as a total salinity group; and for alkalinity groups, all

na

concentration groups were also taken together as a total alkalinity group. The DEGs were then compared among the total salinity group, total alkalinity group and control.

ur

2.9 Functional annotation

GO enrichment analysis of the DEGs was carried out using the GOseq R packages (v1.10.0)

Jo

(parameter: corrected P value < 0.05), implementing a Wallenius non-central hyper-geometric distribution technique(Young et al., 2010). KOBAS software (v2.0.12) (parameter: corrected P value < 0.05) was used to test the statistical enrichment of DEGs in KEGG pathways. 2.10 qRT-PCR validation To validate the RNA-seq results, nine DEGs with high levels of significance were selected for qRT-PCR analysis, using β-actin as the internal control. The primers are listed in Table S1. The PCR was performed with SYBR Green dye using an ABI PRISM 7500 Fast Real-Time PCR machine under the following thermal cycling conditions: 95 °C for 3 min, followed by 40 cycles of 95 °C for 30 s and 55 °C for 1 min, and 72°C for 20 s. The 2-ΔΔCt method was used to analyze the expression levels of the target genes(Livak and Schmittgen, 2001).

Journal Pre-proof

3. Results 3.1 Growth traits after exposure For salinity exposure, the AGR, RGR and IGR in 4.0‰ and 6.0‰ groups were significantly lower than 2.0‰ group and control (P ＜ 0.05).There were no significant differences between 2.0‰ group and control in these values. For alkalinity exposure, the AGR and RGR in 6.9 mmol/L group were significantly higher than 9.9 and 12.9 mmol/L group (P ＜ 0.05), but not

groups, but lower than control (P ＜ 0.05) (Table 1). 3.2 RNA-Seq and de novo assembly

ro of

different from control. The 9.9 mmol/L group had highest IGR in all alkalinity concentration

-p

For 21 samples, a total of 1182224012 raw reads were yielded. After filtering the ambiguous nucleotides and low-quality sequences, a total of 1123448964 clean reads (95.03% of raw reads)

re

were obtained. For each sample, the number of clean bases ranged from 6.2 to 9.65 G, and the

lP

Q20 and GC content ranged from 95.61 to 98.28% and 41.53 to 44.20%, respectively (Table 2). As shown in Table 3, the clean reads were assembled into 177271 unigenes, ranging from 201 to

na

19318 bp in length. The average length was 1703 bp, the median length was 1043 bp, and the N50 and N90 were 2923 bp and 719 bp, respectively. Most unigenes (29.76%) were > 2000 bp,

ur

followed by the 24.89% of unigenes clustered within the 501–1000 bp group. Fewest unigenes (only 5.66%) were in the < 300 bp group (Fig. 1). The transcriptome sequencing data generated in

Jo

this study have been deposited in the NCBI short read archive (SRA) associated with BioProject: PRJNA527357.

3.3 Unigene functional annotation and classification All 177271 of the unigenes from the gill of T. yarkandensis were used as queries for annotation in seven databases. Out of 177271 unigenes, 140923 (79.50%) were annotated in at least one database and 6716 unigenes (3.79%) were annotated in all databases. For databases NR, NT, KO, Swiss-Prot, Pfam, GO, and KOG, the number of annotated genes ranged from 14679 to 111884 (8.28 to 63.11%) (Table 4). A total of 97651 (55.09%) unigenes were assigned to GO terms. These GO terms were grouped into the three major categories, “cellular component,” “biological process,” and “molecular function,” with 56 subcategories (Fig. 2). Within the “cellular component” category,

Journal Pre-proof the “cell” and “cell part” subcategories shared the same frequency of unigenes (34119, 18.53%) and were the most dominant groups, followed by “organelle” (23240, 12.62%) and “macromolecular complex” (23008, 12.49%) groups. Within the “biological process” category, most unigenes were found in the “cellular process” subcategory (59852, 19.76%), followed by the “metabolic process” (51056, 16.85%), “single-organism process” (47539, 15.69%), and “biological regulation” (27455, 9.06%) subcategories. Within the “molecular function” category, 58143 (45.63%) of unigenes were assigned to “binding,” followed by 42119 (333.06%) to “catalytic activity,” 8773 (6.89%) to “transporter activity,” and 4676 (3.67%) to “molecular

ro of

transducer activity.”

The results of the KOG annotation indicated that a total of 50598 unigenes (28.54%) were classified into 26 categories (Fig. 3). Among these categories, the major category was “signal

-p

transduction mechanisms” (7853, 15.52%), followed by “general function prediction” (7291,

re

14.41%), “posttranslational modification, protein turnover, chaperones” (5565, 10.99%), and

lP

“translation, ribosomal structure, and biogenesis” (4494, 8.88%). The KEGG BLAST search results showed that 14679 unigenes (8.28%) were assigned to 231 pathways (Table S2). These pathways were classified into five hierarchies: “cellular processes,”

na

“environmental information processing,” “genetic information processing,” “metabolism,” and

ur

“organismal systems.” For each hierarchy, four, three, four, twelve, and nine sub-hierarchies were categorized, respectively (Fig. 4). The number of unigenes annotated in KEGG pathways ranged

Jo

from 1 to 568. “D-arginine and D-ornithine metabolism,” and “C5-branched dibasic acid metabolism” pathways had only one and three annotated unigenes, respectively, while all the other pathways had no less than five unigenes. 3.4 Analysis of DEGs We first compared the DEGs from fish subject to different salinity and alkalinity concentrations. For salinity treatments, the number of upregulated unigenes ranged from 181 to 7452 and the number of downregulated unigenes ranged from 181 to 10961 (determined by pairwise comparison) (Table 5). For alkalinity treatments, 792 to 15270 unigenes were upregulated and 2422 to 15505 unigenes were downregulated (Table 6). When all salinity concentrations were compared to controls, 13526 unigenes were differentially expressed; 6967 unigenes were upregulated and 6559 downregulated (Fig. 5A). When comparing all alkalinity concentrations

Journal Pre-proof with controls, there were 17475 DEGs, 10457 of which were upregulated and 7018 downregulated (Fig. 5B). Comparing alkalinity treatments with salinity treatments, there were 3746 DEGs, 3216 of which were upregulated and 530 downregulated (Fig. 5C). Further analysis showed that only 146 unigenes were both differentially expressed in salinity and alkalinity groups compared to the control (Fig. 6). The profiles of DEGs among salinity, alkalinity, and control groups were presented using a heatmap analysis of hierarchical clustering (Fig. 7). The DEGs were then compared to the KEGG database for pathway enrichment. By considering all salinity concentrations as one group, there were five upregulated pathways:

ro of

“transcriptional misregulation in cancer,” “protein export,” “MAPK signaling pathway,” “collecting duct acid secretion,” and “fatty acid elongation”; and 12 downregulated pathways (Table 7). For alkalinity treatment, 15 pathways were upregulated and 13 were downregulated

-p

(Table 8). Pathways such as “transcriptional misregulation in cancer,” “collecting duct acid

re

secretion,” “fatty acid elongation,” and “protein export” were upregulated in both salinity and alkalinity treatments. Downregulated pathways like “NOD-like receptor signaling pathway,”

lP

“pathways in cancer,” “Ras signaling pathway,” and “viral carcinogenesis” were shared by both

3.5 qRT-PCR validation

na

treatments.

ur

The reliability of the results obtained from RNA-Seq was validated by qRT-PCR. We randomly selected nine DEGs, including fubp3, glud, galphai2, gbas, slco2b1, pdlim5, prkci, zfyve16, and

Jo

abcc1 for PCR amplification. Only one single product was detected from each primer set by melting curve analysis. The results of qRT-PCR were consistent with RNA-Seq (Fig. 8), indicating that the RNA-Seq and Trinity assembly was accurate.

4. Discussion In the present study, we reported the first RNA-Seq study in a non-model fish species, T. yarkandensis. A total of 177271 unigenes were obtained from the gills of T. yarkandensis using RNA-Seq. The length of the 177271 unigenes ranged from 201 to 19318 bp. The N50 was 2923 bp, which is higher in length than that of other fish species (Tian et al., 2019; Xu et al., 2013a; Zhang et al., 2019). Approximately 79.50% of the unigenes were annotated in at least one database, which was higher than the reported results for Micropterus salmoides (52.98%)(Zhang

Journal Pre-proof et al., 2019), Cyprinus carpio (77.3%)(Ji et al., 2012), Sparus aurata (51%)(Calduch-Giner et al., 2013), Leuciscus waleckii (57.6%)(Xu et al., 2013a), and Gymnocypris przewalskii (38.62%)(Tian et al., 2019). To our knowledge, no study has been conducted on the T. yarkandensis genome to date. The findings of the present study provide a large amount of novel transcription data for further genomic study of this species. Under salinity stress, pathways such as “transcriptional misregulation in cancer,” “protein export,” “MAPK signaling pathway,” “collecting duct acid secretion,” and “fatty acid elongation” were upregulated. The “MAPK signaling pathway” is a highly conserved module in all eukaryotes

ro of

and is involved in various cellular functions, including cell proliferation, differentiation, and migration(de Zelicourt et al., 2016; Tian et al., 2019). In fish, the MAPK genes play central roles

-p

in response to salinity stress. These include activation of downstream targets that mediate physiological acclimation via protein phosphorylation cascades, and regulation of intracellular

re

concentrations of inorganic ions and organic osmolytes(Fiol and Kültz, 2007; Kültz and Avila,

lP

2001; Tian et al., 2019; Zhou et al., 2016). In the present study, the expression levels of mapk gene such as mapk1/3, map2k4, map3k5 were significantly upregulated, indicating their potential

na

functions in salinity stress responses in T. yarkandensis. There were 15 upregulated pathways in the alkalinity stress groups. Among them, four were

ur

related to “human disease,” four to “organismal systems,” four to “metabolism,” and three to “genetic information processing.” Interestingly, many V-ATPase genes were upregulated in

Jo

pathways such as “oxidative phosphorylation” and “synaptic vesicle cycle.” V-ATPase is a multi-subunit membrane protein complex that plays an important role in the activation of ion and nutrient transport(Wang et al., 2016). It regulates the pH of intracellular compartments, the cytoplasm, and the extracellular space(Hinton et al., 2009). In fish, response to pH change requires efficient mechanisms to regulate water and ion concentrations and movements (Hwang and Lee, 2007). Long-term regulation of the acid-base balance in fish exposed to alkaline water is dependent upon the independent modulation of Na+ and Cl- movements across the gill epithelium(Wilkie and Wood, 1996). Knockdown of V-ATPase subunit A in zebrafish (Danio rerio) impaired acid secretion and ion balance(Horng et al., 2007). Thus, we propose that T. yarkandensis may use similar mechanisms to other fish species for acid-base balance and it is possible that this mechanism has been conserved through evolution.

Journal Pre-proof For salinity and alkalinity stress, the majority of the downregulated pathways were immunity or disease related. Downregulation of the immune response under alkalinity stress was also reported in medaka (Oryzias latipes) (Yao et al., 2012). It is well known that stress can induce immune dysfunction in animals (Glaser and Kiecolt-Glaser, 2005). In fish, chronic stress has been shown to suppress the immune response, therefore enhancing the opportunity for infection(Bly et al., 1997; Tort, 2011). In the present study, the duration of salinity and alkalinity stress was 30 days, and the salinity and alkalinity concentrations were moderate. This brought chronic stress to T. yarkandensis and therefore suppressed immune-related pathways. However, the detailed

ro of

mechanisms for how salinity and alkalinity downregulate immune responses remain unknown and

-p

require further study.

5. Conclusions

re

In the present study, transcriptomes were generated from the gills of T. yarkandensis that were

lP

subjected to different concentrations of salinity and alkalinity. The DEGs and enriched pathways were identified. To our knowledge, this is the first transcriptome study on T. yarkandensis. The

na

transcriptomic data presented here will provide further understanding of the fish’s response to

ur

salinity and alkalinity stress, and further insight into its genome.

Conflict of interest

Jo

The authors have no conflicts of interest to declare.

Acknowledgements

This work was financially supported by National Natural Science Foundation of China (31360635, 31800391), Special Grant of Tarim University (TDZKGG201604), and the Basic Science and Technology Program of Xinjiang Production and Construction Corps (2017DB003).

References Bly, J.E., Quiniou, S.M., Clem, L.W., 1997. Environmental effects on fish immune mechanisms. Dev. Biol. Stand. 90, 33–43. Buchfink, B., Xie, C., Huson, D.H., 2015. Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60. Calduch-Giner, J.A., Bermejo-Nogales, A., Benedito-Palos, L., Estensoro, I., Ballester-Lozano, G., Sitjà-Bobadilla, A., Pérez-Sánchez, J., 2013. Deep sequencing for de novo construction of a marine fish (Sparus aurata) transcriptome database with a large coverage of protein-coding

Journal Pre-proof

Jo

ur

na

lP

re

-p

ro of

transcripts. BMC Genomics 14, 178. Chen, K., Li, E., Li, T., Xu, C., Wang, X., Lin, H., Qin, J.G., Chen, L., 2015. Transcriptome and molecular pathway analysis of the hepatopancreas in the pacific white shrimp Litopenaeus vannamei under chronic low-salinity stress. PLOS ONE 10, e0131503. Chen S.A., Wang Z.C., Song Y., Chen Y., Lei M.H., 2010. Ecological research about Triplophysa (Hedinichthgs) yarkandensis (Day) in Tarim River. J. Northeast Agr. Univ. 41, 90–93. Chen, S.A., Cheng, Y., Fan, Z.M., Niu, Y.J., Xie, C.X., Lei, H., 2014. Spatial-temporal characteristics of water quality in the upper reaches of the Tarim river in Alar. J. Hydroecol. 35, 15–21. Choi, C.Y., An, K.W., An, M.I., 2008. Molecular characterization and mRNA expression of glutathione peroxidase and glutathione S-transferase during osmotic stress in olive flounder (Paralichthys olivaceus). Comp. Biochem. Phys. A: Mol. Integr. Physiol. 149, 330–337. Davidson, N.M., Oshlack, A., 2014. Corset: enabling differential gene expression analysis for de novoassembled transcriptomes. Genome Biol. 15, 410. de Zelicourt, A., Colcombet, J., Hirt, H., 2016. The role of MAPK Modules and ABA during abiotic stress signaling. Trends Plant Sci. 21, 677–685. Fiol, D.F., Kültz, D., 2007. Osmotic stress sensing and signaling in fishes. FEBS J. 274, 5790– 5798. Gilmour, K.M., Euverman, R.M., Esbaugh, A.J., Kenney, L., Chew, S.F., Ip, Y.K., Perry, S.F., 2007. Mechanisms of acid-base regulation in the African lungfish Protopterus annectens. J. Exp. Biol. 210, 1944–1959. Glaser, R., Kiecolt-Glaser, J.K., 2005. Stress-induced immune dysfunction: implications for health. Nat. Rev. Immunol. 5, 243–251. Götz, S., García-Gómez, J.M., Terol, J., Williams, T.D., Nagaraj, S.H., Nueda, M.J., Robles, M., Talón, M., Dopazo, J., Conesa, A., 2008. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 36, 3420–3435. Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q., Chen, Z., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N., Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotech. 29, 644–652. Hinton, A., Bond, S., Forgac, M., 2009. V-ATPase functions in normal and disease processes. Pflug. Arch. Eur. J. Phy. 457, 589–598. Horng, J.-L., Lin, L.-Y., Huang, C.-J., Katoh, F., Kaneko, T., Hwang, P.-P., 2007. Knockdown of V-ATPase subunit A (atp6v1a) impairs acid secretion and ion balance in zebrafish (Danio rerio). Am. J. Physiol-Reg. I. 292, R2068–R2076. Hwang, P.-P., Lee, T.-H., 2007. New insights into fish ion regulation and mitochondrion-rich cells. Comp. Biochem. Phys. A: Mol. Integr. Physiol. 148, 479–497. Ji, P., Liu, G., Xu, J., Wang, X., Li, J., Zhao, Z., Zhang, X., Zhang, Y., Xu, P., Sun, X., 2012. Characterization of common carp transcriptome: sequencing, De novo assembly, annotation and comparative genomics. PLoS ONE 7, e35152. Kültz, D., Avila, K., 2001. Mitogen-activated protein kinases are in vivo transducers of osmosensory signals in fish gill cells. Comp. Biochem. Phys. B: Biochem. Mol. Biol. 129, 821–829.

Journal Pre-proof

Jo

ur

na

lP

re

-p

ro of

Langmead, B., Salzberg, S.L., 2012. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. Li, B., Dewey, C.N., 2011. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12, 323. Li, E., Wang, S., Li, C., Wang, X., Chen, K., Chen, L., 2014. Transcriptome sequencing revealed the genes and pathways involved in salinity stress of Chinese mitten crab, Eriocheir sinensis. Physiol. Genomics. Livak, K.J., Schmittgen, T.D., 2001. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 25, 402–408. Love, M.I., Huber, W., Anders, S., 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. Lushchak, V.I., 2011. Environmentally induced oxidative stress in aquatic animals. Aquat. Toxicol. 101, 13–30. McNamara, J.C., Rosa, J.C., Greene, L.J., Augusto, A., 2004. Free amino acid pools as effectors of osmostic adjustment in different tissues of the freshwater shrimp macrobrachium olfersii (crustacea, decapoda) during long-term salinity acclimation. Mar. Freshw. Behav. Phy. 37, 193–208. Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A.C., Kanehisa, M., 2007. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 35, W182–W185. Nie, Z., Wu, H., Wei, J., Zhang, X., Ma, Z., 2013. Length-weight relationship and morphological studies in the Kashgarian loach Triplophysa yarkandensis (Day, 1877) from the Tarim River, Tarim River basin, north-west China. Indian J. Fish. 60, 15–19. Tian, F., Liu, S., Shi, J., Qi, H., Zhao, K., Xie, B., 2019. Transcriptomic profiling reveals molecular regulation of seasonal reproduction in Tibetan highland fish, Gymnocypris przewalskii. BMC Genomics 20, 2. Tian, Y., Wen, H., Qi, X., Zhang, X., Li, Y., 2019. Identification of mapk gene family in Lateolabrax maculatus and their expression profiles in response to hypoxia and salinity challenges. Gene 684, 20–29. Tort, L., 2011. Stress and immune modulation in fish. Dev. Comp. Immunol. 35, 1366–1375. Wang, F., Wang, C., Sun, Y., Wang, N., Li, X., Dong, Y., Yao, N., Liu, X., Chen, H., Chen, X., Wang, Z., Li, H., 2016. Overexpression of vacuolar proton pump ATPase (V-H+-ATPase) subunits B, C and H confers tolerance to salt and saline-alkali stresses in transgenic alfalfa (Medicago sativa L.). J. Integer. Agr. 15, 2279–2289. Wang, J., Han, H.D., Zhao, Q.D., Xu J.L., 2013. Study on hydrochemical components of river water in the Tarim River basin, Xinjiang, China. Arid Zone Res. 30, 10–15. Wang, J., Li, B., Meng, Y., Ma, X., Lai, Y., Si, E., Yang, K., Ren, P., Shang, X., Wang, H., 2015. Transcriptomic profiling of the salt-stress response in the halophyte Halogeton glomeratus. BMC Genomics 16, 169. Wilkie, M.P., Wood, C.M., 1996. The adaptations of fish to extremely alkaline environments. Comp. Biochem. Phys. B: Biochem. Mol. Biol.113, 665–673. Wilkie, M.P., Wood, C.M., 1994. The effects of extremely alkaline water (pH 9.5) on rainbow trout gill function and morphology. J. Fish Biol. 45, 87–98. Wilkie, M.P., Wright, P.A., Iwama, G.K., Wood, C.M., 1994. The physiological adaptations of the lahontan cutthroat trout (Oncorhynchus clarki henshawi) following transfer from well water

Journal Pre-proof

Jo

ur

na

lP

re

-p

ro of

to the highly alkaline waters of Pyramid lake, Nevada (pH9.4). Physiol. Zool. 67, 355–380. Xu, J., Ji, P., Wang, B., Zhao, L., Wang, J., Zhao, Z., Zhang, Y., Li, J., Xu, P., Sun, X., 2013a. Transcriptome sequencing and analysis of wild amur ide (Leuciscus waleckii) inhabiting an extreme alkaline-saline lake reveals insights into stress adaptation. PLOS ONE 8, e59703. Xu, J., Li, J.T., Jiang, Y., Peng, W., Yao, Z., Chen, B., Jiang, L., Feng, J., Ji, P., Liu, G., Liu, Z., Tai, R., Dong, C., Sun, Xiaoqing, Zhao, Z.-X., Zhang, Y., Wang, J., Li, S., Zhao, Y., Yang, J., Sun, Xiaowen, Xu, P., 2017. Genomic basis of adaptive evolution: the survival of amur ide (leuciscus waleckii) in an extremely alkaline environment. Mol. Biol. Evol. 34, 145–159. Xu, J., Li, Q., Xu, L., Wang, S., Jiang, Y., Zhao, Z., Zhang, Y., Li, J., Dong, C., Xu, P., Sun, X., 2013b. Gene expression changes leading extreme alkaline tolerance in Amur ide (Leuciscus waleckii) inhabiting soda lake. BMC Genomics 14, 682. Yao, Z.L., Wang, H., Chen, L., Zhou, K., Ying, C.Q., Lai, Q.F., 2012. Transcriptomic profiles of Japanese medaka (Oryzias latipes) in response to alkalinity stress. Genet. Mol. Res. 11, 2200–2246. Young, M.D., Wakefield, M.J., Smyth, G.K., Oshlack, A., 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol. 11, R14. Yin, M.C., 1995. Fish ecology. China Agriculture Press, Beijing. Zhang, W., Liu, K., Tan, B., Liu, H., Dong, X., Yang, Q., Chi, S., Zhang, S., Wang, H., 2019. Transcriptome, enzyme activity and histopathology analysis reveal the effects of dietary carbohydrate on glycometabolism in juvenile largemouth bass, Micropterus salmoides. Aquaculture 504, 39–51. Zhao, Y., Wang, J., Thammaratsuntorn, J., Wu, J.W., Wei, J.H., Wang, Y., Xu, J.W., Zhao, J.L., 2015. Comparative transcriptome analysis of nile tilapia (Oreochromis niloticus) in response to alkalinity stress. Genet. Mol. Res. 14, 17916–17926. Zhou, X., Naguro, I., Ichijo, H., Watanabe, K., 2016. Mitogen-activated protein kinases as key players in osmotic stress signaling. B.B.A.-Gen. Subjects 1860, 2037–2052.

Journal Pre-proof Figure captions: Fig.1. Overview of unigene length distributions of the Triplophysa yarkandensis transcriptome. The x-axis indicates unigene size and the y-axis indicates the number of unigenes with different lengths. Fig.2. Gene Ontology classification of assembled unigenes. Fig.3. Histogram representing the KOG classification of 50598 unigenes. The letters on the x-axis indicate different KOG categories. Fig.4. KEGG pathway enrichment of all annotated unigenes.

ro of

Fig.5. Differential expression analyses of unigenes between fish subject to salinity/control (A), alkalinity/control (B), and alkalinity/salinity (C) treatments. Each dot represents one unigene, red

-p

dots represent the significantly upregulated unigenes, and green dots represent the significantly downregulated unigenes. Blue dots represent unigenes that are not significantly differentially

re

expressed.

lP

Fig.6. Pairwise comparisons of differentially expressed unigenes from fish subject to salinity/control, alkalinity/control, and alkalinity/salinity treatments.

na

Fig.7. Heatmap analysis of hierarchical clustering of differentially expressed unigenes from fish subject to salinity, alkalinity, and control treatments. Each column represents a salinity, alkalinity,

ur

or control sample and each row represents a unigene. Different colors indicate differences in expression level. Positive numbers indicate upregulation and negative numbers indicate

Jo

downregulation.

Fig.8. Validation of differentially expressed unigenes by qRT-PCR. The green columns indicate the results of RNA-Seq and the red columns represent the results of qRT-PCR.

Journal Pre-proof

Table 1. The absolute, relative and instantaneous growth rate after 30 days of different concentration of salinity and alkalinity exposure in T. yarkandensis Salinity (‰)

Alkalinity (mmol/L)

Control 2.0‰ a

0.0438±0.0001

4.0‰ a

6.0‰

0.0168±0.0003

b

6.9 mmol/L

0.0160±0.0024

b

0.0605±0.0014c

0.0419±0.00001

Relative growth rate

10.014±1.9781a

10.00±2.2130a

9.8124±0.8181b

9.2162±0.3547b

10.1048±1.8028a

9.9801±0.5629b

9.3288±1.8345b

Instantaneous growth rate

0.1602±0.00005a

0.1352±0.00023a

0.0045±0.0001b

0.0038±0.0002b

0.01206±0.0012c

0.0044±0.0023b

0.0151±0.0001c

Note: The different letters on the parameters in one row mean significant difference (P ＜ 0.05)．

o r p

l a

o J

n r u

f o

e

0.0029±0.0001

12.9 mmol/L b

Absolute growth rate

r P

0.0727±0.0032

9.9 mmol/L a

Journal Pre-proof

Raw Reads

Clean reads

Percentage kept (%)

Clean bases

Q20(%)

Q30(%)

GC(%)

con_1

62788100

59939302

95.46

8.99G

97.49

92.79

44.20

con_2

53969806

49749842

92.18

7.46G

95.61

89.06

43.89

con_3

56981682

54993590

96.51

8.25G

97.37

92.52

43.69

A6.9_1

55294032

51737494

93.57

7.76G

96.9

91.52

41.79

A6.9_2

59958414

57083816

95.21

8.56G

97.62

93.03

41.82

A6.9_3

56958394

54847934

96.29

8.23G

97.57

92.95

41.61

A9.9_1

61126418

58089710

95.03

8.71G

97.73

93.31

42.74

A9.9_2

56972204

54187558

95.11

8.13G

97.55

92.86

42.49

A9.9_3

56008932

53359258

95.27

8.00G

92.89

42.82

A12.9_1

45761802

41881380

91.52

6.28G

97.23

92.32

41.67

A12.9_2

56826052

52854164

93.01

7.93G

97.61

93.07

42.14

A12.9_3

64290386

60095218

93.47

9.01G

96.85

91.47

42.51

S2_1

42807696

41349814

96.59

ro of

97.55

6.20G

98.28

94.62

43.09

S2_2

58670078

57204764

97.50

8.58G

98.26

94.56

42.48

S2_3

56186980

53859666

95.86

8.08G

98.23

94.52

42.25

S4_1

58365672

55872240

95.73

8.38G

98.22

94.45

42.20

S4_2

55834148

53881960

96.50

8.08G

98.15

94.28

41.98

S4_3

48362816

46221814

95.57

6.93G

97.41

92.56

41.53

S6_1

50840188

48028872

94.47

7.20G

97.29

92.36

42.66

S6_2

67971244

64301866

94.60

9.65G

97.04

91.83

42.64

S6_3

56248968

53908702

8.09G

97.69

93.19

42.72

lP

re

-p

Sample

na

Table 2. Summary of RNA-seq data.

Jo

ur

95.84

Journal Pre-proof

Table 3. Summary of assembled results of RNA-seq. Length (bp)

Min Length

201

Mean Length

1703

Median Length

1043

Max Length

19318

N50

2923

N90

719

Total Nucleotides

301806400

Jo

ur

na

lP

re

-p

ro of

Item

Journal Pre-proof

Table 4. Statistics of success rate of gene annotation. Percentage (%)

Annotated in NR

100243

56.55

Annotated in NT

111884

63.11

Annotated in KO

14679

8.28

Annotated in SwissProt

103326

58.29

Annotated in PFAM

97651

55.09

Annotated in GO

97651

55.09

Annotated in KOG

50598

28.54

Annotated in all Databases

6716

3.79

Annotated in at least one Database

140923

79.50

Total Unigenes

177271

100

Jo

ur

na

lP

re

-p

ro of

Number of Genes

Journal Pre-proof

Table 5. The different expressed unigenes in different salinity concentration gradients control

2.0‰

4.0‰

6.0‰

control

--

up: 3876 down: 6378

2.0‰

--

--

up: 5850 down: 10961 up: 259 down: 181

4.0‰

--

--

--

up: 7452 down: 6552 up: 181 down: 439 up: 327 down: 936

6.0%

--

--

--

Jo

ur

na

lP

re

-p

ro of

--

Journal Pre-proof

Table 6. The different expressed unigenes in different alkalinity concentration gradients control

6.9 mmol/L

9.9 mmol/L

12.9 mmol/L

control

--

up: 14126 down: 15046

6.9 mmol/L

--

--

up: 10695 down: 8287 up: 792 down: 2422

9.9 mmol/L

--

--

--

up: 15270 down: 15505 up: 11324 down: 9788 up: 7596 down: 6306

12.9 mmol/L

--

--

--

Jo

ur

na

lP

re

-p

ro of

--

Journal Pre-proof

Table 7. Up and down regulated pathways in the gill of salinity treated T. yarkandensis. Unigene

P value

Up regulated

Transcriptional misregulation in cancer Protein export MAPK signaling pathway Collecting duct acid secretion Fatty acid elongation

27 7 17 8 7

0.00098018 0.00603407 0.02083553 0.02430592 0.0362882

Down regulated

Amphetamine addiction Pathways in cancer Tight junction Dopaminergic synapse Cocaine addiction Viral carcinogenesis Ras signaling pathway NOD-like receptor signaling pathway Endocytosis Hippo signaling pathway Alcoholism Other types of O-glycan biosynthesis

12 47 23 18 8 27 26 12 36 19 14 7

0.00170881 0.00434816 0.00775285 0.01099711 0.01292317 0.02052732 0.02397617 0.02712124 0.03080031 0.03409204 0.04366263 0.04666189

Jo

ur

na

lP

re

-p

ro of

Pathway

Journal Pre-proof Table 8. Up and down regulated pathways in the gill of alkalinity treated T. yarkandensis. Unigene

P value

Up regulated

Transcriptional misregulation in cancer Epithelial cell signaling in Helicobacter pylori infection Synaptic vesicle cycle Vibrio cholerae infection Protein export Spliceosome Oxidative phosphorylation Collecting duct acid secretion Rheumatoid arthritis Fatty acid elongation Glyoxylate and dicarboxylate metabolism Circadian rhythm Valine, leucine and isoleucine degradation Fat digestion and absorption Ribosome biogenesis in eukaryotes

32 23 17 20 8 28 31 10 18 9 11 6 15 6 14

0.00150513 0.00392623 0.00662983 0.00668037 0.00726663 0.01012418 0.01561612 0.0165923 0.02057665 0.02104461 0.0210955 0.02234435 0.03143308 0.03394339 0.04203459

Down regulated

Herpes simplex infection Viral carcinogenesis Ras signaling pathway B cell receptor signaling pathway RIG-I-like receptor signaling pathway NOD-like receptor signaling pathway Inflammatory bowel disease (IBD) Graft-versus-host disease Allograft rejection NF-kappa B signaling pathway Pancreatic cancer Pathways in cancer Viral myocarditis

34 33 32 15 13 14 11 8 8 15 14 48 13

0.00800874 0.00839901 0.00917293 0.0128086 0.02002194 0.0213271 0.02224531 0.02591424 0.02836096 0.03286906 0.03982375 0.0421838 0.04349501

Jo

ur

na

lP

re

-p

ro of

Pathway

Journal Pre-proof

Jo

ur

na

lP

re

-p

ro of

Highlights: 1. We obtained the transcriptome profile of Triplophysa yarkandensis for the first time. 2. We studied the responding to different levels of salinity and alkalinity stress of T. yarkandensis on the transcriptome level. 3. Both salinity and alkalinity stress downregulated immunity or disease related genes and pathways.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8