Three adults of T. yarkandensis were collected from the reach of Shule river basin (96.71 E, 40.57 N) at an altitude of 1,300 m. A total of fifteen tissues were collected from the three fish specimens including heart, liver, brain, spleen and kidney. These tissues were stored in ultra-low temperature liquid nitrogen at −80°C for standby. Liquid nitrogen was added to grind the tissues into powder, which allowed to extracted total RNA from the mixed tissues using the standard Trizol kit (Promega SV Total RNA Isolation system), and detect the quality of total RNA using agarose gel electrophoresis. This process allowed, the detection of the degradation degree, purity, quality, concentration and integrity of total RNA by Nanodrop spectrophotometer (Thermo Fisher Scientific, MA, United State), Qubit fluorescence quantitative analyzer (Thermo Fisher Scientific, MA, Unite State) and Agilent 2,100 biological analyzer (Agilent, CA, United State), to ensure that the samples met the requirements for transcriptome sequencing. After detection, the mRNA within the total RNA was broken into specific short fragments, and the double-terminal cDNA library of the transcriptome was amplified by RT-PCR. Based on sequencing by synthesis (SBS) technology, the cDNA library was sequenced by Illumina Hiseq high-throughput sequencing platform (Illumina, CA, United State), which could produce a large number of high-quality reads (raw data). After detecting quality and filtering, clean reads were produced.

The clean reads were employed for de novo assembly using Trinity (v2.5.1) (Grabherr et al., 2011). The CD-HIT-EST program was utilized to remove the redundant sequences, and the parameters were set as the similarity of 0.95 and the sequence length of 10 bases (Li and Godzik, 2006). Annotation information were finally obtained by selecting BLAST (Altschul et al., 1997) parameter E-value ≤ 1e-5 and Hmmer (Finn et al., 2014) parameter E-value ≤ 1e-10. The annotation information of unigenes were derived from nine databases, including NCBI non-redundant protein sequence (NR) (Deng et al., 2006), Swiss prot (Rolf et al., 2004), Clusters of Orthologous Groups (COG) (Tatusov et al., 2000), EuKaryotic Orthologous Groups (KOG) (Koonin et al., 2004), Evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) (Jaime et al., 2016), Gene Ontology (GO) (Ashburner et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Minoru et al., 2004), Protein family (Pfam) (Finn et al., 2014), and Translation of EMBL (TrEMBL) (Amos and Rolf, 2000).

Bowtie (Langmead, 2009) was applied to compare the sequenced reads with the unigene library, and the expression level was estimated in combination with RSEM (Dewey and Li, 2011). The expression abundance of corresponding unigene was expressed by FPKM (Fragments Per Kilobase of transcript per Million mapped reads) value (Trapnell et al., 2010). Differential expression analysis was performed using DeSeq2 (v1.6.3) (Love et al., 2014). Differential expression results were obtained with FDR (False Discovery Rate) > .01 and the difference multiple FC (fold change) ≥ 2 were considered as the screening standards. The differentially expressed genes of five tissues were analyzed by GO functional annotation and KEGG pathway enrichment analysis.

In Shule river, the cations were predominantly Ca²⁺, Mg²⁺, Na⁺ and the anions were mainly HCO₃ ⁻, SO₄ ^2-, Cl⁻, and the pH was 7–8 (Zhou et al., 2004). The water was weakly saline and alkaline. Studies have indicated that T. yarkandensis had a certain tolerance under precondition of salinity and alkalinity stress (Yao et al., 2018). In the current study, six genes were identified which participated in the regulation process of salinity-alkalinity adaptation, these genes were mapk1 (mitogen-activated protein kinase 1), abcc1 (ATP-binding cassette sub-family C member 1), gpx1 (glutathione peroxidase 1), gpx4 (glutathione peroxidase 4), cat (catalase) and aqp1 (aquaporin-1) (Li et al., 2016; Chen et al., 2020; Zhang, 2020; Wang et al., 2022c). From the differential genes obtained in this study, the six genes were searched and the differences between tissues were analyzed. Additionally, according to the results of KEGG enrichment, the pathway map involving these genes was obtained. Finally, mapk1 gene was selected for bioinformatics analysis, including sequence characteristics and important amino acid sites. Firstly, FGENESH (http://www.softberry.com/berry.phtml?topic=fgenesh&group = programsandsubgroup = gfind) was used to analyze its sequence characteristics. Moreover, a search at the BLAST tool at NCBI was conducted for similar sequences of mapk1 gene, and important amino acid sites were compared by using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/).

Water temperature is one of the crucial environmental factors not only affecting almost all life activities of fish, but also affecting the physiological and biochemical reactions of fish. Fish can gradually adapt low temperature environment by regulating the concentration of enzymes in different metabolic pathways. A previous report not only indicated that low temperature heightened lipid synthesis in the body of Trilophysa bleekeri, but also showed the content of polyunsaturated fatty acids in the body was up-regulated in the low temperature group (Yuan, 2021). Through fatty acid metabolism, the content of long-chain unsaturated fatty acids in phospholipid of cell membrane increased and the cell membrane would maintain stability and fluidity. In the current study, six genes was identified, which are involved in fatty acid metabolism pathway, these genes were acox1 (acyl-CoA oxidase 1, palmitoyl), acox3 (acyl-CoA oxidase 3, pristanoyl), acaa2 (acetyl-CoA acyltransferase 2), elovl4 (elongation of very long chain fatty acids elongase gene 4), acaca (acetyl-CoA carboxylases alpha) and fasn (fatty acid synthase) (Xu et al., 2020; Yang et al., 2020; Wang et al., 2022d; Liu et al., 2022). From the differential genes obtained in this study, the six genes were searched and the differences among tissues were analyzed. Moreover, according to the results of KEGG enrichment, the pathway map involving these genes was obtained. Finally, acox1 gene was selected for bioinformatics analysis, including sequence characteristics and important amino acid sites. Firstly, FGENESH (http://www.softberry.com/berry.phtml ?topic = fgeneshandgroup = programsandsubgroup = gfind) was used to analyze its sequence characteristics. Moreover, a search at the BLAST tool at NCBI was conducted for similar sequences of acox1 gene, and important amino acid sites were compared by using Clustal Omega (https://www.ebi.ac.uk/Tools/msa/clustalo/).

The cDNA library was sequenced using Illumina Hiseq high-throughput sequencing platform. After sequencing quality control, a total of 102.11 Gb clean data were obtained. Overall, the percentage of Q30 bases of all products exceeded 91.28%, and the GC content exceeded 44.90%. The sample sequencing data evaluation statistics information was shown in the Supplementary Table S1. A total of 100,112 unigenes were obtained by assembly, the N50 of transcript and unigene were 3,938 bp and 2,579 bp, with high assembly integrity (Supplementary Table S2, Supplementary Figure S1).

In the current study, 23,008 unigenes were annotated to nine databases, including 5,207, 19,535, 19,637, 15,404, 18,696, 13,941, 19,259, 18,421 and 19,270 unigenes in COG, GO, KEGG, KOG, Pfam, Swiss prot, TrEMBL, eggNOG and NR databases, respectively (Supplementary Table S3). Based on the BLASTx similarity analysis against the NR database, T. yarkandensis unigenes had the highest number of hits compared to other fish species such as to Carassius auratus (45.37%), followed by Anabarilius grahami (18.62%), Danio rerio (13.12%) (Figure 1A). A total of 19,637 assembled unigenes were assigned to 48 terms for GO classification in the biological process (25 terms), cellular component (6 terms), and molecular function (17 terms) categories (Figure 1B). Cellular anatomical entity (15,116 unigenes) showed the most annotated unigenes, followed by cellular process (13,879 unigenes), binding (12,320 unigenes), and intracellular (10,866 unigenes). In KOG database (Figure 1C), unigenes have been divided into 25 functional categories, and among them, the largest was signal transduction mechanisms, with 3,151 unigenes, accounting for 20.46%, followed by nothing but general function prediction (2,858 unigenes, 18.55%), posttranslational modification, protein turnover, chaperones (1,435 unigenes, 9.32%). No unknown genes were found from the KOG database annotation.

We compared DEGs in five tissues: brain, heart, liver, kidney and spleen. Results of differentially expressed genes between each two tissues were shown in Supplementary Table S4. TYbrain_vs_TYliver screened the largest differentially expressed genes number, with a 9,814 DEGs, of which 3,818 genes were up-regulated and 5,996 genes were down-regulated (Figure 2A). Hierarchical clustering analysis was performed on the selected differentially expressed genes between brain and liver, and the genes with the same or similar expression behavior were clustered (Figure 2B). GO annotation results and KEGG enrichment of DEGs between brain and liver were shown in Figures 2C, D. The three biological duplicate samples of the same tissue could cluster, indicating that the biological duplicate of the samples used in this study was favorable. The GO annotation results of ten groups were similar. Most of the terms included “cellular process”, “single-organism process”, “biological regulation”, “metabolic process”, “cell”, “cell part”, “binding” and “catalytic activity”. In KEGG pathway enrichment analysis, DEGs were enriched including MAPK signaling pathway, calcium signaling pathway and neuroaction ligand-receptor interaction.

19,637 unigenes in KEGG database correspond to 330 pathways, and the top 20 pathways with the largest number of sequences were shown in the Supplementary Table S5. Unigenes with the largest three pathways annotations (neuroactive ligand-receptor interaction, calcium signaling pathway and MAPK signaling pathway) were associated with environmental information processing. In addition, eight of the top 20 pathways were correlated with signaling molecules and interaction, signal transduction. A previous study showed that different selenium sources could regulate the adaptation of yellow catfish to low temperature stress, neuroactive ligand-receptor interaction pathway significant enrichment (Hu et al., 2021). In the current study, the number of unigenes annotated in the neuroactive ligand-receptor interaction pathway was the largest, which was speculated to be associated with the low temperature environment adaptation of T. yarkandensis.

In the current study, six candidate genes in the salinity-alkalinity adaptation: mapk1, abcc1, cat, gpx1, gpx4 and aqp1 were screened, and the expression differences of six genes among five tissues of T. yarkandensis were analyzed. The experimental results indicated that these genes were differentially expressed among the five tissues. As shown in Table 1, the expression of cat, abcc1 and gpx4 genes in liver was higher than that in brain, heart, kidney and spleen, while the expression of gpx1 and aqp1 genes in liver, heart, kidney and spleen was higher than that in brain. In addition, the expression of mapk1 genes in liver and kidney was higher. A previous research revealed that salinity stress could affect the changes of metabolic pathways (Liu et al., 2020). The six genes involved in the pathway were shown in Figure 3. Six genes in five tissues of T. yarkandensis were differentially expressed and all were highly expressed in the liver, which it was suspected to be related to the metabolism. Finally, mapk1 gene for bioinformatics analysis were selected, and the result showed that the sequence length was 2,424 bp, and there was one coding region of the predicted gene structure. Besides, divergent amino acid changes (Val-Ile at 321, 325) by comparing amino acid sequences with seven other fish were found (Figure 4A).

In this study, six genes involved in fatty acid metabolism pathway, elovl4, acaca, fasn, acaa2, acox1, acox3, were selected to analyze the expression differences of six genes among liver and other tissues of T. yarkandensis. The results demonstrated that the expression of elovl4, acaca, fasn, acox1 and acox3 in liver were higher than that in brain, heart, kidney and spleen. The expression of acaa2 gene in liver were higher than that in brain, heart and spleen. The expression of six genes in fatty acid metabolism pathway differed among different tissues (Table 2). These genes involved in fatty acid metabolism pathway were shown in Figure 5. The expression of six genes in liver were higher than that in other tissues, indicating that T. yarkandensis had relatively fast fatty acid metabolism, which might be related to its low temperature adaptation. Finally, acox1 gene for bioinformatics analysis were selected, and the result showed that the sequence length was 5,277 bp, and there were two coding regions of the predicted gene structure. Besides, divergent amino acid changes (His-Gln at 402) by comparing amino acid sequences with seven other fish were found (Figure 4B).