Illumina sequenced read data were deposited in the NCBI Sequence Read Archive (accession \# SRP071094).

Introduction {#sec001}
============

Fish survival is affected by abiotic stresses such as extreme temperatures, salinity, hypoxia and chemical toxicity \[[@pone.0186433.ref001],[@pone.0186433.ref002]\], but water temperature is the chief environmental determinant of development, growth, reproduction, behavior, metabolism and geographical distribution \[[@pone.0186433.ref002]--[@pone.0186433.ref004]\]. As ectotherms, most fish species encounter large temperature variations between day and night and among different seasons, which leads to deleterious consequences when the environmental temperature exceeds species-specific thermal tolerances \[[@pone.0186433.ref005]\]. Heat stress could influence fish activity and reduce aquaculture production \[[@pone.0186433.ref002]\], especially as a consequence of global warming \[[@pone.0186433.ref006]\] which modifies survival, fertility and large-scale migration \[[@pone.0186433.ref007],[@pone.0186433.ref008]\]. Animals respond to such environmental changes variously \[[@pone.0186433.ref006],[@pone.0186433.ref009]--[@pone.0186433.ref012]\], so it is of both biological and ecological interests to understand the mechanisms underlying the response of fish to heat stress.

Schizothoracine fishes are dominant fish species of the Tibetan Plateau (TP), which has more than 60 species distributed throughout TP altitudes \[[@pone.0186433.ref013],[@pone.0186433.ref014]\]. TP aquatic environments are big challenges for aquatic animals and inhospitable due to long-term low temperature, hypoxia, and high salinity \[[@pone.0186433.ref014]--[@pone.0186433.ref016]\]. Unlike other Schizothoracine fishes, Tibetan naked carp *Gymnocypris przewalskii* is one of the best-studied Schizothoracine species, making it a model aquatic animal organism for studying speciation and environmental adaptation \[[@pone.0186433.ref016]--[@pone.0186433.ref020]\]. Lake Qinghai is the largest salt lake of China. It is formed during the creation of the TP, which produces the hostile aquatic environment \[[@pone.0186433.ref013]\]. Consequently, many ancestral endemic fish species were extinct with the exclusive exception of the *G*. *przewalskii*, \[[@pone.0186433.ref013],[@pone.0186433.ref021]\] which has a remarkable adaptive capacity that contributed to its adaptation to the Lake Qinghai \[[@pone.0186433.ref016],[@pone.0186433.ref022]\]. In addition, previous data show that *G*. *przewalskii* could not tolerate the high water temperature exceeded 30°C, which triggered high mortality in field condition \[[@pone.0186433.ref014]\]. While, another evidence showed that this cold-water fish species could survive under severe heat exposure (30°C) and the optimum farming temperature of 24°C contributed to a positive impact on the growth rate in the aquaculture industry \[[@pone.0186433.ref019]\]. Nevertheless, this controversial issue on the thermal acclimation of *G*. *przewalskii* has not been experimentally investigated, and the regulation of mechanism responsible for heat stress in Tibetan Schizothoracine fishes also has been rarely revealed.

Many studies have focused on thermal acclimation in animals \[[@pone.0186433.ref020]--[@pone.0186433.ref025]\], and suggest that larval and later developmental stages are more sensitive to temperature change than adults \[[@pone.0186433.ref026],[@pone.0186433.ref027]\], which allows study of fish responses to heat stress. RNA sequencing technology has been successfully used to reveal the heat (cold) stress response mechanism in larvae of several fish species \[[@pone.0186433.ref005],[@pone.0186433.ref028],[@pone.0186433.ref029]\]. The microRNAs (miRNAs) are a family of endogenous RNAs that regulated gene expression in a sequence-specific manner. The miRNAs are usually 22 nucleotides long with a seed region (2--8 bp) that binds to the 3' untranslated region (UTR) of a targeted mRNA, signaling it for the degradation. The miRNAs have diverse expression patterns and regulate developmental and physiological processes \[[@pone.0186433.ref030]\]. Previous studies have confirmed the mechanisms of miRNAs in animal responses to both biotic and abiotic stresses \[[@pone.0186433.ref031]--[@pone.0186433.ref033]\], while our understanding of the molecular mechanism of thermal acclimation regulated by both mRNA and miRNA in fish remains limited, especially in highland fish species. The potential thermal acclimation of the cold-water fish *G*. *przewalskii* has always been controversial \[[@pone.0186433.ref014],[@pone.0186433.ref023]\]. Here we hypothesize that *G*. *przewalskii* could survive and resist heat stress against the controversial issue. Thus, we used the *G*. *przewalskii* larvae as experiment materials, and assessed mRNA and miRNA transcriptomes to investigate the regulation of thermal acclimation in *G*. *przewalskii*.

Materials and methods {#sec002}
=====================

Animals heat treatment and statistical analysis {#sec003}
-----------------------------------------------

All animal experiments were approved by the Animal Care and Use Committees of the Northwest Institute of Plateau Biology, Chinese Academy of Sciences and the Agriculture Department of Qinghai Province, China. The present study is based on wild-caught *G*. *przewalskii* that were collected in the Lake Qinghai in May 2014 and provided by the National Rescue Center of Qinghai Lake Naked Carp (Xining, Qinghai, China). All the sample were subsequently reared under laboratory conditions. The artificial insemination is based on a wild-caught female and a male of *G*. *przewalskii*. Fertilized eggs were collected and transferred to another re-circulating aerated freshwater system (12 h day/12 h night, 16°C). At 192 hour post fertilization (hpf), all the larvae of *G*. *przewalskii* were hatched out (n = 600, [Fig 1A](#pone.0186433.g001){ref-type="fig"}). This developmental stage allows a reseasonable period of post-hatching and no need for feeding, which provided suitable materials for environmental stress experiment \[[@pone.0186433.ref005],[@pone.0186433.ref029]\]. Then fish larvae at 192 hpf were separated into two groups (300 larvae per group) and immediately transferred into six dishes (100 larvae per dish). Control group (CT) were continued to maintain at 16°C. Heat stress groups (HS) were exposed to 24°C of 12 hours with heating rate of \~ 0.8°C/h ([Fig 1A](#pone.0186433.g001){ref-type="fig"}). No embryonic mortality was observed before the treatment. After heat exposure (216 hpf, [Fig 1A](#pone.0186433.g001){ref-type="fig"}), fish were anaesthetized with MS-222 (2%, dipping treatment, Sigma, St. Louis, MO). To generated three biological replicates, equal number of fish larvae from CT and HS groups were collected thrice. All the larvae samples were immediately frozen in liquid nitrogen prior to total RNA extraction for RNA library preparation followed by subsequent mRNA and small RNA sequencing.

![Establishment of heat stress resistance in *G*. *przewalskii* after exposure to mid-high temperature under further severe heat stress.\
(a) Schematic diagram of heat exposure experiment. *G*. *przewalskii* embryos were incubated at 16°C from fertilization to 192 hpf (hatching out). Larvae at 192 hpf were exposed to heat stress at 24°C (HS group) and the controls (CT group) were maintained at 16°C. Samples for RNA-seq and miRNA-seq were collected at 216 hpf. After a comparative treatment for 12 h, larvae from both groups were further exposed to severe heat stress (30°C) for 48 h. The heating rate is \~ 0.8°C per hour. Images of larvae were taken by using stereomicroscope from Zeiss SteREO DiscoveryV12 with a color CCD camera. Death rates of larvae from both groups were recorded respectively. (b) Bar plot indicates the death rates of larvae from HS and CT groups at different time. HS and CT groups were colored by orange and blue, respectively. Data was shown as mean ± standard error of the mean (n = 3). "\*\*" above error bars indicate *p* \< 0.01. (c) Samples from HS and CT groups were exposed to 30°C for 48 h. Image of larvae were also taken by using Zeiss SteREO DiscoveryV12. Red arrowheads indicate representative malformed and dead larvae.](pone.0186433.g001){#pone.0186433.g001}

Then, larval fish from both HS and CT groups were exposed to severe high temperature (30°C) of 48 hours with a same heating rate of \~ 0.8°C/h ([Fig 1A](#pone.0186433.g001){ref-type="fig"}). The morphological observation of larvae was recorded with a stereomicroscope system (Zeiss SteREO DiscoveryV12, Germany). In addition, the death rate of larvae in both groups was recorded separately at 6, 12, 24, 36, and 48 h. Statistical analysis was performed by SAS v9.1.3 (SAS Institute Inc., USA). Significant differences between samples were analyzed by T test. Differences were accepted as significant when *p* \< 0.05. All data were shown as mean ± standard error of the mean (SEM).

RNA extraction and library construction {#sec004}
---------------------------------------

Total RNA of each sample (n = 2 per group) was isolated using TRIzol (Invitrogen, Carlsbad, CA) and purified using a TruSeq RNA Sample Prep Kit V2 (Illumina Inc., San Diego, CA) according to the manufacturer\'s instructions. The quality and quantity of total RNA samples were assessed with an Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA) and SDS-PAGE. Approximately 35 μg of total RNA were used for mRNA transcriptome library preparation (4 libraries), and were sequenced with an Illumina HiSeq™ 2500 instrument.

Bioinformatics analysis of RNA-seq data {#sec005}
---------------------------------------

RNA-seq raw reads from each library were preprocessed to filter clipped adapter sequences, contaminated sequences, low-quality reads (Q \< 20). High-quality reads were assembled to contigs using the *Trinity* program ([http://trinityrnaseq.sf.net](http://trinityrnaseq.sf.net/)) with default parameters. Contigs of each assembly were performed with TGICL to produce long and complete unigenes \[[@pone.0186433.ref034]\], which were generated with a minimum overlap length of 200 bp. The assembled unigene sequences were aligned with a Blast-X search (cut-off E-value of 1e^-6^) in the public NCBI non-redundant (NR), Swiss-Prot protein and eggNOG databases (<http://eggnog.embl.de/version_3.0/>). Then the unigene was assigned by protein identity with the highest sequence similarity, which was used for functional annotation using the Blast2GO program \[[@pone.0186433.ref035]\]. Finally pathway assignments were generated using the KEGG database, Kyoto Encyclopedia of Genes and Genomes \[[@pone.0186433.ref036]\].

Gene expression analysis {#sec006}
------------------------

Gene expression values were calculated as reads aligned to gene per kilo base of exon per million mapped reads (RPKM) \[[@pone.0186433.ref037]\]. Statistical comparison between two libraries was conducted with *DESeq* package in R software \[[@pone.0186433.ref038]\]. FDR (false discovery rate) \< 0.05 was used as the threshold of *p*-value in multiple tests to measure significant gene expression difference \[[@pone.0186433.ref039]\]. Genes were only considered differently expressed in a given library when the *p* \< 0.05 and a greater than two-fold change (absolute value of log~2~ ratio \> 1) in expression across libraries was observed.

Small RNA sequencing and differential expression analysis {#sec007}
---------------------------------------------------------

Total RNA from larvae (n = 5 per group) was extracted with miRNeasy Mini kit (Qiagen, Germany), and then used for small RNA library construction (4libraries), finally sequenced on an Illumina HiSeq™ 2500 instrument. Notably, the small RNA sequencing was independent of the mRNA transcriptome sequencing and not shared a single lane on the Illumina sequencer.

Similar to the process of RNA-seq, raw reads of small RNA-seq were treated by trimming adapters and removing poor-quality reads. Other RNAs (rRNA, tRNA, snRNA and snoRNA) were removed by blasting against Genbank database ([http://blast.ncbi.nlm.nih.gov](http://blast.ncbi.nlm.nih.gov/)) and Rfam database (<http://sanger.ac.uk/software/Rfam>). Clean reads ranging from 15 to 30 nt were searched against known precursor/mature miRNAs in miRBase 19.0 (<http://www.mirbase.org/>).

Differentially expressed miRNA (DEM) between CT and HS groups were evaluated on a log2-ratio plot. MiRNA expression was first normalized based on total sequencing reads from each library. Criteria for differentially expressed miRNAs were as follows: absolute value of log2 (HS/CT) ≥ 1.0 and having ≥ 1,000 normalized reads in at least one sample.

miRNA target prediction and regulatory network construction {#sec008}
-----------------------------------------------------------

*In silico* analysis of DEM targets were predicted in DEG using miRanda 4.0 algorithm ([microrna.sanger.ac.uk/targets/v4](http://microrna.sanger.ac.uk/targets/v4)) and *TargetScan* ([www.targetscan.org](http://www.targetscan.org/)), based on the complementary region between miRNAs and mRNAs and the thermodynamic stability of the miRNA-mRNA duplex. All the miRNA target mRNAs were calculated and clustered by GO terms and KEGG pathway annotations, respectively. Regulatory network predictions for mRNA and miRNA were performed using Cytoscape 3.3.0 (<http://www.cytoscape.org/>).

Functional characterization of rapid heat stress response features {#sec009}
------------------------------------------------------------------

Real-time quantitative PCR (RT-qPCR) was used to evaluate the expression of selected DEGs and DEMs. Specific primers were designed with Beacon Designer 8.0 software ([S1 Table](#pone.0186433.s004){ref-type="supplementary-material"}). Reversed cDNA and miRNA were synthesized with PrimeScript^TM^ RT Kit (Takara, Dalian, China) and miRcute miRNA Isolation Kit (Tiangen, Beijing, China) following the manufacturer\'s protocols, respectively. RT-qPCR experiments were performed with Power SYBR Green PCR Master Mix Kit (Applied Biosystems, CA, USA) for mRNA and miRcute Plus miRNA qPCR Detection Kit (Tiangen, Beijing, China) for miRNA on ABI ViiA™7 (ABI, CA, USA) instrument. To normalize the expression values, three genes (*β-actin*, *GAPDH* and *EF-1a*) and three miRNA (miR-22a, miR-152 and miR-25) were selected and used as housekeeping control for mRNA and miRNA expression assessments ([S1 Table](#pone.0186433.s004){ref-type="supplementary-material"}) \[[@pone.0186433.ref040],[@pone.0186433.ref041]\]. All samples were analyzed in triplicate as technical replicates and fold changes were calculated using comparative CT method (also known as the 2^−ΔΔCt^ method). Statistical analysis was performed by SAS v9.1.3 (SAS Institute Inc., USA). Significant differences between samples were analyzed by T test. Differences were accepted as significant when *p* \< 0.05. All data are shown as mean ± standard error of the mean (SEM).

Results {#sec010}
=======

Effect of rapid heat stress on Tibetan naked carp larval development {#sec011}
--------------------------------------------------------------------

After larvae exposure to severe high temperature at 30°C ([Fig 1A](#pone.0186433.g001){ref-type="fig"}), the death rates of HS and CT groups had diverged. Statistical analysis showed that the death rates of HS are significantly lower than those of CT at 12, 24, 36 and 48 h (*p* \< 0.01 in all time points, [Fig 1B](#pone.0186433.g001){ref-type="fig"}). The death rates following severe heat exposure for 36 and 48 h reached 93.2% and 100% in the control, while were 33.1% and 54.4% in the HS group. In addition, microscopic observation found that most of larvae in CT group were dead and displayed malformation after exposure to 30°C for 36 h, but the survival larvae in HS group were normal in morphology ([Fig 1C](#pone.0186433.g001){ref-type="fig"}). Thus, these finding suggests that Tibetan naked carp larvae develop a preliminary resistance against severe heat stress (30°C) after exposure to 24°C.

Differential expression of mRNAs and miRNAs {#sec012}
-------------------------------------------

From six larvae RNA-seq datasets (three biological replicates each group) of Tibetan naked carp larvae at 16°C (CT) and 24°C (HS) were obtained after Illumina sequencing. From 92 million raw reads, 77 million clean reads were identified ([Fig 2A](#pone.0186433.g002){ref-type="fig"} and [S2 Table](#pone.0186433.s005){ref-type="supplementary-material"}). High-quality reads were assembled into 30,672 unigenes by *Trinity* \[[@pone.0186433.ref042]\], annotated and mapped to public databases ([S3 Table](#pone.0186433.s006){ref-type="supplementary-material"}). The high Pearson correlation coefficient suggested good repeatability within each group ([S1 Fig](#pone.0186433.s001){ref-type="supplementary-material"}). Then the unigene expression (RPKM) data were calculated between CT and HS with *DESeq* \[[@pone.0186433.ref043]\], which showed that 27,180genes were expressed in both CT and HS as well as 678 genes were specifically expressed in CT ([Fig 3A](#pone.0186433.g003){ref-type="fig"} and [S3 Table](#pone.0186433.s006){ref-type="supplementary-material"}). With heat stress, 2,769 genes are expressed ([Fig 3A](#pone.0186433.g003){ref-type="fig"} and [S3 Table](#pone.0186433.s006){ref-type="supplementary-material"}). The comparison of cumulative RPKM values of all genes in CT and HS groups indicates that HS gene expression is greater than in CT ([Fig 3B](#pone.0186433.g003){ref-type="fig"}). Finally, 324 genes with an absolute log2-ratio value ≥1 were up- (232) or downregulation (92) when HS and CT were compared, which were identified as DEG ([Fig 3C](#pone.0186433.g003){ref-type="fig"}, [S2 Fig](#pone.0186433.s002){ref-type="supplementary-material"} and [S4 Table](#pone.0186433.s007){ref-type="supplementary-material"}).

![Length distribution of contigs, transcripts and unigenes of mRNA and miRNA in Tibetan naked carp.\
(a) Distribution of assemblies (contigs, transcripts and unigenes). The left y-axis and solid lines are the distributions of number (log10-transformed) of assemblies in each 100-bp bin, while the right y-axis and dashed lines are the cumulative curves for each assembly. (b) Distributions of sequence number of clean and unique sequence within miRNA-seq. All the plots were carried out with ggplot2 in R version 3.2.1.](pone.0186433.g002){#pone.0186433.g002}

![Comparisons between biological replication data sets of heat stress group (HS) and control groups (CT).\
(a) Number of expressed genes within HS and CT; (b) RPKM range of specific expressed genes between HS and CT; (c) Significantly up- and down-regulated mRNAs between HS and CT comparisons; (d) Significantly up- and down regulated miRNAs between HS and CT comparisons.](pone.0186433.g003){#pone.0186433.g003}

miRNA sequencing generated more than 60 million raw reads to yield 47 million clean reads which were narrowed to 41,585,187 reads ranging from 15 to 30 nt ([Fig 2B](#pone.0186433.g002){ref-type="fig"} and [S4 Table](#pone.0186433.s007){ref-type="supplementary-material"}). The Pearson correlation analysis also confirmed that the good repeatability among samples within group ([S1 Fig](#pone.0186433.s001){ref-type="supplementary-material"}). Data from CT and HS replicates were represented by more than 10 million clean reads, which was sufficient for the subsequent quantitative analysis of miRNAs. CT and HS shared 1,114 expressed miRNAs, but only 2 miRNAs were upregulated and 10 miRNAs were downregulated when comparing HS and CT groups ([Fig 3D](#pone.0186433.g003){ref-type="fig"}, [S2 Fig](#pone.0186433.s002){ref-type="supplementary-material"} and [S5 Table](#pone.0186433.s008){ref-type="supplementary-material"}) and these were identified as DEMs.

DEGs associated with normal skeletal development {#sec013}
------------------------------------------------

Collagen stabilization and regeneration could contribute to maintain the normal skeletal development in the early development of larval fish \[[@pone.0186433.ref044],[@pone.0186433.ref045]\]. Morphological evidences indicated that severe heat exposure (30°C) had induced malformation of fish larvae as the abnormal skeletal development, and also triggered high death rate in CT group ([Fig 1B and 1C](#pone.0186433.g001){ref-type="fig"}). Intriguingly, we found that two DEGs *collagen alpha-4 (IV) chain* (COL4A4) and *collagen alpha-1(XIV) chain* (COL14A1) that were involved into collagen stabilization and regeneration were significantly up-regulated in response to heat stress ([Table 1](#pone.0186433.t001){ref-type="table"} and [Fig 4A](#pone.0186433.g004){ref-type="fig"}). The malformation observed under heat stress was presumably associated with the altered expression of *COL4A4* and *COL14A1*, whose function was to maintain the normal skeletal development.

![Statistics comparisons of KEGG pathway and validation of RNA-seq data by RT-qPCR.\
(a) Comparisons of changes in miRNAs targets involved in the KEGG pathway. (b) Fold changes of gene expression detected by RNA-seq were plotted against the data of RT-qPCR. The reference line indicates the linear relationship between the results of RNA-seq and RT-qPCR. The correlation between RNA-seq and RT-qPCR data was analyzed by Spearman's rho test (r^2^ = 0.9752). (c) Expression comparisons of selected DEMs target genes according to RNA-seq and RT-qPCR. Bars of RNA-seq and RT-qPCR date were colored by red and purple, respectively. Data of RT-qPCR were shown as mean ± standard error of the mean (n = 3). (d) Comparisons of significant changes in miRNAs between HS and CT groups by RNA-seq and qPCR.](pone.0186433.g004){#pone.0186433.g004}

10.1371/journal.pone.0186433.t001

###### DEGs of pathway within Tibetan naked carp exposure to rapid heat stress.

![](pone.0186433.t001){#pone.0186433.t001g}

  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  Protein                                          Sub-catalog                                                  Log~2~(HS/CT)   Description
  ------------------------------------------------ ------------------------------------------------------------ --------------- ------------------------------------------------------------------
  *Collagen stabilization and regeneration*                                                                                     

  COL4A4                                           collagen catabolic process;\                                 2.487           collagen alpha-4 (IV) chain
                                                   extracellular matrix organization                                            

  COL14A1                                          collagen binding;\                                           2.296           collagen alpha-1(XIV) chain
                                                   extracellular matrix structural constituent                                  

  *Lipid metabolism*                                                                                                            

  NSDHL                                            Steroid biosynthesis                                         -1.079\*        sterol-4alpha-carboxylate 3-dehydrogenase

  CYP27B                                           Steroid biosynthesis                                         -1.128\*        25-hydroxyvitamin D3 1alpha-hydroxylase

  PLB1                                             Glycerophospholipid metabolis; Arachidonic acid metabolism   1.111           phospholipase B1

  SPLA2                                            Glycerophospholipid metabolis                                -1.043\*        secretory phospholipase A2

  ENPP7                                            Sphingolipid metabolism                                      -1.042\*        ectonucleotide pyrophosphatase/phosphodiesterase family member 7

  *Nucleotide metabolism; Amino acid metabolism*                                                                                

  PDE5                                             Purine metabolism                                            1.133           cGMP-specific 3\',5\'-cyclic phosphodiesterase

  MLL3                                             Lysine degradation                                           1.009           histone-lysine N-methyltransferase MLL3

  SUV420H                                          Lysine degradation                                           1.011           histone-lysine N-methyltransferase SUV420H

  NOS2                                             Arginine and proline metabolis                               1.364           nitric-oxide synthase 2

  AANAT                                            Tryptophan metabolism                                        -1.313\*        arylalkylamine N-acetyltransferase

  *Signal transduction*                                                                                                         

  MAP4K3                                           MAPK signaling pathway                                       1.088           mitogen-activated protein kinase kinase kinase kinase 3

  CAMK2                                            Wnt signaling pathway                                        1.104           calcium/calmodulin-dependent protein kinase

  NODAL                                            TGF-beta signaling pathway                                   2.067           nodal

  *Digestive system*                                                                                                            

  ADRB1                                            Salivary secretion                                           1.089\*         adrenergic receptor beta-1

  CELA2                                            Pancreatic secretion                                         1.391\*         pancreatic elastase 2

  APOA1                                            Fat digestion and absorption                                 -1.092\*        apolipoprotein A1

  APOB                                             Fat digestion and absorption                                 -1.034\*        apolipoprotein B

  APOA4                                            Fat digestion and absorption                                 -1.424\*        apolipoprotein A4

  CUBN                                             Vitamin digestion and absorption                             1.777           cubilin
  ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Lipid metabolism and digestive system pathway genes involved in rapid heat stress response {#sec014}
------------------------------------------------------------------------------------------

Lipid metabolism and digestive system pathways are key to fish physiology and maintain health \[[@pone.0186433.ref046]\]. DEGs involved in signal transduction, digestive system and lipid metabolism pathways were significantly higher than others ([Fig 4A](#pone.0186433.g004){ref-type="fig"}). In the tricarboxylic acid (TCA) cycle, the *acylphosphatase* gene involved in pyruvate metabolism, and *phosphatidylinositol-bisphosphatase* gene involved in inositol phosphate metabolism were upregulated after heat stress, but no other gene was modified ([Table 1](#pone.0186433.t001){ref-type="table"} and [S6 Table](#pone.0186433.s009){ref-type="supplementary-material"}), suggesting an overall upregulation of the main pathway. DEGs involved in lipid metabolism, such as *secretory phospholipase A2* and *phospholipase B1* genes engaged in glycerophospholipid metabolism, ether lipid metabolism and arachidonic acid metabolism, *ectonucleotide pyrophosphatase/phosphodiesterase family member 7* in sphingolipid metabolism were significantly upregulated. In addition, DEGs involved in the endocrine system, such as *calcium/calmodulin-dependent protein kinase* (CAMK2) that played an essential role in insulin secretion, GnRH signaling pathway, oxytocin signaling pathway and melanogenesis, were upregulated in HS ([Table 1](#pone.0186433.t001){ref-type="table"} and [S7 Table](#pone.0186433.s010){ref-type="supplementary-material"}). In the digestive system, *adrenergic receptor beta-1* (ADRB1) gene that was involved in salivary secretion and the *pancreatic elastase 2* (CELA2) gene that was involved in protein digestion and absorption were upregulated. Meanwhile, *apolipoprotein A1* (APOA1), *B* (APOB) and *A4* (APOA4) genes participated in vitamin digestion and absorption were downregulated ([Table 1](#pone.0186433.t001){ref-type="table"} and [S8 Table](#pone.0186433.s011){ref-type="supplementary-material"}).

Regulation network of miRNA and target mRNA {#sec015}
-------------------------------------------

Using *in silico* analysis of miRNA target mRNAs with miRanda 4.0 and TargetScan, we obtained a total of 34,258 targeted mRNA and 4,167,259 sites ([S9 Table](#pone.0186433.s012){ref-type="supplementary-material"}). Mapping the DEM and their targeted mRNA yielded more than 3,000 targeted mRNA. Based on the GO and KEGG annotation ([S6](#pone.0186433.s009){ref-type="supplementary-material"} and [S8](#pone.0186433.s011){ref-type="supplementary-material"} Tables), many DEM targeted mRNAs were enriched in the digestive system, signal transduction and lipid metabolism pathways, in a good agreement with previous DEG annotation results ([S3 Fig](#pone.0186433.s003){ref-type="supplementary-material"}). Finally we constructed the network of DEM and their target mRNAs using *Cytoscape* 3.3.0 ([Fig 5](#pone.0186433.g005){ref-type="fig"}).

![Regulation network of DEMs and targeted DEGs in response to rapid heat stress within Tibetan naked carp.\
The red circle represents the DEM, and yellow circle represents the DEG, the line represents that there are links between the DEM and target DEG. All the DEM and target DEG were collected and predicted with *Cytoscape* 3.3.0 software. Image was carried out with *Coreldraw* 14.0 software.](pone.0186433.g005){#pone.0186433.g005}

Validation of RNA-seq data by qPCR {#sec016}
----------------------------------

After comprehensive RT-qPCR validation, expression of 21 selected DEGs and 12 DEMs were measured and compared with the expression profiles from RNA-seq analysis. The expression data for these selected DEGs and DEMs detected by RNA-seq and RT-qPCR were displayed in [S9 Table](#pone.0186433.s012){ref-type="supplementary-material"} and [Fig 4](#pone.0186433.g004){ref-type="fig"}. The data from RT-qPCR and RNA-seq exhibited an excellent agreement on both up- and down-regulated genes ([Fig 4B](#pone.0186433.g004){ref-type="fig"}). The correlation between RNA-seq and RT-qPCR data was analyzed by Pearson's correlation coefficient ([Fig 4B](#pone.0186433.g004){ref-type="fig"}), and a highly statistical significance ([Fig 4C and 4D](#pone.0186433.g004){ref-type="fig"}) was observed.

Discussion {#sec017}
==========

This study represented the first transcriptome-wide analysis of thermal acclimation in Schizothoracine fish species. Numbers of studies had revealed the thermal acclimation in fish species on a large-scale data level \[[@pone.0186433.ref047]--[@pone.0186433.ref050]\], while previous studies of thermal (cold) acclimation in Schizothoracine fishes still based on few candidate genes \[[@pone.0186433.ref051],[@pone.0186433.ref052]\]. Thanks to the rapid development and decreasing costs of RNA sequencing technologies in recent years, we were able to obtain more than 30,000 genes and 1,000 miRNAs in *G*. *przewalskii*. In addition, the combination of mRNA and miRNA transcriptome sequencing enabled us to obtain a larger number of candidates for thermal acclimation. This study also provided valuable genetic resources for future studies on environmental stress response in Schizothoracine fish.

Most TP lakes are frozen for long periods of time with low water temperature for aquatic organisms \[[@pone.0186433.ref053]\], and Schizothoracine fish species endemic to highland lakes could survive under this condition \[[@pone.0186433.ref013],[@pone.0186433.ref014]\]. *G*. *przewalskii* is a commercially important aquaculture species in Northwestern China, and it inhabits in Lake Qinghai with lengthy period of cold water (September to April) \[[@pone.0186433.ref014],[@pone.0186433.ref023],[@pone.0186433.ref053]\]. Water temperature is the chief environmental determinant of the growth performance of fishes \[[@pone.0186433.ref002]--[@pone.0186433.ref004],[@pone.0186433.ref046]\], while it is also the physiological constraint for the ongoing aquaculture industry of this coldwater fish species. In the current study, we focused on the controversial issue of the survival and the growth of *G*. *przewalskii* in high temperature condition \[[@pone.0186433.ref014],[@pone.0186433.ref023]\]. Evidences of earlier studies had revealed that thermal experiences of fish at embryonic stages could have dramatic and persistent effects on the postnatal thermal acclimation capacity \[[@pone.0186433.ref042],[@pone.0186433.ref046]\], and a number of studies had focused on thermal (cold) acclimation in fish larvae \[[@pone.0186433.ref005],[@pone.0186433.ref029],[@pone.0186433.ref048]\]. In this study, we designed a comparative experiment using *G*. *przewalskii* larvae as materials to investigate the potential thermal acclimation in Tibetan Schizothoracine fishes. Data showed that exposure of *G*. *przewalskii* larvae to 24°C for 12 h increased larval survival under further severe heat stress at 30°C, but the direct exposure to 30°C without a mild high temperature transition caused massive malformation and death. The result was in line with previous studies proposing that *G*. *przewalskii* possessed a rapid heat hardening (RHH) ability to establish a heat tolerance under mild heat stress \[[@pone.0186433.ref005],[@pone.0186433.ref029]\]. The current study represented initial evidence of RHH-like responses in *G*. *przewalskii* as well as lays a foundation for studying molecular mechanisms underlying the establishment of RHH in highland larvae in the TP.

How organisms respond to the environment is of increasing importance to biologists and ecologists with the concerning of global climate change \[[@pone.0186433.ref006],[@pone.0186433.ref054]\]. Thus, transcriptome studies to address how environmental stimuli affect gene expression can be used to identify different expressed genes (DEG) between HS and CT by mapping unigenes to a pre-assembly reference transcriptome dataset using a *DESeq* package \[[@pone.0186433.ref043]\]. Specifically, the rapid heat stress triggered and silenced gene expression in *G*. *przewalskii* \[[@pone.0186433.ref029],[@pone.0186433.ref055]\]. Interestingly, species-specific expressed genes were chiefly associated with lipid metabolism, indicating that *G*. *przewalskii* might increase their metabolic rate and lipid consumption upon exposure to heat stress. Also, HS and CT were not different with respect to RPKM values, suggesting that *G*. *przewalskii* acclimate to the high temperature and restore homeostasis eventually.

Using a rigorous set of thresholds, we confirmed that 232 genes were up-regulated and 92 genes were downregulated after rapid heat stress exposure. These DEGs are mainly enriched in energy metabolic and digestive system pathways, corresponding well with previous work in the olive flounder \[[@pone.0186433.ref055]\]. Temperature can affect fish metabolism \[[@pone.0186433.ref002],[@pone.0186433.ref004],[@pone.0186433.ref010]\], and lipids are the chief energy source for most organisms, therefore lipid metabolism is selected to meet energy demands during the larval stage \[[@pone.0186433.ref029],[@pone.0186433.ref047]\].

A striking finding of this study is the observation that the rapid heat stress had significant elevated gene expression of *COL4A4* and *COL14A1*, which functions in collagen stabilization and regeneration. Past evidences reveal that the stabilization and regeneration of collage contributes to maintaining normal skeletal development in larval fish early developmental stages \[[@pone.0186433.ref044],[@pone.0186433.ref045]\]. Notably, there was significant mortality and deformities in the control fish after a rapid heat shock (16 up to 30°C) by microscopic examination. This observation was in agreement with previous studies during the process of thermal acclimation of other fish species \[[@pone.0186433.ref056]--[@pone.0186433.ref058]\]. Additionally, *COL4A4* encodes one of the six subunits of type IV collagen, as the major structural component of basement membranes \[[@pone.0186433.ref059]\], involved into normal skeletal development \[[@pone.0186433.ref044]\]. Another evidence shows that an increase of *COL14A1* expression contributes to maintain the normal development of mouse axial skeleton \[[@pone.0186433.ref060]\]. Thus, the heat stress altered expression of these gene groups that could potentially be associated with the developmental issues observed. These findings also indicate a link between heat stresses alter collagen stabilization and regeneration related gene expression and skeletal deformity.

Only two miRNAs in the HS group were upregulated and 10 miRNAs were downregulated. Data show that miRNAs can regulate cell differentiation, organogenesis, development and growth \[[@pone.0186433.ref030],[@pone.0186433.ref061]\]. Two *let-7* family members, *let-7d-3p* and *let-7c-1-3p* were downregulated and they are involved in responses to rapid heat stress, which is similar to *let-7c* and *let-7d* expression pattern in response to heat stress in rat small intestine \[[@pone.0186433.ref062]\]. *Let-7* was subsequently identified as the first known *C elegans* miRNA \[[@pone.0186433.ref063]\]. *let-7* and its family members are highly conserved across species in sequence and function, such as cell proliferation and growth pathways \[[@pone.0186433.ref064]\]. In addition, *miR-145* and *miR-125* were uniquely expressed between HS and CT groups and they have significant roles in cell response to stresses such as heat shock and oxidative stress \[[@pone.0186433.ref065]\]. Furthermore, the function of the remaining eight DEMs are not clearly illustrated, therefore we speculate that these miRNAs are presumably involved in the rapid heat stress response regulation in Tibetan naked carp. Most DEGs described here are involved in metabolic pathways, and DEM target mRNAs are also related to metabolic pathways, suggesting that an important role in the heat stress response.

We predicted and constructed a network diagram of miRNAs and their targets mRNAs by *Cytoscape*. We initially identified more than 100 genes which were potentially targeted by 12 differentially expressed miRNAs according to *in silico* analysis, but most of target mRNAs were not significantly expressed according to HS and CT comparisons. Finally, 56 DEGs targeted by 12 DEMs were identified, including *hsp70*, *hsp90a* and *hsp90b* that were regulated by DEMs in this network, suggesting that miRNAs significantly regulated Tibetan naked carp responses to heat stress \[[@pone.0186433.ref066]\]. Myosins are a family of ATP-dependent motor proteins with roles in muscle contraction and other motility processes \[[@pone.0186433.ref067]\], and *myosin 8a*, *-11*and *-15* are in this network, suggesting a link with heat stress response regulation to muscle movement.

Supporting information {#sec018}
======================

###### Pearson correlation between samples (n = 2) in each groups (HS and CT).

(PDF)

###### 

Click here for additional data file.

###### Volcano plot of DEGs and DEMs of Tibetan naked carp response to rapid heat stress.

(PDF)

###### 

Click here for additional data file.

###### Gene ontology classification of the DEGs and DEMs target genes of Tibetan naked carp response to rapid heat stress.

(PDF)

###### 

Click here for additional data file.

###### Primer list of selected DEGs and DEMs in RT-qPCR verification experiments.

(XLS)

###### 

Click here for additional data file.

###### Overview of sequencing, assembly and analysis of Tibetan naked carp transcriptome.

(XLS)

###### 

Click here for additional data file.

###### Annotation of all Tibetan naked carp unigenes.

(XLS)

###### 

Click here for additional data file.

###### List of statistically significant different expression genes (DEGs).

(XLS)

###### 

Click here for additional data file.

###### The result of conservative miRNA annotation analysis.

(XLS)

###### 

Click here for additional data file.

###### The significantly up/downregulated miRNAs in Tibetan naked carp response to rapid heat stress.

(XLS)

###### 

Click here for additional data file.

###### Enriched KEGG pathways of DEGs in Tibetan naked carp response to rapid heat stress.

(XLS)

###### 

Click here for additional data file.

###### GO terms enriched in up-regulated and down-regulated genes.

(XLS)

###### 

Click here for additional data file.

###### The result of DEMs and targeted DEMs analysis.

(XLS)

###### 

Click here for additional data file.

We thank colleagues from the National Rescue Center of Qinghai Lake Naked Carp for providing the technical assistance in fieldwork. The authors are also grateful to Drs. Renyi Zhang, Yongtao Tang, Chenguang Feng and Kunyuan Wanghe from our research group for facilitating field surveys and laboratory assistance.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.
