Integrated mRNA and microRNA transcriptome analyses reveal regulation of thermal acclimation in Gymnocypris przewalskii: A case study in Tibetan Schizothoracine fish

Environmental acclimation is important episode in wildlife occupation of the high-altitude Tibetan Plateau (TP). Transcriptome-wide studies on thermal acclimation mechanism in fish species are rarely revealed in Tibetan Plateau fish at high altitude. Thus, we used mRNA and miRNA transcriptome sequencing to investigate regulation of thermal acclimation in larval Tibetan naked carp, Gymnocypris przewalskii. We first remodeled the regulation network of mRNA and miRNA in thermal acclimation, and then identified differential expression of miRNAs and target mRNAs enriched in metabolic and digestive pathways. Interestingly, we identified two candidate genes contributed to normal skeletal development. The altered expression of these gene groups could potentially be associated with the developmental issues of deformity and induced larval death. Our results have three important implications: first, these findings provide strong evidences to support our hypothesis that G. przewalskii possess ability to build heat-tolerance against the controversial issue. Second, this study shows that transcriptional and post-transcriptional regulations are extensively involved in thermal acclimation. Third, the integrated mRNA and microRNA transcriptome analyses provide a large number of valuable genetic resources for future studies on environmental stress response in G. przewalskii and as a case study in Tibetan Schizothoracine fish.


Introduction
Fish survival is affected by abiotic stresses such as extreme temperatures, salinity, hypoxia and chemical toxicity [1,2], but water temperature is the chief environmental determinant of development, growth, reproduction, behavior, metabolism and geographical distribution [2][3][4]. As 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 wildcaught 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). This developmental stage allows a reseasonable period of post-hatching and no need for feeding, which provided suitable materials for environmental stress experiment [5,29]. 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). No embryonic mortality was observed before the treatment. After heat exposure (216 hpf, Fig 1A), 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. 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.
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). 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
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 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) with default parameters. Contigs of each assembly were performed with TGICL to produce long and complete unigenes [34], 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 nonredundant (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 [35]. Finally pathway assignments were generated using the KEGG database, Kyoto Encyclopedia of Genes and Genomes [36].

Gene expression analysis
Gene expression values were calculated as reads aligned to gene per kilo base of exon per million mapped reads (RPKM) [37]. Statistical comparison between two libraries was conducted with DESeq package in R software [38]. FDR (false discovery rate) < 0.05 was used as the threshold of p-value in multiple tests to measure significant gene expression difference [39]. 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 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.
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
In silico analysis of DEM targets were predicted in DEG using miRanda 4.0 algorithm (microrna.sanger.ac.uk/targets/v4) and TargetScan (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
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). 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) [40,41]. 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).

Effect of rapid heat stress on Tibetan naked carp larval development
After larvae exposure to severe high temperature at 30˚C (Fig 1A), 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). 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). 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
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 and S2 Table). High-quality reads were assembled into 30,672 unigenes by Trinity [42], annotated and mapped to public databases (S3 Table). The high Pearson correlation coefficient suggested good repeatability within each group (S1 Fig). Then the unigene expression (RPKM) data were calculated between CT and HS with DESeq [43], which showed that 27,180genes were expressed in both CT and HS as well as 678 genes were specifically expressed in CT (Fig 3A and S3 Table). With heat stress, 2,769 genes are expressed (Fig 3A and S3 Table). 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). 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,  S2 Fig and S4 Table).
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 and S4  Table). The Pearson correlation analysis also confirmed that the good repeatability among samples within group (S1 Fig). 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, S2 Fig and  S5 Table) and these were identified as DEMs.

DEGs associated with normal skeletal development
Collagen stabilization and regeneration could contribute to maintain the normal skeletal development in the early development of larval fish [44,45]. 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). 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 and Fig 4A). 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.
Lipid metabolism and digestive system pathway genes involved in rapid heat stress response Lipid metabolism and digestive system pathways are key to fish physiology and maintain health [46]. DEGs involved in signal transduction, digestive system and lipid metabolism pathways were significantly higher than others (Fig 4A). In the tricarboxylic acid (TCA) cycle, the acylphosphatase gene involved in pyruvate metabolism, and phosphatidylinosit ol-bisphosphatase gene involved in inositol phosphate metabolism were upregulated after heat stress, but no other gene was modified (Table 1 and S6 Table), 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 and S7 Table). 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 and S8 Table).

Regulation network of miRNA and target mRNA
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). Mapping the DEM and their targeted mRNA yielded more than 3,000 targeted mRNA. Based on the GO and KEGG annotation (S6 and S8 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). Finally we constructed the network of DEM and their target mRNAs using Cytoscape 3.3.0 (Fig 5).

Validation of RNA-seq data by qPCR
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 and Fig 4. The data from RT-qPCR and RNA-seq exhibited an excellent agreement on both up-and down-regulated genes (Fig 4B). The correlation between RNA-seq and RT-qPCR data was analyzed by Pearson's correlation coefficient (Fig 4B), and a highly statistical significance (Fig 4C and 4D) was observed.

Discussion
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 [47][48][49][50], while previous studies of thermal (cold) acclimation in Schizothoracine fishes still based on few candidate genes [51,52]. 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 [53], and Schizothoracine fish species endemic to highland lakes could survive Integrated mRNA and microRNA transcriptome analyses of thermal acclimation in Gymnocypris przewalskii under this condition [13,14]. 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) [14,23,53]. Water temperature is the chief environmental determinant of the growth performance of fishes [2][3][4]46], 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 [14,23]. 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 [42,46], and a number of studies had focused on thermal (cold) acclimation in fish larvae [5,29,48]. 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 [5,29]. 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 [6,54]. 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 [43]. Specifically, the rapid heat stress triggered and silenced gene expression in G. przewalskii [29,55]. 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 [55]. Temperature can affect fish metabolism [2,4,10], and lipids are the chief energy source for most organisms, therefore lipid metabolism is selected to meet energy demands during the larval stage [29,47].
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 [44,45]. 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 [56][57][58]. Additionally, COL4A4 encodes one of the six subunits of type IV collagen, as the major structural component of basement membranes [59], involved into normal skeletal development [44]. Another evidence shows that an increase of COL14A1 expression contributes to maintain the normal development of mouse axial skeleton [60]. 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 [30,61]. 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 [62]. Let-7 was subsequently identified as the first known C elegans miRNA [63]. let-7 and its family members are highly conserved across species in sequence and function, such as cell proliferation and growth pathways [64]. 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 [65]. 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 [66]. Myosins are a family of ATP-dependent motor proteins with roles in muscle contraction and other motility processes [67], and myosin 8a, -11and -15 are in this network, suggesting a link with heat stress response regulation to muscle movement.