Introduction

The recent advent of next generation RNA sequencing (RNA-Seq) and publication of reference genomes for many organisms have allowed researchers to study transcription profiles. The majority of the mammalian genome (up to 80%) is transcribed, of which only 2–3% are protein-coding RNAs (mRNAs) and the rest are noncoding RNAs. Most transcribed ncRNAs are larger than 200 nucleotides and defined as long non-coding RNAs (lncRNA) with some further processed to generate small RNAs1. LncRNA have fewer, longer exons when compared to coding genes and exhibit cell-type specific expression. These long non-coding RNAs are reported to play critical roles in various biological processes, including chromatin modification, regulation of transcription, influence of nuclear architecture and regulation of gene expression at post-transcriptional and post-translational levels, and are also reported to interact with DNA, RNA and proteins during these processes2. Availability of sequencing techniques led to the discovery of non-coding transcripts involved in regulation of gene expression, as well as novel protein coding genes.

Studies in humans (Homo sapiens) and mice (Mus musculus) support that lncRNA are involved in transcriptional regulation of skeletal muscle development and differentiation2. For instance, ChiP-Seq studies identified two different lncRNAs (DRRRNA and CERNA) that transcribe from the enhancer region of myogenic transcription factor MyoD2,3 to increase expression of MyoD and myogenin. An additional study identified a number of muscle specific lncRNAs regulated by Yin Yang 1 protein that function as both positive or negative regulators of muscle differentiation4. Expanding the investigation of lncRNAs as myogenic regulators in additional species is critical to both 1) establish whether previously reported mechanisms are conserved and 2) identify novel mechanisms within a species.

Rainbow trout (Oncorhynchus mykiss) is a fish species valued as a research model, food source, and recreational fish. This species demonstrates a pattern of indeterminate growth, or the ability to grow in length and continually accumulate muscle throughout its lifespan5. Muscle regulatory factors (MRF) such as MyoD and myogenin that regulate hyperplasia and hypertrophy in mammalian muscle are also important for these processes in salmonids6. However, unique MRF expression patterns occur during differentiation of rainbow trout myogenic precursor cells (MPCs)7, suggesting that these mechanisms are significant for continued hyperplasic capacity and indeterminate growth.

Determining how lncRNAs are regulated by biological factors that affect muscle growth is central for identifying lncRNA-related mechanisms important for myogenesis. In salmonids, estradiol (E2) is a biological factor that negatively effects muscle growth8,9, partially through direct effects on MPC proliferation10,11 and protein turnover12,13. Plasma E2 levels increase during sexual maturation and its catabolic effects in muscle are important for directing nutrients away from muscle growth and toward gonad development8,14, an energy intensive process that compromises muscle quality15,16. Our previous findings support a role for noncoding RNAs in the short-term E2-response (24 hr post-injection), specifically through regulation of miRNAs that directly affect expression of genes related to proteolysis, cell differentiation, and cell proliferation11. The present study aims to identify the differentially expressed lncRNAs and mRNAs in skeletal muscle exposed to E2 and create a network of the mRNAs and lncRNAs involved in the process. Here we use an RNA-Seq approach to comprehensively investigate temporal gene expression differences in skeletal muscle of rainbow trout after a single injection of E2. We also construct a functional lncRNA-mRNA regulatory network that is associated with skeletal muscle response to E2 treatment, and conduct a pathway enrichment analysis. This comprehensive analysis helps identify novel lncRNA and their related mechanisms with emphasis on regulation of muscle growth in rainbow trout under the influence of E2.

Results

Generation of a muscle transcriptome reference

To provide a comprehensive understanding of the effects of E2 treatment on muscle physiology, the lncRNA and mRNA transcriptomes were analyzed in skeletal muscle from E2-treated and control rainbow trout 24 and 72 hours post-injection. A total of 789,485,036 paired-end raw reads were generated from 16 samples with 101-bp read length (n = 4 per time and treatment combination). The number of sequences from each sample ranged from 36.2 to 63.8 million. After quality control including removal of ambiguous nucleotides, low-quality bases and ribosomal RNA sequences, a total of 749,490,934 cleaned reads (94.9%) were harvested for further analysis. The number of cleaned reads of each sample ranged from 34.5 to 60.9 million (Table 1).

Table 1 Summary of samples and RNA-Seq data. CTRL indicates control group. E2 represents estradiol treatment group.
Full size table

The cleaned reads were pooled and assembled by Trinity17. Then CD-HIT-EST18 was used to remove the redundancy. As shown in Table 2, the transcriptomes were assembled into 243,509 contigs (203,148 Trinity components), ranging from 201 to 20,635 bp in length. To provide an expressed transcriptome reference and filter out transcripts that had very low read counts, EdgeR19 was used to remove the transcripts with count-per-million (CPM) less than 1. This produced 63,181 contigs (31,419 Trinity components), again ranging from 201 to 20,635 bp in length. The average length was 1,466 bp, N50 length was 1,982 bp and median length was 1,189 bp (Table 2). All cleaned reads were mapped to the expressed transcriptome reference by using the ultrafast short read aligner Bowtie20. The mapping ratio ranged from 84.3% to 89.8% with an average of 87.7% (Table 1).

Table 2 Statistics of transcriptome assembly.
Full size table

Identification of differentially expressed genes

The pipeline reported in our previous study21 was used to identify all lncRNAs from the expressed transcriptome reference. Two biological replicates (CTRL3 & EST4) at 24 hour time point were identified as outliers based upon multidimensional distance scaling analysis and therefore were excluded. Finally, differential expression analyses were performed of coding genes and lncRNAs in E2 treated fish at 24 hours and 72 hours compared with the respective control fish. Given a false discovery rate (FDR) of 5%, 479 (226 lncRNA and 253 mRNA) differentially expressed genes (DEGs) were detected between E2 and control fish at 24 hours while only 19 DEGs (9 lncRNA and 10 mRNA) were detected at 72 hours (Supplementary Table S1 & S2). Only one lncRNA and one mRNA were differentially expressed at both time points (Fig. 1A,B). Due to the transient nature of the response, an extensive analysis of DEGs were completed only on differentially expressed lncRNA and mRNA detected 24 hour post-injection.

Figure 1
figure 1

The number of differentially expressed genes in rainbow trout treated with estradiol (E2). (A) Venn diagram of common differentially expressed mRNA in E2 treated fish between 2 time points (24 hours post-injection versus 72 hours post-injection). (B) Venn diagram of common differentially expressed lncRNA in E2 treated fish between 2 time points. (C) A volcano plot of differentially expressed transcripts (lncRNAs and mRNAs) between control and E2 treated group. (D–E) Unsupervised hierarchical clustering of the expression profiles of differentially expressed mRNAs (D) and lncRNAs (E) both distinguish E2 treated and control group. Up- and down-regulated genes are indicated respectively by red and blue colors.

Full size image

Transcriptome analysis identified 102 up-regulated and 151 down-regulated mRNAs and 119 up-regulated and 117 down-regulated lncRNAs in fish 24 hours post-E2 injection (Fig. 1C–E). Heat maps were generated using differentially expressed mRNAs (Fig. 1D) and lncRNAs (Fig. 1E) that self-segregated into control (CTRL) and EST (E2 treated) clusters in an unsupervised hierarchical clustering analysis. There were 36 genes that were differentially expressed more than 16-fold between control and E2 treated samples, of which 14 were up-regulated and 12 were down-regulated by E2.

Experimental validation of lncRNA and mRNA

To confirm the sequencing results, seven lncRNAs and five mRNAs that were differentially expressed were randomly selected for validation using real time RT-PCR. Each transcript was confirmed as differentially expressed in E2-treated fish in the same direction as indicated by the RNAseq analysis (Fig. 2A,B).

Figure 2
figure 2

Differentially expressed lncRNAs and mRNAs validated by qRT-PCR. Comparison between RNA-seq and qRT-PCR validation results. X-axis shows lncRNAs (A) and mRNAs (B) validated in this study, relative expression normalized to β-actin; Y-axis shows log2Ratio of expression of E2 versus control. Stars indicate significant difference based on t-test (p < 0.05).

Full size image

Differentially expressed genes represent several important pathways

To obtain the functions of the differentially expressed mRNAs and connections among them, we performed GO term and KEGG pathway enrichment analyses. To comprehensively interpret the GO term enrichment result, treemap was constructed based on the result of semantic similarity analysis. In each treemap, related terms were joined into a ‘supercluster’ with the same color and the most significant term as the representative of the group. As shown in Fig. 3a,b, in response to E2, skeletal muscle tissue showed notable up-regulations related to oxidative stress response, hormone response, protein ubiquitination, cysteine biosynthesis and DNA repair, and down-regulations mostly related to mesenchyme morphogenesis and cAMP biosynthesis in biological process. Notably, the term of response to leptin (GO:0044321) was identified to be enriched (Supplementary Fig. S1). Up-regulated genes in molecular function relate to kinase activity, hormone receptor activity, amino-acyl transferase activity and ubiquitin conjugating enzyme activity. Down regulated genes related to molecular function were related to muscle structural protein binding, motor activity, pyruvate carboxykinase activity and glucose-6-phosphate activity (Supplementary Fig. S2a,b). The enrichment of up-regulated genes in cellular component largely related to chromosome and ubiquitin conjugating enzyme complex, and down-regulated genes related to skeletal muscle structural protein complex (Supplementary Fig. S3a,b).

Figure 3
figure 3

GO analysis of differentially expressed mRNAs. (a) Treemap of up-regulated representative GO terms in biological processes. (b) Treemap of down-regulated representative GO terms in biological processes. Each rectangle is a single cluster representative. The representative GO terms were joined into ‘superclusters’ of loosely related terms with same color.

Full size image

The KEGG pathway analyses mapped the DEGs to KEGG reference pathways to infer systemic biological behaviors. Comparing to the control group, we observed significant KEGG pathway enrichment of DEGs in skeletal muscle in response to E2 stress. Of the top 20 over-represented KEGG pathways (Fig. 4), four were involved in signal transduction, including Rap1 signaling pathway, Jak-STAT signaling pathway, calcium signaling pathway and AMPK signaling pathway. An additional four pathways were associated with the endocrine system, which are PPAR signaling pathway, insulin resistance, estrogen signaling pathway and adipocytokine signaling pathway. Finally, the remaining 12 over-represented KEGG pathways were classified into different functional groups, including amino acid and carbohydrate metabolism, cellular community, protein synthesis and digestion.

Figure 4
figure 4

KEGG pathway analysis of differentially expressed genes. The rich factor was calculated using the gene count divided by the expected gene count.

Full size image

Construction of lncRNA-mRNA co-expression network

To construct the differentially expressed lncRNA-mRNA co-expression network, the normalized expression values of the lncRNAs and mRNAs were obtained by using DESeq222. Subsequently, Pearson’s correlation coefficient (PCC) was calculated between the normalized expression values of each of the lncRNA-mRNA pairs. The lncRNA-mRNA pairs with a PCC greater than 0.99 and FDR less than 0.05 were selected for network construction. The co-expression network consisted of 681 co-expression relationships between 164 lncRNAs and 201 mRNAs (Fig. 5). We then considered the node degree of the network, as a higher degree indicated that the nodes were likely to be hubs and therefore involved in more competing interactions. As a more stringent threshold, the top 5% (top 20) of the nodes were defined as hubs (Table 3). These 20 hub nodes that contain 14 lncRNAs co-expressed with more than 40% of the nodes in the network, implying the centrality of these nodes.

Figure 5
figure 5

lncRNA-mRNA co-expression network. The differentially expressed lncRNA-mRNA co-expression network consisted of 681 co-expression relationship between 164 lncRNAs and 201 mRNAs. Blue circles represent lncRNAs, and red circles denote mRNAs.

Full size image
Table 3 Top 20 nodes with the highest degree in differential lncRNA-mRNA co-expression network.
Full size table

A lncRNA-pathway network was constructed to identify the possible functions of lncRNAs, in which nodes represent lncRNAs or over-represented pathways. Connections between nodes were made if mRNAs that co-express with a lncRNA associate with a common pathway, thus supporting that the pathway is regulated by the corresponding lncRNAs (Fig. 6). Sixty-five lncRNAs (Supplementary Table S3) were linked with 20 enriched pathways (Table 4, Supplementary Table S4) in the lncRNA-pathway network, suggesting these lncRNAs play a central role in regulating the E2 response.

Figure 6
figure 6

lncRNA-pathway network. Blue circles represent lncRNAs, and orange circles denote pathways.

Full size image
Table 4 Key pathways in E2 treated skeletal muscle.
Full size table

Discussion

In rainbow trout, E2 is a maturation-related signal that negatively affects muscle growth. The present study identified lncRNAs and mRNAs regulated by E2 in rainbow trout skeletal muscle, thereby identifying mechanisms contributing to E2-induced muscle catabolism. Association analysis suggests that lncRNAs are central regulators of multiple components within a single pathway, supporting lncRNAs as important mechanisms for muscle physiology and the E2 response.

LncRNAs are emerging as regulators of diverse biological functions23. We identified 226 lncRNAs that were differentially expressed, suggesting that lncRNAs are regulated by E2 and likely contribute to the catabolic effects of this steroid in skeletal muscle. Only three of the differentially expressed lncRNAs (lnc-OM206, 3272 & 7282) were previously reported in a rainbow trout lncRNA reference transcriptome pooled from 15 tissues, including muscle21. This lack of overlap suggests that lncRNAs are temporospatially expressed, introducing the need for greater coverage of the rainbow trout lncRNA transcriptome across various tissues, periods of development, and treatment conditions. In particular, the sensitivity of lncRNAs (and mRNAs) to treatment conditions are supported by this study with the dramatic reduction in number of DEGs between 24 and 72 hr post-injection.

The functional roles of lncRNAs are poorly characterized. Co-expression models assist in large-scale prediction of lncRNA functions through integrating expression profiles of protein-coding genes and lncRNAs (coding-non-coding gene co-expression network)24,25. In this study a similar co-expression analysis was used to predict lncRNA functions based on their association with co-expressed genes, based on the fact that genes with similar expression patterns have a greater tendency to be involved in the same pathways26. The resultant co-expression network contained a total of 164 differentially expressed lncRNAs and 201 differentially expressed mRNAs. This network provides a global view of possible lncRNA-coding gene expression associations that are influenced by E2 treatment in salmonids. In this network, the top 5% (20) large degree nodes were all lncRNAs, of which, 14 lncRNAs co-expressed with more than 40% of the nodes. These observations suggest that lncRNAs have functional significance in the E2-induced reduction of muscle growth in rainbow trout.

To examine the key lncRNAs and their potential functions, lncRNA-pathway network was constructed based on pathway enrichment analysis. The lncRNA-pathway network contains 65 lncRNAs linked with 20 significantly enriched pathways. The top pathway containing the most lncRNAs linked was focal adhesion (FA). Focal adhesions are integrin-containing, multi-protein structures that mediate the regulatory effects of a cell in response to extracellular matrix (ECM) adhesion27 which also contribute to ECM degradation28. Our observation that E2 induces regulation of lncRNAs associating with FAs may contribute to changes in the structure and functional capacity of the skeletal muscle ECM. In mammals skeletal muscle ECM adapts to conditions that affect muscle growth such as injury or disease29, conditions that produce similar negative effects on muscle growth induced by E2 in salmonids8,30. Additionally, spawning female rainbow trout exhibiting muscle atrophy display reduced expression of ECM and FA-related genes31, supporting that regulation of related proteins represents muscle restructuring.

In addition to FA, we also found that calcium signaling pathway was linked with 15 lncRNAs in our lncRNA-pathway network. Calcium ions are important for cellular signaling, which plays a pivotal role in almost all cellular processes. Calpain is a protein belonging to the family of non-lysosomal cysteine proteases that is activated by calcium ions. Calpains function in various biological processes in skeletal muscle, including apoptosis32, proteolysis of myofibrillar proteins33, and myogenic events34. Increased caplain expression occurs concurrently with elevated E2 during spawning in rainbow trout35, however, regulation of calpastatins, the calpain inhibitors, also occurs during this period36, which is under direct effect of E212. Therefore, regulation of lncRNAs affecting calcium signaling capacity is likely significant for E2-induced elevations in muscle protein degradation. Probably further affecting calcium signaling was down-regulation of calsequestrin, a calcium-binding protein that regulates calcium homeostasis37. Our previous findings suggest that target genes of miR-23a-3p are involved in mitochondrial outer membrane permeability11; calcium signaling also plays a vital role in mitochondrial permeability which results in cell death. Additional lncRNAs regulated by E2 in skeletal muscle have functional roles in nutrient metabolism, including amino acids, fatty acids, and pyruvate, as well as protein degradation. As E2 shifts nutrient partitioning from muscle accretion in support of gonad development, the energy requirements of muscle are expected to change. These data support a role for lncRNA as regulators of these processes. Although the results of the present study require further experimental verification, they provide insight into the significant role that lncRNAs have in the physiological and biochemical response to E2 in skeletal muscle of rainbow trout.

Various genes functioning as molecular factors, cellular components and those involved in biological processes were up-regulated in rainbow trout skeletal muscle under the influence of E2. These proteins function to regulate different signaling pathways including JAK/STAT pathway, response to oxidative stress, cell-cell signaling and intracellular estrogen receptor signaling pathways with positive regulation of various receptor biosynthesis and their binding activity (insulin receptor binding, insulin like growth factor binding, integrin binding etc). Besides, important protein cross-linking molecules involved in histone ubiquitination, protein autoubiquitination, and their responsible enzymes like ubiquitin conjugating enzyme and its complex were also up-regulated. Elevated expression of important metabolic enzymes, including aldehyde dehydrogenase and succinate dehydrogenase, important for energy synthesis was observed.

Estradiol induces the activation of estrogen receptors that are functional through genomic and non-genomic pathways38. Besides activation of estrogen receptors, E2 is recognized as a regulator of the IGF pathway in salmonids39. IGFBP-5, an IGF binding protein, was reported to be highly conserved among humans, rats (Rattus norvegicus), mice and frog (Xenpus laevis), and is abundantly expressed in skeletal muscle during myoblast differentiation in mouse40. This response inhibits muscle cell differentiation by binding to the IGF receptor and impeding IGF activity41. In the current study, real time PCR and RNAseq analysis both support increased IFGBP-5 expression after 24 hours of E2 treatment, and a similar response was observed in a previous rainbow trout E2 injection study14. However, long-term E2 exposure (30 days) in rainbow trout had the opposite effect8. The differential time-dependent response could reflect the complex and dynamic nature of the IGF endocrine system, in which relative concentrations of both IGF and IGF binding proteins regulate IGF binding to receptors. Collectively these observations support that regulation of IGFBP-5 in muscle is important for effects of E2 on muscle cell differentiation.

Production of reactive oxygen species leading to oxidative stress as a response to E2 was reported in Japanese medaka (Oryzias latipes) and in Japanese sea bass (Lateolabrax japonicus)42,43,44. DNA damage induced by reactive oxygen species was reported in fathead minnows (Pimephales promelas) and Japanese sea bass exposed to E244,45. Reactive oxygen species (ROS) contribute to autophagy46 and muscle atrophy by activating muscle degrading pathways involving calpains, caspase-3, ATG4b, MuRF-1 and atrogin-147. Supporting this concept is our previous report of increased expression of caspase-9, caspase-3 and atrogin-1 by E2 in rainbow trout11. In vivo and in vitro experiments in rainbow trout and rainbow trout primary myocytes showed that exposure to E2 increased expression of atrogin-1/fbxo32 and MURF genes and increased protein degradation11,48. Atrogin-1 is a ubiquitin ligase responsible for ubiquitinating specific proteins such as MyoD and myogenin for proteasomal degradation49,50, and increased atrogin-1 expression is typically associated with catabolic conditions in both fish and mammals51,52. Supporting E2-induced reductions in muscle protein retention in rainbow trout was increased expression of autophagy related 4b cysteine peptidase (ATG4b) and previous literature reports a similar response for additional components of the autophagy-lysosome system8,12. Collectively, these findings indicate the regulation of E2 on the skeletal myogenic program and mechanisms regulating protein retention and muscle accretion.

Characterized GO terms that showed reduced expression in E2-treated skeletal muscle after 24 hours are mainly involved in structural make up of skeletal muscle. Significantly down regulated proteins are primarily involved in actin cytoskeleton, actin binding, motor activity and mesenchymal morphogenesis. Additionally, GO terms related to cAMP biosynthesis, cristae formation, regulation of smoothened signaling pathway, and RNA mediated DNA polymerase activity were differentially regulated. In humans and mouse mesenchymal stromal cells are multipotent cells that differentiate to different cell lineages including mesodermal cells considered as precursors for myoblasts or MPCs53,54. The bone marrow derived mesenchymal cells with myogenic markers migrate to sites of muscle regeneration55,56. In the current study, the observed decrease in the genes responsible for mesenchymal cell migration and morphogenesis may have negative effects on muscle development. Wnt signaling in rat is necessary for differentiation of mesenchymal cells to myogenic origin57, which also inhibits their differentiation to adipogenic origin by decreasing the expression of CCAAT enhancer binding protein alpha and peroxisome proliferator-activated receptor gamma (PPARγ)57. Peroxisome proliferator-activated receptor gamma coactivator 1 alpha (PPARGC1A) is a coactivator of PPARγ and highly expressed in mouse skeletal muscle58,59. Studies also indicate that PPARGC1A and PPARGC1B play prominent roles in estrogen receptor signaling pathway by acting as cofactors to enhance transactivation of estrogen receptor α60,61,62. In our study, PPARGC1A exhibited increased expression in skeletal muscle within 24 hours of E2 exposure, supporting that effects of E2 in muscle extend to regulating mechanisms affecting mesenchymal cell differentiation and fiber type.

In summary, our comprehensive analyses provide novel knowledge of mRNAs and lncRNAs at the transcriptomic level regarding the influence of E2 on rainbow trout skeletal muscle within a relatively short period after steroid injection (24 hr). These results and conclusions will serve as important resources for future experiments that further investigate the role and regulation of lncRNAs in rainbow trout. The dynamic expression response observed in the current study may be driven by a dramatic reduction in plasma E2 between 24 and 72 hrs post-injection due to high E2 clearance. A previous study using an identical injection approach reported a greater than 3-fold reduction in plasma E2 between these sampling periods (121 ng/ml vs 36 ng/ml)14. Therefore, an experimental model that enables sustained elevated plasma E2 would be especially valuable to mimic the steady-state elevation in E2 production during sexual maturation and further coordinate the gene expression data with changes in body weight and dynamics of muscle accretion.

Methods

Ethics statement

All animal experiments in this study were performed at the USDA/ARS National Center for Cool and Cold Water Aquaculture (NCCCWA) and approved by the NCCCWA Institutional Animal Care and Use Committee (protocol #50).

Experimental design

A total of forty fish weighing approximately 40 g were randomly assigned to four experimental tanks (n = 10 fish per tank). The study consisted of two treatments, including intraperitoneal injections of E2 and the delivery vehicle to serve as the control. Treatments were applied to fish in duplicate tanks. Estradiol was resuspended (10 μg/μL) in 95% ethanol and diluted to 2.5 μg/μL with vegetable oil. The control treatment contained an equal ratio of ethanol: vegetable oil as compared to E2 suspension. Estradiol and the vehicle injection methodology was adapted from previously published procedures used in tilapia63. Feed was withheld the day of E2 injection and throughout the study period. Fish were anesthetized with tricaine methanesulphonate (MS-222, 100 mg/l), weighed, and received intraperitoneal injections (2.0 μl/g body weight) of E2 (5.0 μg/g body weight) or the vehicle. Fish in one tank per treatment were each harvested 24 and 72 hours post-injection using lethal dose of MS-222 (300 mg/l). Skeletal muscle was dissected and immediately frozen using liquid nitrogen for further processing.

RNA extraction and quality control

Total RNA was extracted from skeletal muscle samples of fish treated with E2 and controls at 24 and 72 hours (6 biological replicates each) using TRIZOL reagent (Invitrogen, Carlsbad, CA) per the manufacturer’s suggested protocol. Quality and quantity of RNA was estimated using the A260:A280 ratio. Integrity and size distribution of RNA was evaluated using a Bioanalyzer 2100 (Agilent technologies, Santa Clara, CA). Four replicates of each treatment were subject to RNAseq.

cDNA library construction and sequencing

Library construction (poly(A)-enriched) and sequencing was completed by the Johns Hopkins University Genetic Resources Core facility High Throughput Sequencing Center (Baltimore, MD). cDNA libraries were constructed using TruSeq RNA library preparations and high throughput sequencing was completed using the Illumina HiSeq platform (100 bp single reads).

Sequence data processing, de novo assembly and differential expression analysis

Adaptor sequences were trimmed and ambiguous and low quality bases were removed, then read lengths less than 50 were removed. TRINITY17 was used to assemble all cleaned reads with default parameters. CD-HIT-EST was used to remove the shorter redundant transcripts when they were 100% covered by other transcripts with more than 99% identity18. All the cleaned reads were mapped to the assembled transcriptome by Bowtie20. RSEM was used to estimate and quantify gene expression levels from RNA-Seq data64. The final counts matrix file was used as input for the R package edgeR19 The exact test in edgeR was completed to discover the DEGs between the treatment groups. False discovery rate was used for multiple test correction. Any genes with a fold change greater than 2.0 and FDR of less than 0.05 were defined as DEG.

Validation of sequencing results

Sequencing results were validated using real time RT-PCR. RNA from 24 hour samples (6 replicates from each treatment) was used for cDNA synthesis using miScript II (Qiagen, Valencia, CA). Expression of mRNA and lncRNA was normalized to β-actin, which was not identified as a DEG by RNAseq analysis. Real time RT-PCR followed by melt curve analysis was performed to determine the specificity of the amplicons. Standard curves with 10-time serial dilutions of a pooled cDNA sample were generated to calculate the efficiency of qPCR. Five μl of diluted cDNA (1:1024) was used with iQ™ SYBR® Green Supermix (Bio-Rad) and 300 nM of forward and reverse primer in a final volume of 25 μl for the reaction. Cycle conditions were the same for all the primers except the annealing temperatures for different primers. Initial denaturation at 95 °C for 3 minutes followed by a denaturation at 95 °C for 40 sec, annealing for 30 sec, extension at 72 °C for 30 seconds for 40 cycles and a final extension for 10 min at 72 °C. Melt curve analysis was performed with an increase of 0.5 °C increase every cycle. Quantification cycle (Cq) values were used for quantification of expression using the log-linear equation of standard curve. Relative fold changes were calculated by setting the values of the controls to 1.0 and comparing the respective treatment groups. Statistical analysis was performed using Student’s t-test and those with a P-value < 0.05 were considered statistically significant.

LncRNA identification and GO and KEGG enrichment analysis of mRNAs

All DEGs were mapped to the rainbow trout genome65 using BLAT66. We used the pipeline reported in our previous study21 for lncRNA identification to detect differentially expressed lncRNAs. All differentially expressed mRNAs were subjected to similarity search against NCBI non-redundant (nr) protein database using BLASTx67 with an e-value cutoff of 1e-10. Gene names and GI were assigned to each mRNA based on the BLASTx result. ID mapping was performed using our in-house script to extract all associated GO terms for each mRNA. KEGG pathways were assigned to each mRNA using the online KEGG Automatic Annotation Server (KASS)68. The R package GOstats was used to run hypergeometric testing on GO and KEGG terms69. Redundant GO terms were removed by REVIGO, a Web server that summarizes long, unintelligible lists of GO terms by finding a representative subset of terms using semantic similarity measurement based clustering algorithm70.

Co-expression analysis

To identify co-expressed lncRNA-mRNA pairs, Pearson’s correlation coefficients were calculated based on the normalized expression value between every differentially expressed lncRNA and mRNA pair. Only the strong correlations (≥0.99) were selected to construct the network. The threshold of FDR was < 0.05. Cytoscape was used to construct the co-expression network71.

Data availability

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.