Abstract
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 (; ). This regulatory mechanism significantly expands transcript variability and proteome diversity, as well as provides additional information in eukaryotic cells (). 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 (; ). AS participates in organogenesis (), development (Yamaguchi and Iwasa, 2018), tissue identity (; ), disease occurrence (; ) and heavy metals, cold and heat, pathogen attacks (; ; ; ) in eukaryotes.
The flexibility of AS contributes to environmental adaption and phenotypic plasticity in organism (). 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 (). 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 (). Stress-related AS events have been identified in many species, such as mild thermal, oxidative and heavy metal stress in yeast (), temperature, salinity and air exposure stress in Pacific oyster (; ), chemical stress or thermal stress in Drosophila (; ; ; ), IFN antiviral response in zebrafish (Zhang et al., 2018), heat shock in pigs () and stress treatment in human cells (), 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 (; ; ; ; Wang et al., 2018; Yang et al., 2018; Ye et al., 2018; Yu et al., 2018; Yue et al., 2018). In teleost, some studies highlighted the important role of AS in some genes (; ; ). 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 (). However, only few studies paid attention to the AS responses to stress in fish (Zwollo, 2011; ; ; Wan and Su, 2015; ; 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.
Tilapia is one of the most important agricultural fish worldwide (). Currently, Nile tilapia (Oreochromis niloticus), blue tilapia (Oreochromis aureus), Mozambique tilapia (Oreochromis mossambicus) and various hybrids including red tilapia are the major cultured species/strains worldwide (Xia et al., 2015; ). The gene expression changes in response to different stressors including cold stress (), alkalinity stress (Zhao et al., 2015), salinity adaptation (Xu et al., 2015), hypoxia stress (; Wang et al., 2017; Xia et al., 2018), cortisol stress (Tsalafouta et al., 2018), thermal stress (Ventura et al., 2018; Yue et al., 2018), oxidative stress (), and immune response () have been explored by RNA-seq in tilapia. Genome-wide identification of copy number variation (), QTL intervals for hypoxia and salt tolerance traits (; ) and long non-coding RNAs (Yue et al., 2018) have been performed. In tilapia, AS events were explored for hypoxia stress based on RNA-seq datasets (; Xia et al., 2018) and nutrition status based on gene cloning () recently. No systematic investigation of AS events among different tissues in response to stressors is available for tilapia.
Cold is a severe stressor to fish in aquatic environment. Most of cold-related AS studies are focused on plants (; ), insects () and mammals (; ). 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 (). Moreover, AS and promoter switching of genes are found to be regulated under cold stress (). 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) (). 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 (Xu et al., 2015) 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.
Materials and Methods
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. (Yue et al., 2018). 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 TruSeqTMRNA 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 database1. The 30 RNA-seq datasets (Supplementary Additional File 1) were generated from 10 tissues (testis, brain, kidney, eye, 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 (). 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 () 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 (). 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 (). DAS events between groups were then analyzed using JunctionSeq (). 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 () and using the software FasParser (Sun, 2017), and protein structure of AS genes was plotted with online tool SMART ().
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 website3. We carried out overrepresentation test (Released 20181113) for both DAS and DEG datasets by using the online tool Panther database4 (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 Metascape5 (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 tool6. 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 PrimeScriptTMRT 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 detected by running 2.5% agarose electrophoresis. All reactions were conducted in triplicate.
Results
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.
FIGURE 1
The AS events in transcriptome were classified into 12 types according to their structures (). 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).
FIGURE 2
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) (; 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.
FIGURE 3
Reactome pathway enrichment analysis found that DAS genes in mRNA Splicing – Major Pathway (FDR = 2.73E-02, fold enrichment = 4.54), mRNA splicing pathway (FDR = 3.07E-02, fold enrichment = 4.33) and Processing of Capped Intron-Containing Pre-mRNA (FDR = 5.37E-03, fold enrichment = 4.42) were overrepresented in brain. However, in heart, DAS genes were overrepresented in much more pathways, e.g., metabolism of RNA (FDR = 1.76E-04, fold enrichment = 4.54), translation (FDR = 2.16E-05, fold enrichment = 8.17), eukaryotic translation initiation (FDR = 1.76E-04, fold enrichment = 4.54), major pathway of rRNA processing in the nucleolus and cytosol (FDR = 1.26E-05, fold enrichment = 10.17) (Table 1).
TABLE 1
| Tissue value | Reactome pathways | Gene number | Fold enrichment | FDR |
| brain | mRNA Splicing – major pathway | 11 | 4.54 | 2.73E-02 |
| Processing of capped intron-containing pre-mRNA | 14 | 4.42 | 5.37E-03 | |
| mRNA splicing | 11 | 4.33 | 3.07E-02 | |
| heart | RHO GTPases activate PAKs | 3 | 24.41 | 3.49E-02 |
| SRP-dependent cotranslational protein targeting to membrane | 10 | 20.1 | 3.38E-07 | |
| Non-sense mediated decay (NMD) independent of the exon junction complex (EJC) | 10 | 19.2 | 2.56E-07 | |
| Formation of a pool of free 40S subunits | 10 | 17.62 | 3.70E-07 | |
| Non-sense mediated decay (NMD) enhanced by the exon junction complex (EJC) | 10 | 16.27 | 5.65E-07 | |
| Non-sense-mediated decay (NMD) | 10 | 16.27 | 4.52E-07 | |
| GTP hydrolysis and joining of the 60S ribosomal subunit | 10 | 15.12 | 7.29E-07 | |
| L13a-mediated translational silencing of ceruloplasmin expression | 10 | 15.12 | 6.25E-07 | |
| Cap-dependent translation initiation | 10 | 14.12 | 1.01E-06 | |
| Eukaryotic translation initiation | 10 | 14.12 | 8.98E-07 | |
| Major pathway of rRNA processing in the nucleolus and cytosol | 10 | 10.17 | 1.51E-05 | |
| rRNA processing | 10 | 10.17 | 1.37E-05 | |
| rRNA processing in the nucleus and cytosol | 10 | 10.17 | 1.26E-05 | |
| Translation | 11 | 8.17 | 2.16E-05 | |
| Metabolism of RNA | 15 | 4.54 | 1.76E-04 |
Reactome pathway overrepresentation analysis of DAS events in brain and heart of Nile tilapia under cold stress.
Pathways with FDR < 0.05 were showed in the table. Gene number, fold enrichment and FDR values were given.
Panther protein classification found that the DAS genes in brain were enriched in multiple categories including mRNA processing factor (FDR = 3.93E-04, fold enrichment = 5.95), RNA binding protein (FDR = 3.06E-02, fold enrichment = 2.4) and mRNA processing factor (FDR = 3.93E-04, fold enrichment = 5.95), while the DAS genes in heart were enriched in the category of ribosomal protein (FDR = 2.91E-02, fold enrichment = 7.02), RNA binding protein (FDR = 6.24E-05, fold enrichment = 4.55), and actin family cytoskeletal protein (FDR = 3.66E-02, fold enrichment = 4.48) (Supplementary Additional File 5).
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, Log10(P) = −4.47). Gene number and fold enrichment of every GO term were shown in the Figure 4.
FIGURE 4
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 exon-included 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.
FIGURE 5
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.
FIGURE 6
Different Expression Genes Under Cold Stress
Alternative splicing events affected expression of downstream genes (Zhou Z. et al., 2018). To investigate the effect to downstream genes, we further carried out a DEG analysis of our transcriptome datasets. There are 2,918 DEGs in brain and 1,673 DEGs in heart in response to cold stress (15°C) (Supplementary Additional File 6).
In brain, significant overrepresentation of DEGs was found in 13 Reactome pathways, e.g., metabolism (FDR = 1.62E-07, fold enrichment = 1.58), mRNA splicing (FDR = 1.21E-02, fold enrichment = 2.27), cell cycle (FDR = 2.82E-02, fold enrichment = 1.67), metabolism of lipids (FDR = 3.30E-04, fold enrichment = 1.81), and metabolism of protein (FDR = 9.47E-03, fold enrichment = 1.39) (Supplementary Additional File 7). There were 11 protein class categories were significantly overrepresented including RNA binding protein, transaminase, isomerase, oxidoreductase (Supplementary Additional File 8).
In heart, DEGs were enriched in 33 Reactome pathways, e.g., metabolism (FDR = 5.83E-12, fold enrichment = 1.97), signal transduction (FDR = 5.31E-03, fold enrichment = 1.45), transport of small molecules (FDR = 1.63E-04, fold enrichment = 2.06), gene expression (transcription) (FDR = 4.93E-02, fold enrichment = 1.61), and cell cycle (FDR = 4.81E-02, fold enrichment = 1.82) pathway in the reactome pathway category (Supplementary Additional File 6). Seven protein classes were enriched, mainly including transaminase, metalloprotease, oxidoreductase, and RNA binding protein (Supplementary Additional File 7).
It was observed that DEGs were significantly enriched in 15 categories of biological processes, including carboxylic acid metabolic process (GO:0019752; log10(P) = −5.98), monocarboxylic acid metabolic process (GO:0032787; log10(P) = −4.82), thrombocyte differentiation (GO:0002574; log10(P) = −4.41), circadian regulation of gene expression (GO:0032922; log10(P) = −4.41) and neuron projection development (GO:0031175; log10(P) = −4.18) in brain (Figure 7A).
FIGURE 7
In heart, DEGs were significantly enriched in 14 categories of biological processes. The most enriched pathways were carboxylic acid metabolic process (GO:0019752; log10(P) = −5.68), nucleotide metabolic process (GO:0009117; log10(P) = −5.64), heart development (GO:0007507; log10(P) = −5.26), generation of precursor metabolites and energy (GO:0006091; log10(P) = −5.18) and striated muscle tissue development (GO:0014706; log10(P) = −5.05) (Figure 7B). The most DEG genes were found in the oxidation-reduction process (GO:0055114) for both brain (111 genes) and heart tissues (76 genes) (Figures 7A,B).
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 of 11 genes were shared among 4 datasets, including ENSONIG00000005637, ENSONIG00000008480 ENSONIG 00000004328, ENSONIG00000017643, ENSONIG00000018614, ENSONIG00000019320, ENSONIG00000015833, ENSONIG0 0000014686 ENSONIG00000011498, ENSONIG00000011813, and ENSONIG0000000410. The information on overlapping between DAS and DEG datasets were presented in Figure 8.
FIGURE 8
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 CDS (Figure 9B). Both AS events caused a deficiency of a part of amino acid sequences in proteins (Figures 9A,B).
FIGURE 9
Discussion
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 (). In teleost, some studies highlighted the important role of AS in a few genes, e.g., thyroid hormone receptors in teleost fish (; ; ). However, few studies carried out genome-wide characterization of AS events in fish (; ). 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 (). 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 (). 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; ). 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 (; Wang et al., 2015). In the DAS gene list of catfish under heat stress, there were five SR genes (srsf2, srsf3, srsf5, srsf7-like, 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 (). 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 (). 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 (). When Drosophila adapted to seasonally cold days, a thermo-sensitive splicing event in 3′ UTR played an important role (). 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 (; ). The timeless (tim) gene in Drosophila melanogaster had two different transcript variants that were thermo-sensitively spliced (; ). In zebrafish, alternative promoter usage was detected for per3 gene after exposure to cold stress (16°C) for 24 h ().
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 (). 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 down-regulated 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 (). 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α (). 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 (; ). 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 (). 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.
Statements
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.
Funding
This work was supported by the Science and Technology Program of Guangzhou, China (Nos. 201804020013 and 201803020043), the National Natural Science Foundation of China (No. 31572612), and Modern Agriculture Talents Support Program (2016-2020).
Acknowledgments
We thank the reviewers for comments on the manuscript.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00244/full#supplementary-material
Abbreviations
- AE
alternative exon end
- AS
alternative splicing
- DAS
differentially alternative splicing
- DEGs
different expression genes
- FDR
false discovery rate
- hnrnp
heterogeneous nuclear RNP
- IR
retention of single intron
- miRNAs
microRNAs
- mRNA
messenger RNA
- MS-222
tricaine mesylate
- MSKIP
continuous exons skipping
- PGRPs
peptidoglycan recognition proteins
- RBPs
RNA-binding proteins
- SKIP
exon skipping
- SR
serine/arginine-rich
- TSS
transcription start site
- TTS
alternative transcription terminal site
- UTR
untranslated region
- XAE
alternative exon ends
- XIR, MIR XMIR
corresponding introns retention events
- XMSKIP
multiple continuous exons skipping
- XSKIP
multiple discontinuous exons skipping.
Footnotes
References
1
AliG. S.ReddyA. S. (2008). Regulation of alternative splicing of pre-mRNAs by stresses.Curr. Top. Microbiol. Immunol.326257–275.
2
AnandA.PatelM.LalremruataA.SinghA. P.AgrawalR.SinghL.et al (2008). Multiple alternative splicing of Dmrt1 during gonadogenesis in Indian mugger, a species exhibiting temperature-dependent sex determination.Gene42556–63. 10.1016/j.gene.2008.08.005
3
BaralleF. E.GiudiceJ. (2017). Alternative splicing as a regulator of development and tissue identity.Nat. Rev. Mol. Cell Biol.18437–451. 10.1038/nrm.2017.27
4
Ben-DovC.HartmannB.LundgrenJ.ValcarcelJ. (2008). Genome-wide analysis of alternative Pre-mRNA splicing.J. Biol. Chem.2831229–1233.
5
BoothroydC. E.WijnenH.NaefF.SaezL.YoungM. W. (2007). Integration of light and temperature in the regulation of circadian gene expression in Drosophila.PLoS Genet.3:e54. 10.1371/journal.pgen.0030054
6
BurrisT. P. (2008). Nuclear hormone receptors for heme: REV-ERBalpha and REV-ERBbeta are ligand-regulated components of the mammalian clock.Mol. Endocrinol.221509–1520. 10.1210/me.2007-0519
7
CanapleL.RambaudJ.Dkhissi-BenyahyaO.Rayet, Ba, TanN. S.et al (2006). Reciprocal regulation of brain and muscle Arnt-like protein 1 and peroxisome proliferator-activated receptor α defines a novel positive feedback loop in the rodent liver circadian clock.Mol. Endocrinol.201715–1727.
8
ChanC. B.ChengC. H. K. (2004). Identification and functional characterization of two alternatively spliced growth hormone secretagogue receptor transcripts from the pituitary of black seabream Acanthopagrus schlegeli.Mol. Cell. Endocrinol.21481–95.
9
ChangM. X.ZhangJ. (2017). Alternative Pre-mRNA splicing in mammals and teleost fish: a effective strategy for the regulation of immune responses against pathogen infection.Int. J. Mol. Sci.18:1530. 10.3390/ijms18071530
10
ChatchaiphanS.SrisapoomeP.KimJ. H.DevlinR. H.Na-NakornU. (2017). De Novo transcriptome characterization and growth-related gene expression profiling of diploid and triploid bighead catfish (Clarias macrocephalus Gunther, 1864).Mar. Biotechnol.1936–48. 10.1007/s10126-017-9730-3
11
ChenC. H.LiB. J.GuX. H.LinH. R.XiaJ. H. (2019). Marker-assisted selection of YY supermales from a genetically improved farmed tilapia-derived strain.Zool. Res.40108–112. 10.24272/j.issn.2095-8137.2018.071
12
ChengY.GvakhariaB.HardinP. E. (1998). Two alternatively spliced transcripts from the Drosophila period gene rescue rhythms having different molecular and behavioral characteristics.Mol. Cell. Biol.186505–6514.
13
ChristensenA. H.SharrockR. A.QuailP. H. (1992). Maize polyubiquitin genes: structure, thermal perturbation of expression and transcript splicing, and promoter activity following transfer to protoplasts by electroporation.Plant Mol. Biol.18675–689.
14
CongM.WuH. F.CaoT. F.LvJ. S.WangQ.JiC. L.et al (2018). Digital gene expression analysis in the gills of Ruditapes philippinarum exposed to short- and long-term exposures of ammonia nitrogen.Aquat. Toxicol.194121–131. 10.1016/j.aquatox.2017.11.012
15
ContrerasA. V.TorresN.TovarA. R. (2013). PPAR-alpha as a key nutritional and environmental sensor for metabolic adaptation.Adv. Nutr.4439–452. 10.3945/an.113.003798
16
DioumE. M.RutterJ.TuckermanJ. R.GonzalezG.Gilles-GonzalezM. A.McKnightS. L. (2002). NPAS2: a gas-responsive transcription factor.Science2982385–2387.
17
ElliottD. J.GrellscheldS. N. (2006). Alternative RNA splicing regulation in the testis.Reproduction132811–819.
18
El-SayedA. (2006). Environmental Requirements in Tilapia Culture.Wallingford: CABI Publishing.
19
FillatreauS.SixA.MagadanS.CastroR.SunyerJ. O.BoudinotP. (2013). The astonishing diversity of ig classes and B cell repertoires in teleost fish.Front. Immunol.4:28. 10.3389/fimmu.2013.00028
20
FloreaL.SongL.SalzbergS. L. (2013). Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues.F1000Res.2:188. 10.12688/f1000research.2-188.v2
21
FujikakeN.NagaiY.PopielH. A.KanoH.YamaguchiM.TodaT. (2005). Alternative splicing regulates the transcriptional activity of Drosophila heat shock transcription factor in response to heat/cold stress.FEBS Lett.5793842–3848.
22
GanH. M.AustinC.LintonS. (2018). Transcriptome-guided identification of carbohydrate active enzymes (cazy) from the christmas island red crab, Gecarcoidea natalis and a vote for the inclusion of transcriptome-derived crustacean cazys in comparative studies.Mar. Biotechnol.20654–665. 10.1007/s10126-018-9836-2
23
GemignaniF.SazaniP.MorcosP.KoleR. (2002). Temperature-dependent splicing of beta-globin pre-mRNA.Nucleic Acids Res.304592–4598.
24
GiannettoA.MaisanoM.CappelloT.OlivaS.ParrinoV.NatalottoA.et al (2017). Effects of oxygen availability on oxidative stress biomarkers in the Mediterranean mussel Mytilus galloprovincialis.Mar. Biotechnol.19614–626. 10.1007/s10126-017-9780-6
25
GuX. H.JiangD. L.HuangY.LiB. J.ChenC. H.LinH. R.et al (2018). Identifying a major QTL associated with salinity tolerance in Nile Tilapia using QTL-Seq.Mar. Biotechnol.2098–107. 10.1007/s10126-017-9790-4
26
HartleyS. W.MullikinJ. C. (2015). QoRTs: a comprehensive toolset for quality control and data processing of RNA-Seq experiments.BMC Bioinformatics16:224. 10.1186/s12859-015-0670-5
27
HartleyS. W.MullikinJ. C. (2016). Detection and visualization of differential splicing in RNA-Seq data with JunctionSeq.Nucleic Acids Res.44:e127. 10.1093/nar/gkw501
28
HeA. Y.LiuC. Z.ChenL. Q.NingL. J.ZhangM. L.LiE. C.et al (2014). Identification, characterization and nutritional regulation of two isoforms of acyl-coenzyme A oxidase 1 gene in Nile tilapia (Oreochromis niloticus).Gene54530–35. 10.1016/j.gene.2014.05.010
29
HopfN.Plesofsky-VigN.BramblR. (1992). The heat shock response of pollen and other tissues of maize.Plant Mol. Biol.19623–630.
30
HuB.JinJ. P.GuoA. Y.ZhangH.LuoJ. C.GaoG. (2015). GSDS 2.0: an upgraded gene feature visualization server.Bioinformatics311296–1297. 10.1093/bioinformatics/btu817
31
HuP.LiuM.LiuY.WangJ.ZhangD.NiuH.et al (2016). Transcriptome comparison reveals a genetic network regulating the lower temperature limit in fish.Sci. Rep.6:28952. 10.1038/srep28952
32
HuangB. Y.ZhangL. L.TangX. Y.ZhangG. F.LiL. (2016). Genome-wide analysis of alternative splicing provides insights into stress adaptation of the pacific oyster.Mar. Biotechnol.18598–609.
33
ImS.LeeH. N.JungH. S.YangS.ParkE. J.HwangM. S.et al (2017). Transcriptome-based identification of the desiccation response genes in marine red algae Pyropia tenera (Rhodophyta) and enhancement of abiotic stress tolerance by ptdrg2 in chlamydomonas.Mar. Biotechnol.19232–245. 10.1007/s10126-017-9744-x
34
JacobA. G.SinghR. K.ComiskeyD. F.RouhierM. F.MohammadF.BebeeT. W.et al (2014). Stress-induced alternative splice forms of mdm2 and mdmx modulate the p53-pathway in distinct ways.PLoS One9:e104444. 10.1371/journal.pone.0104444
35
JaksicA. M.SchlottererC. (2016). The interplay of temperature and genotype on patterns of alternative splicing in Drosophila melanogaster.Genetics204315–325.
36
JamesA. B.SyedN. H.BordageS.MarshallJ.NimmoG. A.JenkinsG. I.et al (2012). Alternative splicing mediates responses of the Arabidopsis circadian clock to temperature changes.Plant Cell24961–981. 10.1105/tpc.111.093948
37
JohnsonJ. M.CastleJ.Garrett-EngeleP.KanZ. Y.LoerchP. M.ArmourC. D.et al (2003). Genome-wide survey of human alternative pre-mRNA splicing with exon junction microarrays.Science3022141–2144.
38
JonssonV.OsterlundT.NermanO.KristianssonE. (2016). Statistical evaluation of methods for identification of differentially abundant genes in comparative metagenomics.BMC Genomics17:78. 10.1186/s12864-016-2386-y
39
JuX. H.XuH. J.YongY. H.AnL. L.XuY. M.JiaoP. R.et al (2014). Heat stress upregulates the expression of TLR4 and its alternative splicing variant in Bama miniature pigs.J. Integr. Agric.132479–2487.
40
JuZ.DunhamR. A.LiuZ. (2002). Differential gene expression in the brain of channel catfish (Ictalurus punctatus) in response to cold acclimation.Mol. Genet. Genomics26887–95.
41
KerstenS. (2014). Integrated physiology and systems biology of PPAR alpha.Mol. Metab.3354–371.
42
KiddP. B.YoungM. W.SiggiaE. D. (2015). Temperature compensation and temperature sensation in the circadian clock.Proc. Natl. Acad. Sci. U.S.A.112E6284–E6292. 10.1073/pnas.1511215112
43
KimY. H.KwonD. H.AhnH. M.KohY. H.LeeS. H. (2014). Induction of soluble AChE expression via alternative splicing by chemical stress in Drosophila melanogaster.Insect Biochem. Mol. Biol.4875–82. 10.1016/j.ibmb.2014.03.001
44
KucherenkoM. M.ShcherbataH. R. (2018). miRNA targeting and alternative splicing in the stress response - events hosted by membrane-less compartments.J. Cell Sci.131:jcs202002. 10.1242/jcs.202002
45
LeppekK.DasR.BarnaM. (2018). Functional 5 ’ UTR mRNA structures in eukaryotic translation regulation and how to find them.Nat. Rev. Mol. Cell Biol.19158–174.
46
LetunicI.DoerksT.BorkP. (2015). SMART: recent updates, new developments and status in 2015.Nucleic Acids Res.43D257–D260. 10.1093/nar/gku949
47
LiB. J.LiH. L.MengZ.ZhangY.LinH.YueG. H.et al (2017a). Copy number variations in tilapia genomes.Mar. Biotechnol.1911–21.
48
LiH. L.GuX. H.LiB. J.ChenX.LinH. R.XiaJ. H. (2017b). Characterization and functional analysis of hypoxia-inducible factor HIF1alpha and its inhibitor HIF1alphan in tilapia.PLoS One12:e0173478. 10.1371/journal.pone.0173478
49
LiH. L.GuX. H.LiB. J.ChenC. H.LinH. R.XiaJ. H. (2017c). Genome-wide QTL analysis identified significant associations between hypoxia tolerance and mutations in the GPR132 and ABCG4 genes in Nile Tilapia.Mar. Biotechnol.19441–453. 10.1007/s10126-017-9762-8
50
LiY. Q.LaiS. M.WangR. J.ZhaoY. C.QinH.JiangL. X.et al (2017d). RNA-Seq analysis of the antioxidant status and immune response of Portunus trituberculatus following aerial exposure.Mar. Biotechnol.1989–101. 10.1007/s10126-017-9731-2
51
LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al (2009). The sequence alignment/map format and SAMtools.Bioinformatics252078–2079. 10.1093/bioinformatics/btp352
52
LiuM.GuoX. M. (2017). A novel and stress adaptive alternative oxidase derived from alternative splicing of duplicated exon in oyster Crassostrea virginica.Sci. Rep.7:10785. 10.1038/s41598-017-10976-w
53
LongY.SongG. L.YanJ. J.HeX. Z.LiQ.CuiZ. B. (2013). Transcriptomic characterization of cold acclimation in larval zebrafish.BMC Genomics14:612. 10.1186/1471-2164-14-612
54
LuJ.PeatmanE.WangW.YangQ.AbernathyJ.WangS.et al (2010). Alternative splicing in teleost fish genomes: same-species and cross-species analysis and comparisons.Mol. Genet. Genomics283531–539. 10.1007/s00438-010-0538-3
55
MajercakJ.ChenW.-F.EderyI. (2004). Splicing of the period gene 3’-terminal intron is regulated by light, circadian clock factors, and phospholipase C.Mol. Cell. Biol.243359–3372.
56
MajercakJ.SidoteD.HardinP. E.EderyI. (1999). How a circadian clock adapts to seasonal decreases in temperature and day length.Neuron24219–230.
57
ManiatisT.TasicB. (2002). Alternative pre-mRNA splicing and proteome expansion in metazoans.Nature418236–243.
58
MarchandO.SafiR.EscrivaH.Van RompaeyE.PrunetP.LaudetV. (2001). Molecular cloning and characterization of thyroid hormone receptors in teleost fish.J. Mol. Endocrinol.2651–65.
59
Martinez-MontielN.Rosas-MurrietaN. H.Anaya RuizM.Monjaraz-GuzmanE.Martinez-ContrerasR. (2018). Alternative splicing as a target for cancer treatment.Int. J. Mol. Sci.19:545.
60
MastrangeloA. M.MaroneD.LaidòG.De LeonardisA. M.De VitaP. (2012). Alternative splicing: enhancing ability to cope with stress via transcriptome plasticity.Plant Sci.18540–49. 10.1016/j.plantsci.2011.09.006
61
MazzucotelliE.MastrangeloA. A.CrosattiC.GuerraD.StancaA. M.CattivelliL. (2008). Abiotic stress response in plants: when post-transcriptional and post-translational regulations control transcription.Plant Sci.174420–431. 10.3390/ijms19123737
62
MelangathliG.SenT.KumarR.BawaP.SrinivasanS.VijayraghavanU. (2017). Functions for fission yeast splicing factors SpSlu7 and SpPrp18 in alternative splice-site choice and stress-specific regulated splicing.PLoS One12:e0188159. 10.1371/journal.pone.0188159
63
MontelliS.MazzottaG.VaninS.CaccinL.CorraS.De PittaC.et al (2015). period and timeless mRNA splicing profiles under natural conditions in Drosophila melanogaster.J. Biol. Rhythms30217–227. 10.1177/0748730415583575
64
MooreL. J.SomamotoT.LieK. K.DijkstraJ. M.HordvikI. (2005). Characterisation of salmon and trout CD8 alpha and CD8 beta.Mol. Immunol.421225–1234.
65
NilsenT. W.GraveleyB. R. (2010). Expansion of the eukaryotic proteome by alternative splicing.Nature463457–463. 10.1038/nature08909
66
PoseD.VerhageL.OttF.YantL.MathieuJ.AngenentG. C.et al (2013). Temperature-dependent regulation of flowering by antagonistic FLM variants.Nature503414–417. 10.1038/nature12633
67
RuanK.ZhuY.LiC.BrazillJ. M.ZhaiR. G. (2015). Alternative splicing of Drosophila Nmnat functions as a switch to enhance neuroprotection under stress.Nat. Commun.6:10057. 10.1038/ncomms10057
68
SanoY.ShiinaT.NaitouK.NakamoriH.ShimizuY. (2015). Hibernation-specific alternative splicing of the mRNA encoding cold-inducible RNA-binding protein in the hearts of hamsters.Biochem. Biophys. Res. Commun.462322–325. 10.1016/j.bbrc.2015.04.135
69
ScottiM. M.SwansonM. S. (2016). RNA mis-splicing in disease.Nat. Rev. Genet.1719–32.
70
SeokH.ChaB. S. (2013). Refocusing peroxisome proliferator activated receptor-α: a new insight for therapeutic roles in diabetes.Diabetes Metab. J.37326–332.
71
ShiY.LiuW. G.HeM. X. (2018). Proteome and transcriptome analysis of ovary, intersex gonads, and testis reveals potential key sex reversal/differentiation genes and mechanism in scallop Chlamys nobilis.Mar. Biotechnol.20220–245. 10.1007/s10126-018-9800-1
72
SorekR. (2007). The birth of new exons: mechanisms and evolutionary consequences.RNA131603–1608.
73
SteriM.IddaM. L.WhalenM. B.OrruV. (2018). Genetic variants in mRNA untranslated regions.Wiley Interdiscip. Rev. RNA9:e1474. 10.1002/wrna.1474
74
SunY. B. (2017). FasParser: a package for manipulating sequence data.Zool. Res.38110–112. 10.24272/j.issn.2095-8137.2017.017
75
TakedaY.ChumanY.ShirasuN.SatoS.MatsushimaA.KanekiA.et al (2004). Structural analysis and identification of novel isoforms of the circadian clock gene period in the silk moth Bombyx mori.Zool. Sci.21903–915.
76
TanS.WangW.TianC.NiuD.ZhouT.JinY.et al (2019). Heat stress induced alternative splicing in catfish as determined by transcriptome analysis.Comp. Biochem. Physiol. Part D Genomics Proteomics29166–172. 10.1016/j.cbd.2018.11.008
77
TanS. X.WangW. W.ZhongX. X.TianC. X.NiuD. H.BaoL. S.et al (2018). Increased alternative splicing as a host response to Edwardsiella ictaluri infection in catfish.Mar. Biotechnol.20729–738. 10.1007/s10126-018-9844-2
78
TangL. T.RanX. Q.MaoN.ZhangF. P.NiuX.RuanY. Q.et al (2018). Analysis of alternative splicing events by RNA sequencing in the ovaries of Xiang pig at estrous and diestrous.Theriogenology11960–68. 10.1016/j.theriogenology.2018.06.022
79
TrapnellC.HendricksonD. G.SauvageauM.GoffL.RinnJ. L.PachterL. (2013). Differential analysis of gene regulation at transcript resolution with RNA-seq.Nat. Biotechnol.3146–53.
80
TrapnellC.PachterL.SalzbergS. L. (2009). TopHat: discovering splice junctions with RNA-Seq.Bioinformatics251105–1111. 10.1093/bioinformatics/btp120
81
TrapnellC.RobertsA.GoffL.PerteaG.KimD.KelleyD. R.et al (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks.Nat. Protoc.7562–578. 10.1038/nprot.2012.016
82
TrapnellC.WilliamsB. A.PerteaG.MortazaviA.KwanG.van BarenM. J.et al (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation.Nat. Biotechnol.28511–515. 10.1038/nbt.1621
83
TripathiS.PohlM. O.ZhouY. Y.Rodriguez-FrandsenA.WangG. J.SteinD. A.et al (2015). Meta- and orthogonal integration of influenza “OMICs” data defines a role for UBR4 in virus budding.Cell Host Microbe18723–735. 10.1016/j.chom.2015.11.002
84
TsalafoutaA.SarropoulouE.PapandroulakisN.PavlidisM. (2018). Characterization and expression dynamics of key genes involved in the gilthead sea bream (Sparus aurata) cortisol stress response during early ontogeny.Mar. Biotechnol.20611–622. 10.1007/s10126-018-9833-5
85
VenturaP.ToullecG.FricanoC.ChapronL.MeunierV.RottingerE.et al (2018). Cnidarian primary cell culture as a tool to investigate the effect of thermal stress at cellular level.Mar. Biotechnol.20144–154. 10.1007/s10126-017-9791-3
86
VornanenM.HassinenM.KoskinenH.KrasnovA. (2005). Steady-state effects of temperature acclimation on the transcriptome of the rainbow trout heart.Am. J. Physiol. Regul. Integr. Comp. Physiol.289R1177–R1184.
87
WanQ. Y.SuJ. G. (2015). Transcriptome analysis provides insights into the regulatory function of alternative splicing in antiviral immunity in grass carp (Ctenopharyngodon idella).Sci. Rep.5:12946. 10.1038/srep12946
88
WangE. T.SandbergR.LuoS. J.KhrebtukovaI.ZhangL.MayrC.et al (2008). Alternative isoform regulation in human tissue transcriptomes.Nature456470–476. 10.1038/nature07509
89
WangL.BaiB.HuangS. Q.LiuP.WanZ. Y.YeB. Q.et al (2017). QTL mapping for resistance to iridovirus in Asian seabass using genotyping-by-sequencing.Mar. Biotechnol.19517–527. 10.1007/s10126-017-9770-8
90
WangY.LiuJ.HuangB.XuY. M.LiJ.HuangL. F.et al (2015). Mechanism of alternative splicing and its regulation.Biomed. Rep.3152–158.
91
WangZ. C.CuiJ.SongJ.WangH. Z.GaoK. L.QiuX. M.et al (2018). Comparative transcriptome analysis reveals growth-related genes in juvenile Chinese sea cucumber, Russian sea cucumber, and their hybrids.Mar. Biotechnol.20193–205. 10.1007/s10126-018-9796-6
92
XiaJ. H.BaiZ. Y.MengZ. N.ZhangY.WangL.LiuF.et al (2015). Signatures of selection in tilapia revealed by whole genome resequencing.Sci. Rep.5: 14168. 10.1038/srep14168
93
XiaJ. H.LiH. L.LiB. J.GuX. H.LinH. R. (2018). Acute hypoxia stress induced abundant differential expression genes and alternative splicing events in heart of tilapia.Gene63952–61. 10.1016/j.gene.2017.10.002
94
XuZ.GanL.LiT.XuC.ChenK.WangX.et al (2015). Transcriptome profiling and molecular pathway analysis of genes in association with salinity adaptation in Nile Tilapia Oreochromis niloticus.PLoS One10:e0136506. 10.1371/journal.pone.0136506
95
YamaguchiS.IwasaY. (2018). Temperature-dependent sex determination, realized by hormonal dynamics with enzymatic reactions sensitive to ambient temperature.J. Theor. Biol.453146–155. 10.1016/j.jtbi.2018.05.023
96
YangC. G.JiangM.WenH.TianJ.LiuW.WuF.et al (2015). Analysis of differential gene expression under low-temperature stress in Nile tilapia (Oreochromis niloticus) using digital gene expression.Gene564134–140. 10.1016/j.gene.2015.01.038
97
YangX. L.IkhwanuddinM.LiX. C.LinF.WuQ. Y.ZhangY. L.et al (2018). Comparative transcriptome analysis provides insights into differentially expressed genes and long non-coding RNAs between ovary and testis of the mud crab (Scylla paramamosain).Mar. Biotechnol.2020–34. 10.1007/s10126-017-9784-2
98
YeH.XiaoS. J.WangX. Q.WangZ. Y.ZhangZ. S.ZhuC. K.et al (2018). Characterization of spleen transcriptome of Schizothorax prenanti during Aeromonas hydrophila infection.Mar. Biotechnol.20246–256. 10.1007/s10126-018-9801-0
99
YeoG.HolsteD.KreimanG.BurgeC. B. (2004). Variation in alternative splicing across human tissues.Genome Biol.5:R74.
100
YinL.WuN.LazarM. A. (2010). Nuclear receptor Rev-erbalpha: a heme receptor that coordinates circadian rhythm and metabolism.Nucl. Recept. Signal.8:e001. 10.1621/nrs.08001
101
YuL. Y.XuD. D.YeH.YueH. M.OokaS.KondoH.et al (2018). Gonadal transcriptome analysis of pacific abalone Haliotis discus discus: identification of genes involved in germ cell development.Mar. Biotechnol.20467–480. 10.1007/s10126-018-9809-5
102
YueC. Y.LiQ.YuH. (2018). Gonad transcriptome analysis of the pacific oyster Crassostrea gigas identifies potential genes regulating the sex determination and differentiation process.Mar. Biotechnol.20206–219. 10.1007/s10126-018-9798-4
103
ZambonA. C.McDearmonE. L.SalomonisN.VranizanK. M.JohansenK. L.AdeyD.et al (2003). Time- and exercise-dependent gene regulation in human skeletal muscle.Genome Biol.4:R61.
104
ZhangQ. M.ZhaoX.LiZ.WuM.GuiJ. F.ZhangY. B. (2018). Alternative splicing transcripts of zebrafish LGP2 gene differentially contribute to IFN antiviral response.J. Immunol.200688–703. 10.4049/jimmunol.1701388
105
ZhaoY.WangJ.ThammaratsuntornJ.WuJ. W.WeiJ. H.WangY.et al (2015). Comparative transcriptome analysis of Nile tilapia (Oreochromis niloticus) in response to alkalinity stress.Genet. Mol. Res.1417916–17926. 10.4238/2015.December.22.16
106
ZhouT.GuiL.LiuM.LiW.HuP.DuarteD. F. C.et al (2018). Transcriptomic responses to low temperature stress in the Nile tilapia, Oreochromis niloticus.Fish Shellfish Immunol.841145–1156. 10.1016/j.fsi.2018.10.023
107
ZhouZ.LiM.ChengH.FanW. L.YuanZ. R.GaoQ.et al (2018). An intercross population study reveals genes associated with body size and plumage color in ducks.Nat. Commun.9:2648.
108
ZwolloP. (2011). Dissecting teleost B cell differentiation using transcription factors.Dev. Comp. Immunol.35898–905. 10.1016/j.dci.2011.01.009
Summary
Keywords
alternative splicing, tilapia, cold stress, RNA-seq, genome
Citation
Li BJ, Zhu ZX, Qin H, Meng ZN, Lin HR and Xia JH (2020) Genome-Wide Characterization of Alternative Splicing Events and Their Responses to Cold Stress in Tilapia. Front. Genet. 11:244. doi: 10.3389/fgene.2020.00244
Received
26 November 2019
Accepted
28 February 2020
Published
18 March 2020
Volume
11 - 2020
Edited by
Shikai Liu, Ocean University of China, China
Reviewed by
Chen Jiang, Dalian Ocean University, China; Yulin Jin, Emory University, United States
Updates
Copyright
© 2020 Li, Zhu, Qin, Meng, Lin and Xia.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jun Hong Xia, xiajunh3@mail.sysu.edu.cn; flyxia936@163.com
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.