Crosstalk among lncRNAs, microRNAs and mRNAs in the muscle ‘degradome’ of rainbow trout

In fish, protein-coding and noncoding genes involved in muscle atrophy are not fully characterized. In this study, we characterized coding and noncoding genes involved in gonadogenesis-associated muscle atrophy, and investigated the potential functional interplay between these genes. Using RNA-Seq, we compared expression pattern of mRNAs, long noncoding RNAs (lncRNAs) and microRNAs of atrophying skeletal muscle from gravid females and control skeletal muscle from age-matched sterile individuals. A total of 852 mRNAs, 1,160 lncRNAs and 28 microRNAs were differentially expressed (DE) between the two groups. Muscle atrophy appears to be mediated by many genes encoding ubiquitin-proteasome system, autophagy related proteases, lysosomal proteases and transcription factors. Transcripts encoding atrogin-1 and mir-29 showed exceptional high expression in atrophying muscle, suggesting an important role in bulk muscle proteolysis. DE genes were co-localized in the genome with strong expression correlation, and they exhibited extensive ‘lncRNA-mRNA’, ‘lncRNA-microRNA’, ‘mRNA-microRNA’ and ‘lncRNA-protein’ physical interactions. DE genes exhibiting potential functional interactions comprised the highly correlated ‘lncRNA-mRNA-microRNA’ gene network described as ‘degradome’. This study pinpoints extensive coding and noncoding RNA interactions during muscle atrophy in fish, and provides valuable resources for future mechanistic studies.

factor 12 . MicroRNA, mir-182 indirectly regulates the expression of key muscle atrophy genes including atrogin-1, cathepsins and autophagy related genes by targeting transcription factor FOXO3 13 . Muscle specific microRNA, mir-1, regulates dexamethasone mediated muscle atrophy by targeting heat shock protein 70 (HSP70) 14 . However, microRNAs that regulate muscle atrophy in salmonids have not been previously investigated.
LncRNAs are a recently discovered class of noncoding RNAs with critical gene regulatory roles 15 . LncRNAs are known to regulate genes by direct interaction with microRNAs, mRNAs and proteins. Several lncRNAs bind microRNAs by sequence complementarity, and this interaction leads to cellular sequestration of microRNA (sponge effect) and lncRNA-mRNA competition for microRNA binding 16 . For example, lncRNA H19 binds and sponges away let-7 family microRNAs from repressing its protein-coding targets 17 . Similarly, muscle specific lncRNA, linc-MD1, competes with MAML1 and MEF2C to bind microRNAs mir-135 and mir-133, respectively 18 . LncRNA, MALAT1, modulates mir-133 mediated downregulation of serum response factor (SRF) by sharing mir-133 binding site 19 . LncRNAs can directly bind or physically interact with mRNA leading to mRNA decay 20 and translation suppression 21 . Some lncRNAs hybridize with the 3′ UTR of target mRNA and facilitate Staufen-1 mediated mRNA decay 20 . On the other hand, lincRNA-p21 directly binds JUNB and CTNNB1 mRNAs and suppresses their translation 21 . LncRNA's physical interaction with proteins modulates the stability 22 , cellular availability (sequestration) 23 , activity 24 and cellular localization 24 of proteins. For example, lncRNA, UPAT1, binds to UHRF1 protein and interferes with its ubiquitination and subsequent degradation 22 . LncRNA, MALAT1, binds to SR splicing factors and regulates their phosphorylation and hence cellular localization 24 . Both 'lncRNA-microRNA' and 'lncRNA-protein-coding genes' interactions are known to regulate development 18 , disease 25 and cancers 26,27 ; however, their involvement in skeletal muscle atrophy remains unknown.
To identify coding and noncoding genes involved in muscle atrophy associated with sexual maturation, we sequenced mRNAs, lncRNAs, and microRNAs from atrophying skeletal muscle of gravid fish and normal skeletal muscle of sterile fish. We subsequently performed differential gene expression of the two groups. In addition, we investigated functional interactions between DE lncRNAs, microRNAs and protein-coding genes in terms of expression correlation, genome co-localization and physical interaction to investigate gene-regulatory circuits during muscle atrophy. This study provides the first genome-wide lncRNA-mRNA-microRNA interaction network describing fish muscle degradation, and defining these interactions will clarify how energetic demand at sexual maturation triggers skeletal muscle atrophy.

Result and Discussion
Characteristics of atrophying and normal skeletal muscle. Differential expression of coding and noncoding genes was performed in atrophying skeletal muscle from diploid (2N), gravid fish in comparison to non-atrophying muscle of sterile triploid (3N) fish. This same set of skeletal muscle samples was used in our laboratory for several previous studies 1, 28,29 . Compared to sterile fish, fertile females yielded less separable muscle per whole body weight (49.9% ± 6.7% vs 62.6% ± 2.2%, p = 0.01), muscle protein (16.9% ± 0.7% vs ~19.1% ± 0.7%, p = 0.01) and muscle shear force (178 ± 19 gram/gram vs 240 ± 18 gram/gram, p = 0.01). On the other hand, atrophied skeletal muscle from 2N females had a higher moisture content (80.3% ± 0.7% vs 77.2% ± 0.6%) and pH (6.61 ± 0.03 vs 6.41 ± 0.04) (Fig. 1). Atrophied muscle also had numerically lower crude fat content than normal muscle, but the difference was not statistically significant (p = 0.30). These textural and compositional difference between two groups of muscle result from extensive muscle atrophy in gravid fish triggered by the energetic demand of sexual maturation. Differential expression of mRNAs, lncRNAs and microRNAs in atrophying muscle. To identify genes likely involved in skeletal muscle atrophy during sexual maturation, we performed deep lncRNA, mRNA and microRNA sequencing, and quantified DE genes between atrophying skeletal muscle of gravid fish and non-atrophying skeletal muscle from sterile fish. A total of 852 mRNAs, 1,160 lncRNAs and 28 microRNAs were DE between these two groups (FDR-p-value < 0.01, fold change: > 3 or <−3) (Fig. 2, Table 1 and Supplementary dataset 1). A total of 1,025 transcripts (352 mRNAs, 661 lncRNAs and 12 microRNAs) were upregulated and 1,015 transcripts (500 mRNAs, 499 lncRNAs and 16 microRNAs) were downregulated in atrophying muscle. Real time PCR validation of 4 transcripts from each DE list of lncRNAs, microRNAs and mRNAs is provided in Supplementary dataset 2. Previously, a microarray based approach performed on the same set of muscle samples identified only 82 upregulated and 120 downregulated protein-coding genes 1 , suggesting identification, in this study, of a large number of additional candidate genes involved in muscle atrophy. DE protein-coding genes, lncRNAs and microRNAs are described in separate sections hereafter.
Protein-coding genes. Many genes that promote proteolysis were significantly upregulated in atrophying skeletal muscle. At least 37 genes involved in protein ubiquitination, 22 genes involved in autophagy-related proteolysis, and 15 lysosomal and other proteases (cathepsin D, cathepsin B, cathepsin L and cathepsin Z) showed upregulation in atrophying muscle (Table 2 and Supplementary dataset 1). On the other hand, genes that negatively regulate the ubiquitin-proteasome system (ubiquitin carboxyl-terminal hydrolase 10, ubiquitin-like domain-containing CTD phosphatase 1 and uridine-cytidine kinase 2) and autophagy (CDGSH iron-sulfur domain-containing protein 2) were downregulated. Amino acid and fat biosynthetic genes were downregulated while genes involved in amino acid catabolism and transport were highly upregulated (Supplementary dataset 1). Similarly, genes associated with muscle sarcomere and extracellular matrix were downregulated, consistent with the loss of muscle mass and shear force during atrophy. As an example, 47 collagen-related genes and 24 non-collagen, extracellular matrix protein genes were significantly downregulated. A previous study also showed similar expression pattern of genes involved in protein ubiquitination and associated with the autophagy-lysosome system, extracellular matrix and sarcomere structure during muscle atrophy in mammals 30  or regulators, 28 were upregulated and 25 were downregulated (Supplementary dataset 1). While the proteolytic role for the majority of the TFs was unknown, some transcription regulators, such as zinc finger and BTB domain-containing protein 16 and ddb1-and cul4-associated factor 6, had known function in protein catabolism. On the other hand, development related TFs like myoD were downregulated. These findings suggest that muscle atrophy is triggered by upregulation of proteolytic and catabolic genes with concomitant downregulation of muscle sarcomere, extracellular matrix, muscle development and biosynthetic genes. The ubiquitin proteasome system appeared to be the major proteolytic system governing muscle atrophy. F-box only protein 32 (FBXO32) (atrogin-1), an E3 ubiquitin ligase, was the most highly upregulated genes in atrophying muscle suggesting that it might be the major player of muscle proteolysis during atrophy. Atrogin-1 genes transcripts, GSONMT00016768001 and GSONMT00031929001, exhibited 378-and 152-fold upregulation, respectively (Table 2). Their expression was validated by real time PCR (Supplementary dataset 2). Overexpression of atrogin-1 during starvation induced skeletal muscle atrophy has been reported previously in rainbow trout 5 , Atlantic salmon 6 and mammals 31 .
As fish progress from pre-spawning through spawning, severity of skeletal muscle atrophy increases as indicated by loss of muscle mass and muscle protein, and a reduction in muscle shear force, 28 as a measure of ultrastructural changes associated with muscle breakdown. To further investigate the potential contribution of DE genes in sexual maturation associated muscle atrophy, we looked at the trends of expression pattern of DE genes over 4 months during pre-spawning (July, November) and spawning (December and January) using RNA sequencing data from our previously described source 3 . Transcript abundance of ubiquitin-proteasome system genes and autophagy-related proteolytic genes remained constant in July and November, sharply increased in December, and then declined in January (Fig. 3a,b). Expression level of proteases positively correlated with severity of muscle atrophy in the aforementioned timeframes. Late December represents the time of peak sexual maturation associated muscle atrophy. Expression of genes coding for different atrogin-1 isoforms and cathepsin D was highest in December and then declined in January (Fig. 3c,d). Cathepsin D is involved in sexual maturation associated muscle atrophy, and we found an increase in Cathepsin D transcript level. Nonetheless, previously we did not observe a significant change in catalytic activity of cathepsin D during atrophy 28 . Extracellular matrix protein genes and development related genes showed opposite expression trends (Fig. 3e,f), consistent with the loss of muscle firmness and development during atrophy. These findings suggest that DE genes may serve as reliable candidate(s) critical to sexual maturation associated muscle atrophy in fish. However, the trends of the gene expression showed in Fig. 3 should be taken with caution because a single RNA-Sea library from 10 pooled fish at each time point was used and no replicates were available to run statistical analyses". Long noncoding RNAs (LncRNAs). For differential lncRNA expression, we initially performed genome-wide discovery of lncRNA transcripts using RNA-Seq reads sequenced from skeletal muscle of gravid and sterile fish. Approximately 15,000 lncRNA transcripts were identified from this assembly and merged with our previously published lncRNA reference 32 ; and both sources were used as a reference for gene expression analysis. A total of 1,160 lncRNAs were DE between atrophying and non-atrophying skeletal muscle. Of 1,160 DE lncRNAs, 225 and 10 lncRNAs had sequence homology with lncRNAs from Atlantic salmon and zebrafish, respectively (sequence identity: >80%, E value: <E-10, query cover: >50 nucleotides) (Supplementary dataset 3), but their functional annotation was not available in any species. Like protein-coding mRNAs, expression level of DE lncRNAs correlated with the severity of muscle atrophy during pre-spawning and spawning. Transcript abundance of upregulated lncRNAs remained constant during pre-spawning, but drastically increased in December and then declined in January (Fig. 3g). On the other hand, transcript abundance of downregulated lncRNAs showed an opposite trend (Fig. 3h). These findings suggest that expression of these DE lncRNAs may be involved in sexual maturation associated muscle atrophy in rainbow trout.

MicroRNAs.
A total of 28 microRNAs were DE between skeletal muscle of gravid and sterile fish. Of them, differential expression of mir-1, mir-133, and mir-29 during mammalian muscle atrophy has been previously reported 11,33 , but the remainder of DE microRNAs was reported for the first time (Table 1). A total of 665 unique mRNA genes were predicted as potential target genes of these 28 DE microRNAs; of these mRNA genes, 17 were also DE. Some of these DE microRNAs and their predicted DE mRNA targets showed reciprocal differential expression (Supplementary dataset 4). As an example, mir-29a predicted target mRNAs encoding collagen alpha, ATP binding cassette subfamily f member 3, alanine tRNA synthase, scavenger receptor class b member 1, and fk506-binding protein 2 were downregulated while mir-29a was upregulated in atrophying muscle. Similarly, mir-125b-1 was downregulated, and its predicted mRNA targets encoding CCAAT enhancer-binding protein delta and pancreatic progenitor cell differentiation and proliferation factor A were upregulated. Out of the DE microRNAs, sixteen were downregulated microRNAs, potentially targeting a total of 206 different protein-coding genes. Twenty-six of the predicted target genes were proteolytic enzymes, including genes involved in ubiquitin-proteasome and autophagy-lysosome mediated proteolysis (Supplementary dataset 4). Twelve upregulated microRNAs were predicted to target 468 different protein-coding genes. Consistent with  Table 1. DE microRNAs between atrophying muscle of gravid fish and normal skeletal muscle of sterile fish. Positive and negative value of fold change represent upregulation and downregulation respectively in atrophying skeletal muscle of gravid fish. Fold change was considered significant at cutoff: > 3 or <−3, FDRp-value < 0.01. Several isoforms of let-7 were downregulated and several isoforms of mir-29 were upregulated. their upregulated expression during atrophy, 101 predicted target genes were directly involved in extracellular matrix, muscle structure, or development (Supplementary dataset 4). Aforementioned findings suggest that some genes involved in muscle atrophy may not be necessarily regulated at transcription level, and its fate is determined post-transcriptionally by regulated expression of microRNA. Let-7j was the most highly downregulated microRNA (−1056 × fold) in atrophying muscle (Table 1). It was predicted to target 63 different protein-coding genes that account for a wide range of functions (Supplementary dataset 4). Consistent with its downregulation in atrophying muscle, some of its predicted targets were proteolytic genes such as E3 ubiquitin-protein ligase NRDP1, ubiquitin-conjugating enzyme E2 and protein VPRBP (Supplementary dataset 4). Conversely, six different isoforms of mir-29 were highly upregulated in atrophying muscle (Table 1). Mir-29a, the most highly upregulated microRNA in atrophying muscle (11.7 × fold) was predicted to target 78 genes; the highest number of predicted target genes among DE microRNAs (Supplementary dataset 4). Consistent with its upregulation in atrophying muscle, predicted target genes of mir-29a included genes involved in muscle differentiation (IGF-BP 5), muscle sarcomere (e. g. myosin) structure, extracellular matrix (e. g. collagen), fat biosynthesis (e. g. long-chain-fatty-acid-ligase acsbg2 and acyl-coenzyme a thioesterase 11), protein synthesis (e. g. 60 s ribosomal protein l7) and development (e. g. prospero homeobox protein 1). These findings suggest that DE microRNAs may contribute to muscle atrophy by regulating proteolysis and other genes during muscle atrophy.
Tissue specific and temporal expression of DE lncRNAs and mRNAs. LncRNAs show strict spatial (tissue specific) and temporal (time dependent) expression patterns 34 . To investigate tissue specificity of DE genes, we studied their expression pattern across 13 vital tissues including red and white muscle (see method section for classification of tissue specific genes). About 40% (462/1,160) of DE lncRNAs and approximately 41% (348/852) of DE mRNAs were 'specific' to red or white muscle (Fig. 4). These specificities indicated about 2.5-fold enrichment of muscle specific lncRNAs, and about 5.5-fold enrichment of muscle specific mRNAs in  Table 2. Selected proteolytic genes highly upregulated in atrophying skeletal muscle of gravid female rainbow trout relative to non-atrophying muscle of same-aged sterile rainbow trout. Fold change was considered significant at cutoff: > 3 or <−3, FDR-p-value < 0.01. Interestingly, a majority of the most highly upregulated mRNAs were muscle 'specific' . As an example, 47 out of 61 mRNAs, with fold change greater than 15, had muscle restricted expression patterns that included mRNA encoding atrogin-1. Muscle specific expression of atrogin-1 has also been previously reported in mammals 31 . In addition, some muscle specific lncRNAs, such as linc-MD1, are known to play an important role in regulation of muscle specific genes 18 . In addition to tissue specific expression, ~38% (442/1,160) of DE lncRNAs and ~40% (342/852) of DE mRNAs exhibited temporal (time dependent) expression patterns during pre-spawning and spawning (Fig. 4). These findings suggest that a significant proportion of the DE transcriptome in atrophying muscle is comprised of muscle specific gene expression, exhibiting a responsive expression pattern of muscle atrophy during sexual maturation.
Genomic co-localization of DE lncRNA and mRNA genes. Functionally related lncRNAs and mRNAs physically co-localize in the genome 35 Supplementary dataset 6, tab 3). The difference was statistically significant (Chi square p-value < 0.001) suggesting that co-localized lncRNA and protein-coding genes tend to be more

DE lncRNAs acting as microRNA sponges or microRNA precursors. Direct 'lncRNA-microRNA'
binding has important functional consequences including lncRNA mediated sponging of microRNA and 'lncRNA-microRNA' competition for mRNA binding 16 . To identify DE lncRNAs that potentially interfere with microRNA mediated gene regulation, we searched for high confidence microRNA binding sites in DE lncR-NAs and mRNAs. A total of 134 trout microRNAs had binding sites in DE lncRNAs as well as DE mRNAs (Supplementary dataset 7). Some of these microRNAs, such as mir-133, mir-214, and mir-221, are known to regulate muscle atrophy or proteolysis 38,39 . DE lncRNAs shared microRNA binding sites with mRNAs encoding important proteolytic proteins such as atrogin-1, cathepsins, serine proteases, and several enzymes in the ubiquitin proteasome system (Table 3 and Supplementary dataset 7). For example, atrogin-1 (GSONMT00016768001), the most highly upregulated gene in atrophying muscle, shared mir-22-3p binding sites with Omy200063021 and Omy400145202. Expression patterns of a subset of LncRNAs and mRNAs that shared microRNA binding sites and were strongly correlated is shown in Table 3. In addition to acting as a microRNA sponge, some lncRNAs serve as precursors of microRNAs and other classes of small noncoding RNA (sRNAs) 34,40 . DE lncRNA, Omy400148395, harbored mir-27 loci, and Omy200105075 harbored aly-miR-398c-like loci. Expression of mir-27 and aly-miR-398c-like microRNAs was positively correlated with their potential host lncRNAs; correlation R values were 0.60 and 0.81, respectively. This observation suggests that these microRNAs could be generated by post-transcription processing of lncRNA transcripts. Together, these findings suggest that some DE lncRNAs may sequester or generate microRNAs involved in muscle atrophy.

Physical interaction between lncRNAs and protein-coding genes.
Direct 'lncRNA-mRNA' physical interactions lead to mRNA decay 20 and translation suppression 21 . To investigate potential existence of lncRNA-mRNA physical interactions, we used the IntaRNA tool that considers site accessibility and user defined seed requirement to predict the interaction 41 . At an interaction energy threshold <−100 Kcal/mole, 1,151 DE lncRNA-DE mRNA pairs exhibited potential physical interactions (Table 4 and Supplementary dataset 8). Interestingly, lncRNA-mRNA pairs, showing evidence of physical interaction, were more strongly correlated in expression (R > 0.85) than lncRNA-mRNA pairs without the potential physical interaction (Table 4 and Supplementary dataset 8).
In addition to interacting with mRNAs, several lncRNAs showed evidence of physical interaction with proteins of DE mRNAs. 'LncRNA-protein' physical interaction was computed using CatRapid Omics tool 42 . A total of 14,602 'DE lncRNA-protein' pairs showed evidence of physical interaction at an interaction strength ≥96% and a discriminative power ≥96% (Table 4 and Supplementary dataset 8). Interestingly, 'DE lncRNA-DE protein' pairs showing evidence of direct physical interactions were more strongly correlated in expression (at transcript level) (R > 0.85) than 'DE lncRNA-DE protein-coding gene' pairs without the evidence of physical interaction (Chi square p-value < 0.001). Atrogin-1, cathepsin, and several enzymes of the ubiquitin-proteasome system were among the proteins with potential interaction with lncRNAs. 'LncRNA-protein' binding regulates protein's stability 22 , availability 23 , activity 24 and cellular localization 24 , suggesting that DE lncRNAs may have an important role in determining fate of DE protein-coding genes during muscle atrophy.
'LncRNA-mRNA-microRNA' interactome in atrophying muscle. LncRNAs, microRNAs, and mRNAs comprise interacted gene regulatory networks in the cell 43 , probably due to mutual regulation between microRNA and lncRNA 16 . To investigate existence of such gene interaction during muscle atrophy, we computed lncRNA-mRNA-microRNA interaction networks based on their expression pattern across 30 RNA-Seq datasets. At a correlation threshold R > 0.97 or <−0.97, about 50% (1,584) of DE transcripts were components of strongly correlated gene networks (Fig. 5). Interestingly, a majority of the correlated transcripts clustered in one of two major networks; the first comprised of downregulated transcripts and the second comprised of upregulated transcripts. The first network consisted of 430 transcripts (137 lncRNAs, 219 mRNAs and 74 microRNAs). The second   network consisted of 960 transcripts (559 lncRNAs, 235 mRNAs and 166 microRNAs). The second network appeared to be an interacted gene "reactome" regulating muscle proteolysis because almost all upregulated proteolytic genes, including enzymes of ubiquitin proteasome system, autophagy related proteolytic genes and other proteases such as cathepsins, were in this network. Similarly, a majority of the upregulated microRNAs including mir-29, let-7, let-7d and mir-132 were in the network. Interestingly, atrogin-1 transcripts (the most highly upregulated transcripts) and mir-29a (the most highly upregulated microRNA) were in the center of this network suggesting a key role in the network. A crucial role of atrogin-1 in muscle atrophy has repeatedly been reported in fish and mammals 5,6,31 . We identified this network as 'the rainbow trout muscle degradome' because it appeared to be a regulatory network of muscle degrading coding and noncoding genes. The network was largely comprised of the genes that showed evidence of physical interaction with each other. A sub-network of the degradome shown in Fig. 5 consist of atrogin-1 transcripts, all DE lncRNAs that bind to atrogin-1 mRNAs, and microRNAs that either bind to lncRNAs and/or atrogin-1 mRNAs. The aforementioned observations suggest that coding and noncoding genes involved in muscle atrophy work in the form of a highly interacted gene network.

Conclusion
Sexual maturation associated skeletal muscle atrophy serves as an excellent model to study piscine muscle proteolysis 1,3,29 . Previous efforts to investigate fish muscle proteolysis have provided limited information because these studies relied on individual or a limited set of protein-coding genes 1,5,6 . In the present study, we used deep lncRNA, mRNA, and microRNA sequencing approaches to investigate genes and gene regulatory networks that regulate muscle proteolysis in fish. Through investigation of the atrophying muscle transcriptome, we elucidated that fish muscle atrophy, like mammalian muscle atrophy, is regulated mainly by the ubiquitin-proteasome system. In addition, many autophagy-lysosomal proteases and transcription factors appeared to participate in muscle proteolysis during atrophication. Atrophying muscle exhibited upregulation of proteolytic genes with concomitant downregulation of genes involved in muscle sarcomere, extracellular matrix, protein and fat biosynthesis, and development. This trend in expression pattern of genes in atrophying muscle correlated well with characterizations of the atrophying muscle phenotype (e. g. muscle mass, protein content and muscle shear force), suggesting essential roles for DE genes in muscle atrophy. The present study identified a large number of new candidate coding and noncoding genes in addition to the genes identified by previous microarray and proteomic approaches 1,29 suggesting that the RNA-Seq approach has identified a large number of reliable candidate genes involved in muscle proteolysis/degradation. In present study, we characterized lncRNAs potentially involved in fish muscle proteolysis, and we investigated lncRNA-mRNA, lncRNA-microRNA and mRNA-microRNA interactions that potentially regulate muscle atrophy. A majority of DE lncRNAs and DE mRNA genes were co-localized in the genome and correlated in expression. DE lncRNAs appeared to physically interact extensively with DE protein-coding genes at transcript  and protein levels. Similarly, DE lncRNA also showed potential to bind and sequester cellular microRNAs implicated in muscle proteolysis. LncRNA, mRNA and microRNAs that exhibited the aforementioned interactions are components of a highly correlated gene network in atrophying muscle. This finding indicates that the majority of genes involved in muscle proteolysis are expressed simultaneously by a common gene transcription program. Important to mention, 'DE lncRNA-DE protein-coding' gene pairs that either co-localized in the genome or showed evidence of direct physical interaction or competed for a common microRNA binding, were more frequently correlated in expression than random 'DE lncRNA-DE protein-coding' gene pairs. Perhaps this work is the first genome wide study that provides links between expression correlation and potential functional interactions between lncRNA and mRNA in fish. The present study has investigated potential coding and noncoding RNA interactions during muscle atrophy and contributes to our understanding of how energetic demand of sexual maturation triggers skeletal muscle atrophy in fish.

Materials and Methods
Ethics statement. Fish muscle tissues were obtained commercially from a private farm; therefore, Institutional Animal Care and Use Committee (IACUC) approval was not required.
Fish population and muscle sampling. We previously described the fish population in another study 1 .
Briefly, mature sterile (3N: triploid) and fertile (2N: diploid) female rainbow trout (about 500 gram) were obtained from Flowing Springs Trout Farm (Delray, WV) during spawning season; these fish were cultured in identical raceways. Water from a common spring was circulated in raceways at temperature 13 ± 3 °C. Both groups of fish were fed ad libitum (Zeiglar Gold; Zeigler Bros., Gardeners, PA) via demand feeder until sampling. At the time of muscle sampling gonado-somatic index (GSI) of fertile fish was 15.8 ± 0.3 (n = 5) compared to 0.3 ± 0.2 (n = 5) in sterile fish confirming the gravid stage of fertile fish. White muscle tissue of 8 fish (4 fertile and 4 sterile) was collected from the dorsal musculature, flash frozen in liquid nitrogen and stored at −80 °C until RNA extraction. Total RNA was extracted from the muscle using TRIzol method (Invitrogen, Carlsbad, CA). For phenotype measurement, boneless and skinless muscle fillet was obtained in a manual filleting procedure. Muscle yield was measured as percentage of whole body weight. A representative portion of the muscle fillet was used for shear force, pH, and proximate analyses (e.g. crude protein content, crude fat content, and moisture content).

Figure 5.
Gene expression network of DE lncRNAs (blue node), DE mRNAs (green node) and microRNAs (pink node) (R > 0.97 or <−0.97). Note that most of the DE genes are clustered in one of the two major networks. The larger network (degradome) comprises of upregulated genes and smaller network comprises of downregulated genes. In the network of upregulated transcripts, 4 atrogin-1 transcripts (red nodes) and mir-29 isoforms (black nodes) are in the center of the network. Sub-network drawn from the larger network contains 4 atrogin-1 transcripts, lncRNAs that bind to atrogin-1, and microRNAs that either bind to the lncRNAs and/or atrogin-1. Note: edges that connect nodes (genes) represent correlated expression at R cutoff 0.97 > or <−0.97; the shorter the length, the stronger the expression correlation.

Discovery of lncRNAs.
LncRNAs from sequencing reads were identified by using the pipeline we described previously 32 . Briefly, reads were mapped to rainbow trout reference genome 44 and assembled using TopHat and Cufflinks, respectively. Transcripts shorter than 200 nucleotides were filtered out. Protein-coding transcripts were removed by their sequence homology with NCBI protein entries. In addition, Coding Potential Calculator (CPC) 45 tool was used to remove any transcripts with protein-coding potential (index value < −0.5).  47 was used to quantify microRNA in Bio Rad CFX96 ™ System. The endogenous controls used for normalization were B-actin for mRNA and lncRNAs, and U6 for microRNA. None of the endogenous control genes was differentially expressed in this study. Fold changes in gene expression was calculated by using ΔΔCt method as described previously 35,48 . Mann-Whitney U test was used to check if the transcript level between atrophying and control muscle was statistically significant (p < 0.05). All 12 transcripts subjected to qPCR validation had Mann Whitney U test p-value < 0.05.
Identification of tissue specific genes. The expression pattern of DE genes was investigated across 13 vital tissues: red muscle, white muscle, spleen, liver, skin, testis, brain, intestine, stomach, kidney, head kidney, gill and fat 44,46,49 . To identify tissue specific genes, we used a statistical approach described by 50 . The normalized expression value (z score) of every gene was calculated in each tissue from TPM counts. Each gene was classified as 'specific' to a tissue if z score was greater than 1.5 in that tissue and less than 1.5 in remaining 12 tissues which is similar to a previous study 50 . We used the same approach to identify genes 'specific' to a particular month during pre-spawning and spawning 3 , as described in the results section.
LncRNA, mRNA and microRNA co-expression. For mRNA-ncRNA correlation analysis, 8 RNA-Seq samples from gravid and sterile fish and 22 RNA-Seq samples sequenced from muscle of selectively bred fish families developed at USDA/NCCCWA were used. Detailed description of these 22 samples, fish population, and sampling procedure has been previously described 51 . TPM values (for lncRNAs and mRNAs) and the total count (for microRNAs) were calculated in all 30 samples (8 samples from 2 ploidy groups and 22 samples from fish families) by mapping reads to corresponding mRNA, lncRNA, and microRNA references as described above. TPM and total count values were normalized by using a scaling method as previously described 52 , and normalized values were used to calculate expression correlation coefficients using Pearson correlations. A lncRNA-mRNA-microRNA expression network was constructed using Expression Correlation in cytoscape 53 .

Identification of cis regulatory promoter motifs.
Promoter regions of all neighboring/overlapping lncRNA-mRNA genes were scanned for putative Transcription factor (TF) binding motifs using ALGGEN PROMO TF motif search tool 54 . Maximum dissimilarity rate was set to 5%, and RE equality/query was set to <0.05. Identification of microRNA-harboring lncRNAs. The rainbow trout genome reference 44 was annotated with rainbow trout lncRNA reference sequences mentioned above. Rainbow trout pre-microRNA sequences (~64 nts long) from a recently published source 55 were aligned to the lncRNA-annotated genome assembly. Pre-microRNA sequences that perfectly align (with no mismatch or gap) within annotated lncRNA loci in the genome were reported.
LncRNA and mRNA targets of microRNAs. In the case of mRNAs, microRNA binding sites were searched in 3′ UTR; whereas for lncRNAs, target sites were searched for throughout the entire sequence length. Three target prediction algorithms: miRanda, PITA and TargetSpy were used to find target genes using sRNAtoolbox 56 . If the same target site is predicted by all 3 tools, it was considered as a potential microRNA target site. For all tools, minimum energy threshold was chosen as −20 Kcal/mole. Threshold scores chosen were 150 for miRanda and 0.99 for TargetSpy.
Prediction of lncRNA-mRNA and lncRNA-protein interaction. LncRNA-mRNA direct physical interaction was predicted using IntaRNA-RNA-RNA interaction tool 41 . All lncRNA-mRNA interactions were recorded at an interaction energy threshold < −100 Kcal/mole. LncRNA-protein interactions were predicted using CatRapid Omics tool 42 . Interaction strength and discriminative power (a measure of predictability of interaction) were set at a ≥96% to consider lncRNA-protein interaction as putative interaction.