Genome-Wide Characterization of Alternative Splicing Events and Their Responses to Cold Stress in Tilapia

Alternative splicing (AS) is an important post-transcriptional regulatory mechanism for cells to generate transcript variability and proteome diversity. No systematic investigation of AS events among different tissues in response to stressors is available for tilapia currently. In this study, AS among different tissues was identified and the cold stress-related AS events were explored in a Nile tilapia (Oreochromis niloticus) line based on 42 RNA-seq datasets using a bioinformatics pipeline. 14,796 (82.76%; SD = 2,840) of the expression genes showed AS events. The two most abundant AS types were alternative transcription start site (TSS) and terminal site (TTS) in tilapia. Testis, brain and kidney possess the most abundant AS gene number, while the blood, muscle and liver possess the least number in each tissue. Furthermore, 208 differentially alternative splicing (DAS) genes in heart and 483 DAS in brain in response to cold stress. The number of AS types for alternative exon end, exon skipping and retention of single intron increased significantly under cold stress. GO enrichment and pathway overrepresentation analysis indicated that many DAS genes, e.g., genes in circadian clock pathway, may influence expression of downstream genes under cold stress. Our study revealed that AS exists extensively in tilapia and plays an important role in cold adaption.


INTRODUCTION
Alternative splicing (AS) generates multiple functional mRNAs and therefore enables different proteins to be synthesized from one single gene (Maniatis and Tasic, 2002;Nilsen and Graveley, 2010). This regulatory mechanism significantly expands transcript variability and proteome diversity, as well as provides additional information in eukaryotic cells (Nilsen and Graveley, 2010). AS events widely exist in eukaryotes (Sorek, 2007). For example, 92-94% of human genes undergo regulation of AS (Wang et al., 2008). Many AS events are well conserved and appear to be related to the phylogenetic distance among species (Mastrangelo et al., 2012;Jonsson et al., 2016).
The flexibility of AS contributes to environmental adaption and phenotypic plasticity in organism (Mastrangelo et al., 2012). When environmental fluctuations are ubiquitous, mRNA can be processed into new isoforms with shifting expression level, so that cells can alter their metabolism and maintain homeostasis (Kucherenko and Shcherbata, 2018). For example, in duck, an AS event in mitf gene causes down-regulation of melanogenesis pathway, which led to plumage color change (Zhou Z. et al., 2018). The different isoforms from peptidoglycan recognition proteins (PGRPs) are required in recognition of different bacteria (Chang and Zhang, 2017). Stress-related AS events have been identified in many species, such as mild thermal, oxidative and heavy metal stress in yeast (Melangathli et al., 2017), temperature, salinity and air exposure stress in Pacific oyster (Huang et al., 2016;Liu and Guo, 2017), chemical stress or thermal stress in Drosophila (Fujikake et al., 2005;Kim et al., 2014;Ruan et al., 2015;Jaksic and Schlotterer, 2016), IFN antiviral response in zebrafish , heat shock in pigs (Ju et al., 2014) and stress treatment in human cells (Jacob et al., 2014), Edwardsiella ictaluri infection in catfish (Tan et al., 2018) and resistance to iridovirus in Asian seabass (Wang et al., 2017).
The high resolution of RNA-seq data allows for exploring not only gene expression but also genome-wide AS events in animals. RNA-seq has been extensively used for measuring gene expression to gain biological insights in aquatic animals recently (Chatchaiphan et al., 2017;Im et al., 2017;Gan et al., 2018;Shi et al., 2018;Wang et al., 2018;Yang et al., 2018;Ye et al., 2018;Yue et al., 2018). In teleost, some studies highlighted the important role of AS in some genes (Marchand et al., 2001;Chan and Cheng, 2004;Moore et al., 2005). AS frequencies ranging from 17% of mapped genes in zebrafish to ∼43% of mapped genes in the pufferfish were estimated using the public EST/cDNA datasets (Lu et al., 2010). However, only few studies paid attention to the AS responses to stress in fish (Zwollo, 2011;Fillatreau et al., 2013;Long et al., 2013;Wan and Su, 2015;Chang and Zhang, 2017;Tan et al., 2019). Further deep mining of transcriptomic data using bioinformatic tools to identify AS events will help clarify gene regulation mechanism under stressors in fish.
Cold is a severe stressor to fish in aquatic environment. Most of cold-related AS studies are focused on plants (James et al., 2012;Pose et al., 2013), insects (Majercak et al., 1999) and mammals (Gemignani et al., 2002;Sano et al., 2015). In zebrafish, differential expression analysis found that RNA splicing is one of the key biological processes; and spliceosome is a highly enriched pathway for genes that up-regulated under cold stress (Long et al., 2013). Moreover, AS and promoter switching of genes are found to be regulated under cold stress (Long et al., 2013). Several studies have paid attention to characterize the transcriptional responses to cold stress in fish brain or heart. For example, differential gene expression was explored in the brain of the channel catfish (Ictalurus punctatus) (Ju et al., 2002). Previous studies also identified temperature acclimation in the heart transcriptome of the rainbow trout (Oncorhynchus mykiss) (Vornanen et al., 2005). There studies indicate the importance of these tissues in the cold stress response in fish.
Although the importance of AS mechanisms in animal stress response has been extensively reported, surprisingly few studies  have surveyed AS events in fish at genome-wide level. Here, we characterized genome-wide AS events among ten different tissues and then identified cold stress-related AS events in brain and heart of a Nile tilapia line based on 42 RNA-seq datasets using a bioinformatics pipeline. Our findings suggested widely existence of AS events in tilapia and their important roles in cold adaptation.

Fish Management, Cold Treatment, and Sample Collection
Thirty-six Nile tilapia individuals (GIFT line; ∼90 dph) were randomly selected from a population raised at the fish facility of School of Life Sciences, Sun Yat-sen University. These fishes were evenly divided into two indoor conical fiber glass tanks (water depth: 70 cm, volume: 500 L) in a recirculating freshwater system for acclimation of 2 days with water temperature at 25-28 • C and dissolved oxygen (DO) > 6 mg L −1 . The fishes were fed with tilapia pellets twice a day and the photoperiod was adjusted to 12 D:12 L in the room.
During low temperature treatment, the fishes in one tank were used as the control and those in the other were considered as test samples. The temperature in the water of test tank was slowly decreased to 15 • C within 12 h using a water cooler. To investigate the AS changes under cold stress, six Nile tilapia in test tank were then cultured in 15 • C for 24 h, other six in control tank were maintained at 26 • C. The detailed information on stress treatment was stated in Li et al. . During treatment, no feed was given to fish. Tricaine mesylate (MS-222) is used for euthanasia of fish. A dose of ∼ 400 ppm for immersion of the fish for ∼10 min is usually applied. Brain and heart tissues were collected from each fish of each tank. The tissue samples were frozen immediately in liquid nitrogen and then stored at −80 • C until RNA isolation.

RNA Isolation and NGS Sequencing
Total RNA from cold-treated samples was isolated using TRIzol reagent (Invitrogen, United Kingdom) according to the manufacturer's protocol. The RNA quality was assessed using the Nanodrop-2000 (Thermo Fisher Scientific, United States) and electrophoresis in 1.5% agarose gel. Total RNA integrity was further evaluated using Bioanalyzer 2100 (Agilent Technologies). Every two samples from the same tissue of the same group were pooled equally. Twelve libraries from heart and brain samples that constructed by TruSeq TM RNA sample preparation kit according to the product instruction (Illumina) were finally sequenced using Illumina HiSeq2500 for pair-end (PE) 150 sequencing.

Identification of AS Events From RNA-seq Datasets
To identify AS events in different tissues of Nile tilapia, we downloaded 30 Nile tilapia RNA-seq data from the NCBI database 1 . The 30 RNA-seq datasets (Supplementary Additional File 1) were generated from 10 tissues (testis, brain, kidney, eye, 1 https://www.ncbi.nlm.nih.gov/ ovary, skin, heart, blood, muscle, and liver). For each tissue, three replicates were available for AS identification.
Both RNA-seq datasets including the 12 samples that generated in our experiment and the 30 samples downloaded from NCBI database were aligned to the tilapia reference genome (Orenile1.0; downloaded from the Ensembl database) 2 using the software TopHat (Trapnell et al., 2009). The read alignments of RNA-seq data were conducted using SAMtools (Li et al., 2009). Aligned reads were then assembled using Cufflinks (Trapnell et al., 2012) and compared to the tilapia reference annotation using Cuffcompare (Trapnell et al., 2010). The "extract-as" program of ASprofile software package was applied to determine AS events among multiple samples (Florea et al., 2013) using the gtf file that generated by cufflinks and the tmap file that generated by cuffcompare. The result of ASprofile was sorted by 12 AS types including SKIP, IR, AE, transcription start site (TSS), alternative transcription terminal site (TTS), MSKIP, XSKIP, XMSKIP, XIR, MIR, XMIR, and alternative exon ends (XAE). Significance test of AS types among 10 tissues was performed using t statistics in Microsoft Office Excel.

AS Events in Response to Cold Stress
TopHat software (Trapnell et al., 2009) can identify novel transcripts from locations of regions and estimate abundance of the transcripts from their depth of coverage in the mapping. Bam files that generated by TopHat were used to identify differential AS (DAS) in response to cold stress. We processed the aligned RNA-seq data using QoRTs package (Hartley and Mullikin, 2015). Once alignment and quality control were completed for the dataset, the coverage counts (for genes, exons, and known/novel splice-junctions) were generated via QoRTs (Hartley and Mullikin, 2015). DAS events between groups were then analyzed using JunctionSeq (Hartley and Mullikin, 2016). False discovery rate (FDR) value of 0.01 was used as the significant threshold for each comparison. Gene structure was displayed using online tool Gene Structure Display Server 2.0 (Hu et al., 2015) and using the software FasParser (Sun, 2017), and protein structure of AS genes was plotted with online tool SMART (Letunic et al., 2015).

Different Expression Gene Analysis
Different expression gene analysis of brain and heart samples were conducted using the cuffdiff program (Trapnell et al., 2013). The generated DEG data were then filtered using the parameters: | log2(fold_change)| > 1 and the summary count values for each gene in two datasets >20 and FDR ≤ 0.05.

GO Annotation and Pathway Enrichment
The Ensembl gene IDs for both DAS and DEG genesets were transferred to Gene symbols in David website 3 . We carried out overrepresentation test (Released 20181113) for both DAS and DEG datasets by using the online tool Panther database 4 (PANTHER version 14.0 released on 2018-12-03) and the zebrafish geneset as the reference list. Several analysis types were applied for the datasets including PANTHER GO-Slim biological process, molecular function, cellular component, PANTHER protein class and Reactome pathways. The Fisher's exact test was performed with FDR correction and a FDR of <0.05.
We also investigated the enrichment of the affiliated genes in the GO Molecular Functions, GO Biological Processes and GO Cellular Components datasets, using the online software Metascape 5 (Tripathi et al., 2015). The genes were annotated using the default resources that provided by Metascape. The p value cutoff was set to 0.01.

RT-PCR Analysis
RT-PCR analysis was performed to verify AS events. Primer pairs (Supplementary Additional File 2) for 10 selected exon-skipped events that predicted by Junctionseq were designed using NCBI primer designing tool 6 . The primers were located in the upper and lower exons of the targeted one. RNA was treated with DNase I to avoid the contamination effect of genomic DNA. Reverse transcription was carried out using the reverse transcription kit (TAKARA PrimeScript TM RT reagent Kit with gDNA Eraser) according to the manufacturer's instructions. PCR amplifications were conducted in a 20 µL volume using a 2 × PCR mastermix (Dongsheng, China) and 10 µg cDNA, 1 µL of forward and reverse primers (10 µM) in a Bio-Rad cycler. The following program was applied including one cycle of 3 min at 94 • C, 35 cycles of 30 s at 94 • C, 30 s at 60 • C and 30 s at 72 • C, followed by a final extension of 10 min at 72 • C. The PCR products were 5 http://metascape.org 6 www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi FIGURE 3 | The number of AS events in brain and heart of tilapia in response to cold stress. We compared the number of AS genes between the control group (26 • C) and the test group (16 • C). (A) * Showed the number of AS events increases significantly in brain (p < 0.01). (B) * Showed the AS events increase significantly in heart (p < 0.01).
detected by running 2.5% agarose electrophoresis. All reactions were conducted in triplicate.

Genome-Wide Identification of AS Events in Different Tissues of Nile Tilapia
To obtain an overview of AS events in tilapia genome, 30 RNA-seq datasets were downloaded from the NCBI database (Supplementary Additional File 1). These datasets consisted of RNA-seq reads from 10 tissues, including testis, brain, kidney, eye, ovary, skin, heart, blood, muscle, and liver. For each tissue, three replicates were available. The total base number of reads in our data sets was 158.8 Gb and the average mapping ratio was 78.9% (70.9-88.0%).
Only successfully mapped high quality reads to the reference genome were used for ab initio transcriptome assembly with cufflinks. From the assembly, cuffcompare and ASprofile were applied to predict AS events in different samples. The average number of expressed genes in 10 tissues was 17,878 (SD = 3,628). Surprisingly, 14,796 (SD = 2,840) of the genes (82.76%) showed AS events (Figure 1). Ovary (91.28%) and muscle (87.64%) had the highest rate of AS, while kidney (78.06%) and liver (79.66%) had the lowest rate, compared to the number of expressed genes in each tissue.
The AS events in transcriptome were classified into 12 types according to their structures (Florea et al., 2013). The numbers of AS events for 10 tissues were calculated based on the 30 NCBI RNA-seq datasets (Figure 2). Testis had the most abundant AS events (51,285), while the number in liver was the lowest (23,246).
TSS and TTS were the most abundant AS types in all tissues (Figure 2). These events were likely to be located in untranslated region (UTR), and might have smaller influence to final protein sequences than other types, e.g., retention of single intron (IR), exon skipping (SKIP), and alternative exon end (AE). However, the 5 UTR and 3 UTR were also essential in biology process, involved into the modulation of messenger RNA (mRNA) transcription, secondary structure, stability, localization, translation, and accessed to regulators like microRNAs (miRNAs) and RNA-binding proteins (RBPs) (Leppek et al., 2018;Steri et al., 2018). The tissue-specific AS events were rare. One tissue specific DAS gene was found for heart, kidney, liver, and ovary respectively; and seven specific AS genes were identified in skin (Supplementary Additional File 3). Most of gene splicing patterns were shared by at least two tissues. However, much more work needs to be done to validate the results of bioinformatics analysis.

AS Events Identified in Tilapia Under Cold Stress
To investigate the AS regulation in response to cold, 12 RNA-seq data was generated for a Nile tilapia line under cold stress. These data were available at the DDBJ database (BioProject Accession: PRJDB6721). In brain, a total of 483 DAS genes were found in response to cold stress, including 426 exons and 527 junctions with significant differentially expressed counts (p < 0.01). In heart, the DAS number was 208, with 166 differentially expressed exons and 178 differentially expressed junctions (p < 0.01) (Supplementary Additional File 4).
In brain, the number of AE, SKIP, IR, continuous exons skipping (MSKIP), alternative exon ends (XAE), multiple discontinuous exons skipping (XSKIP), and corresponding introns retention events (MIR, XMIR) of AS events significantly increased in response to cold stress (p < 0.01) (Figure 3A). In heart, the SKIP, XSKIP and XMSKIP events increased significantly in response to cold stress (p < 0.01) (Figure 3B). These events were likely to cause amino acid changes in proteins.
GO enrichment analysis suggested that many GO categories e.g., mRNA processing (GO:0006397) were significantly enriched in brain. The plasma membrane bounded cell projection organization (GO:0120036) in the biological process category had the most abundant DAS genes. In heart, DAS genes were significantly enriched in multiple categories including protein-containing complex assembly (GO:0065003), muscle cell development (GO:0055001) and heart process (GO:0003015). Protein-containing complex assembly (GO:0065003) had the largest number of DAS genes (20, Log 10 (P) = −4.47). Gene

RT-PCR Validation of DAS Genes
To validate the accuracy and reliability of analysis process, RT-PCR was applied to measure expression level of different transcripts from one gene. For exon skipping or intron retention events, the primer pairs were designed in the upper exon and lower exon of the events, so that nucleotide fragments with different length could be amplified by PCR if there exists AS (Supplementary Additional File 2).
For example, in gene ENSONIG00000015870, an exon skipping event occurred in response to cold stress. The exonincluded PCR product was 378 bp, while the exon-excluded DNA fragment was 180 bp. As seen from Figure 5, the exon-including transcripts were more abundant in cold-treated group (case) than the control group (control). The result was consistent with bioinformatics analysis.
The RT-PCR result of 8 DAS genes in brain and heart revealed that most of AS events (7/8) that predicted in our experiment was credible, except for the gene ENSONIG00000002990 in which an unexpected 450 bp sequence was amplified, but a 211 bp fragment was supposed (Figure 6). This may be originated from individual specific AS events.

Links Between DAS and DEG in Response to Cold Stress
It was observed that many gene overlaps between DAS and DEG datasets. For example, 148 genes were shared between brain DAS and DEG datasets and 64 genes shared between heart DAS and DEG datasets. A total FIGURE 6 | RT-PCR validation of AS events that predicted in RNA-seq data. RT-PCR validation was performed for AS events in eight genes that predicted by bioinformatic analysis. Except for the gene ENSONIG00000002990 producing a ∼450 bp band which was longer than the predicted size 211 bp, the results for other seven genes were as expected.
Frontiers in Genetics | www.frontiersin.org FIGURE 7 | GO enrichment of DEGs that identified in tilapia heart and brain. (A) GO enrichment of DEGs in brain. (B) GO enrichment of DEGs in heart. The horizontal axis showed the gene number in every pathway; and value of log 10 (P) was given after the bar.
There were many DAS and DEG genes, which could be classified into the same pathway categories (Supplementary  Additional File 9). Interestingly, there was a pathway, circadian regulation of gene expression (GO:0032922), with enrichment of both DAS and DEG genes in heart and brain in response to cold stress. Three DAS genes (per3, per2, per1a) and eight DEG genes in brain (clocka, per3, per2, bhlhe40, nfil3-5, nr1d2a, nr1d1, per1a) and heart (clocka, per3, per2, bhlhe40, nfil3-5, nfil3-6, nr1d1, nampta) were classified into this category. The gene list was shown in Supplementary Additional File 8. Among these, per1 and per2 genes showed exon skipping both in brain and heart under cold stress. An 147bp exon in per1 was cut off (Figure 9A). For per2, a 208bp exon was spliced, which leaded to a deletion in FIGURE 8 | Venn diagram of DASs and DEGs in tilapia brain and heart. There are 208 DAS genes in heart, 483 DAS genes in brain, 2,918 DEGs in brain and 1,673 DEGs in heart. Figure 9B). Both AS events caused a deficiency of a part of amino acid sequences in proteins (Figures 9A,B).

Developing a Reference Dataset for Studying on AS Function in Tilapia
Alternative splicing was of great importance in improving regulatory capacities and proteomic complexities in eukaryotes (Ben-Dov et al., 2008). In teleost, some studies highlighted the important role of AS in a few genes, e.g., thyroid hormone receptors in teleost fish (Marchand et al., 2001;Chan and Cheng, 2004;Moore et al., 2005). However, few studies carried out genome-wide characterization of AS events in fish (Lu et al., 2010;Long et al., 2013). In this study, we explored the genome-wide AS events in 10 tissues of Nile tilapia based on 42 RNA-seq datasets. This represented the largest datasets used in the characterization of AS events in tilapia.
In Nile tilapia, the average percentage of AS genes was 82.76% with a range from 78.06% (kidney) to 91.28% (ovary) across 10 tissues. The value was comparable to that in mammal species. For example, 94.4-95.5% of expressed genes in ovary of pig (Tang et al., 2018) and about 70% of genes in human might be alternatively spliced (Johnson et al., 2003). Surprisingly, the value of AS frequency in tilapia was far more than that identified in reported fish species. The AS frequency from 17% of mapped genes in zebrafish to ∼43% of mapped genes in the pufferfish were identified based on a collective EST and CDS datasets from public database (Lu et al., 2010). The low frequencies may be mainly caused from the smaller datasets used in two studies. Further analysis based on larger datasets from different tissues may identify more AS expectedly.
The high abundance of AS numbers in tilapia testis and brain was similar to human (Yeo et al., 2004;Elliott and Grellscheld, 2006). Relatively fewer AS events were found in tilapia liver when compared to the findings in other tissues. The reason behind this might be the well-differentiation and low-complexity of liver tissue (Yeo et al., 2004). We also found 208 DAS genes in heart and 483 DAS genes in brain in response to cold stress. Under cold stress, the number of AE, SKIP and IR increased significantly. Our first investigation of genome-wide AS patterns across 10 tissues of Nile tilapia based on a large dataset provides a basis for study on AS function in tilapia.

mRNA Splicing Might Play an Important Part in Temperature Adaptation
Alternative splicing plays important roles in phenotype determination and biological response to environmental stress. The cold-related DEG analysis had been reported in liver and kidney (Yang et al., 2015;Zhou T. et al., 2018) in previous studies, and heart and brain of tilapia in this study. These studies indicated that cold stress might affect many pathways including metabolism and immunity in tilapia. However, cold-related RNA splicing was not deeply investigated in tilapia.
Splice site selection was under complex regulation by various cis-regulatory elements in pre-mRNAs and trans-regulatory elements. Many RBPs were shown to be critically required for the regulation of AS, such as serine/arginine-rich (SR) and heterogeneous nuclear RNP (hnrnp) protein families (Nilsen and Graveley, 2010;Wang et al., 2015). In the DAS gene list of catfish under heat stress, there were five SR genes (srsf2, srsf3, srsf5, srsf7like, and srsf11-like) and two heterogeneous hnrnp genes (hnrnpm and hnrnpa0) (Tan et al., 2019). In zebrafish, RNA splicing was listed as one of the most highly overrepresented biological processes, while spliceosome was one of the most highly enriched pathways for genes up-regulated under cold stress (Long et al., 2013). These findings suggested that RBPs act in AS-related stress regulation in fish.
In our cold treatment to tilapia, pathway enrichment analysis found overrepresentation of DAS and DEG genes in RNA splicing pathway. In spliceosome pathway, many SR splicing factors were regulated under cold stress. For example, srsf1, srsf3, srsf5 and srpk3 underwent AS regulation. Besides, srsf3, srsf5, srsf6, srpk2 genes showed differential expression level. In this pathway, ddx5 gene was spliced differentially under cold stress in both brain and heart. Its complementary, rbm8a, was up-regulated in brain and heart. Several downstream genes were also influenced, such as prpf and slu7 (Supplementary Additional File 2). Hnrnp protein families (hnrnpc, hnrnpd, hnrnph3, hnrnpll, and hnrnpu) also had differential expression level. Therefore, temperature stress affected the expression and splicing pattern of many genes in the spliceosome pathway of tilapia.

Circadian Clock Pathway Might Play Important Roles in Cold Adaptation
Alternative splicing in per gene family was ubiquitous in many species. In zebrafish, alternative promoter usage was detected for per3 gene under cold stress, which leads to a highly up-regulated transcript encoding a truncated protein lacking the C-terminal domains (Long et al., 2013). Different AS of per also existed in honeybees (Takeda et al., 2004). Per gene was a key gene in the circadian clock pathway, whose protein product combined with CRY protein to suppress the expression of arntl/clock genes (Zambon et al., 2003). ARNTL:CLOCK protein complex went on to regulate the expression of per gene and cry gene. These composed a closed cycle pathway for biological rhythm. In the DEG datasets, the expression level of arntl and clock were both up-regulated. The change in per protein sequences might influence the function of arntl/clock protein, therefore weaked the suppression to gene expression of arntl and clock.
The relationship between circadian clock pathway and temperature was proved in previous study (Kidd et al., 2015). When Drosophila adapted to seasonally cold days, a thermosensitive splicing event in 3 UTR played an important role (Majercak et al., 1999). In Drosophila, per generated two transcript types with an alternative intron in the 3 UTR which was a key regulation on how the Drosophila circadian clock adapted to changes in temperature (Cheng et al., 1998;Majercak et al., 2004). The timeless (tim) gene in Drosophila melanogaster had two different transcript variants that were thermo-sensitively spliced (Boothroyd et al., 2007;Montelli et al., 2015). In zebrafish, alternative promoter usage was detected for per3 gene after exposure to cold stress (16 • C) for 24 h (Long et al., 2013).
The heme biosynthesis was believed to provide a key connection as a signaling molecule to indicate the nutritional status to the circadian oscillatory machinery and allow it to respond appropriately (Burris, 2008). In our DEG dataset, three genes (alas1, cpox, fech) in this pathway had decreased expression levels. According to our result, the rate-limiting enzyme of heme biosynthesis, alas1, had differential AS pattern and downregulated expression level in heart. In organism, heme was a compound that played an important role in oxygen carrying. We speculated that the response of heme biosynthesis was related to oxygen consumption rate under cold.
Genes acting in heme were also regulated. Rev-erbα (also known as NR1D1) was a heme sensor that regulated the metabolic gene pathway, thus serving for coordination of circadian and metabolic pathways (Yin et al., 2010). REV-ERBα expression was regulated by the circadian genes binding to E-boxes within the Rev-erbα promoter. In our research, expression of nr1d1 was significantly up-regulated in brain. nr1d4b and nr1d2a were also up-regulated in both brain and heart. Neuronal PAS domain protein 2 (NPAS2) was a heme-binding protein and played a role in the regulation of circadian rhythms. The DNA binding activity of NPAS2 was related to heme (Dioum et al., 2002). Heme could also combine with NR1D protein to repress expression of arntl and npas2. In our result, the expression of npas2 was decreased in heart.
Rev-erbα response element was also responsive to PPARα (Canaple et al., 2006). PPARα played a pivotal role in the regulation of lipid homeostasis, fatty acid oxidation and acted as a key nutritional and environmental sensor for metabolic adaptation (Contreras et al., 2013;Seok and Cha, 2013). In our result, PPARα had differential AS pattern under cold stress in heart. The role of PPARα was a master regulator of lipid metabolism via regulation of numerous genes (Kersten, 2014). The down-regulation of most genes in cholesterol biosynthesis pathway, which indicated the close relationship between metabolic changes and fluctuations in biological rhythms.

CONCLUSION
Alternative splicing was an essential mechanism in stress response. In this study, we characterized the genome-wide AS events in ten different tissues and identified cold tolerant related AS in heart and brain of Nile tilapia based on 42 RNA-seq datasets. Our findings suggest the widely existence of AS and its important roles in cold adaptation in tilapia. Further functional studies on AS events will help clarifying gene regulation mechanism under stress in tilapia.

DATA AVAILABILITY STATEMENT
All the read data under cold stress were available at the DDBJ database (BioProject Accession: PRJDB6721). All other datasets have been presented within the article/Supplementary Material.

ETHICS STATEMENT
All experiments in this study were approved by the Animal Care and Use Committee of the School of Life Sciences at Sun Yat-sen University and were performed according to the regulations and guidelines established by this committee.

AUTHOR CONTRIBUTIONS
JX and HL contributed to the project conception. BL, ZZ, HQ, and ZM conducted the experiment and data analysis. JX and BL prepared the manuscript. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We thank the reviewers for comments on the manuscript.