RNA-Seq reveals MicroRNA expression signature and genetic polymorphism associated with growth and muscle quality traits in rainbow trout

The role of microRNA expression and genetic variation in microRNA-binding sites of target genes on growth and muscle quality traits is poorly characterized. We used RNA-Seq approach to investigate their importance on 5 growth and muscle quality traits: whole body weight (WBW), muscle yield, muscle crude-fat content, muscle shear force and whiteness. Phenotypic data were collected from 471 fish, representing 98 families (~5 fish/family) from a growth-selected line. Muscle microRNAs and mRNAs were sequenced from 22 families showing divergent phenotypes. Ninety microRNAs showed differential expression between families with divergent phenotypes, and their expression was strongly associated with variation in phenotypes. A total of 204 single nucleotide polymorphisms (SNPs) present in 3′ UTR of target genes either destroyed or created novel illegitimate microRNA target sites; of them, 78 SNPs explained significant variation in the aforementioned 5 muscle traits. Majority of the phenotype-associated SNPs were present in microRNA-binding sites of genes involved in energy metabolism and muscle structure. These findings suggest that variation in microRNA expression and/or sequence variation in microRNA binding sites in target genes play an important role in mediating differences in fish growth and muscle quality phenotypes.

in fish. As an example, mir-143 and mir-203b target myoD, a member of myogenic regulatory factors (MRFs), in Mandarin fish and Nile tilapia, respectively 12,13 . A novel zebrafish microRNA, mir-In300, targets dickkopf-3 (dkk3) gene and hence abolishes the promoter activity of myogenic protein 5 (myf5) 14 . Similarly, zebrafish mir-214 controls hedgehog signaling mediated speciation of muscle cell by regulating suppressor of fused (sufu) gene 15 . Further, Let-7, mir-19 and mir-130 show regulated expression during transition from muscular hyperplasia to hypertrophy in fish 16 . Despite these aforementioned studies, there is still need for complete microRNAome expression and genetic variant profiles to understand the genetic basis of variation in muscle growth and quality traits in food fish production.
Scientists at the USDA, ARS National Center for Cool and Cold Water Aquaculture (NCCCWA) have developed a pedigreed line of rainbow trout selectively bred for 5 generations for improved growth performance to the standard US market body weight and beyond 17 and have characterized muscle yield and quality traits in nucleus families within the line. The objectives of this study were to 1) investigate association of microRNA expression with muscle growth and muscle quality traits in the NCCCWA growth-selected line and 2) investigate effects of single nucleotide polymorphisms (SNPs) in microRNA binding sites on growth and muscle quality traits. Using high throughput deep small RNA sequencing approach, we identified differentially expressed (DE) microRNAs between fish families showing contrasting phenotype for whole-body weight (WBW), muscle yield, crude fat content, shear-force and whiteness of the muscle. We performed 'phenotype-microRNA expression' association using a random sample of 90 fish from 3 rd -generation families (2010 hatch year) to investigate potential impact of microRNA expression to the phenotype in the population. SNPs capable of creating or destroying microRNA binding sites in protein coding target genes were identified and their functional consequence on growth and muscle quality phenotypes was evaluated by 'SNP-phenotype' association analysis using a sample of 786 fish from 3 rd -and 4 th -generation families (2010 and 2012 hatch years).

Result and Discussion
Muscle trait phenotypes and Small RNA sequencing. Phenotypic divergence between the 4 highest-ranked and 4 lowest-ranked families for the 5 traits were: WBW (1221.6 g ± 84.25 vs. 502.1 ± 28.0 g); muscle yield (50.9% ± 1.6 vs. 43 1). Since some fish families were common between the traits, the total number of sequenced families was 22 (Table 1). High throughput small RNA sequencing resulted into mean sequencing depth of 14.5 million reads per sample. After trimming of sequencing adaptor, average length of reads was 22 nucleotides (Supplementary Dataset 1A), which is a typical length average for mature microRNAs. After filtration and adapter trimming, the average number of reads in each sample was 11.9 million. On average, ~0.2 million potential microRNA transcripts were detected from trimmed reads in each sample. Of these potential microRNA transcripts, ~17.5% had sequence homology with mature microRNAs in the miRBase database. From these annotated microRNAs in each sample, different variants of the same microRNA were collectively counted as a single microRNA, which resulted into an average of 1,154 unique microRNA per sample (Table 1).
RNA-Seq based principal component analysis (PCA) of 22 families sequenced for gene expression study showed no cluster suggesting genetic divergence within the population (Supplementary Dataset 1B).
Differentially expressed microRNAs between high and low ranked families for growth and muscle traits. In order to identify microRNAs associated with growth and muscle quality traits, we profiled DE microRNAs between families with divergent phenotypes. None of the microRNAs were DE in association with WBW. However, 90 microRNAs were DE in remaining four muscle traits ( Fig. 2 and Supplementary Dataset 1C). Two of the DE microRNAs were recently discovered microRNAs in rainbow trout; new-miR-66 and new-miR-34 18 . A total of 18, 9, 56 and 77 microRNAs were DE between high and low ranked families for muscle yield, crude fat content, shear force and muscle whiteness, respectively. Most of the DE microRNAs were shared between several phenotypic traits ( Fig. 2   Phenotypic difference for whole body weight (WBW) and 4 muscle quality traits (muscle yield, crude fat content, shear force and whiteness index) of top 4 high ranked and 4 low ranked families (5 fish/family) of selectively-bred trout at ~13-month post-hatch. Phenotypes were statistically different between high ranked and low ranked families (P < 0.01). Error bars on graph show standard deviation.
in association with muscle whiteness. This observation suggests that common mechanisms including coordinated expression of microRNAs may affect muscle yield and quality traits in rainbow trout. This argument is supported by previous report pinpointing to interrelated nature of these muscle traits in fish 19 . Heat map of fold change of differentially expressed (DE) microRNAs between high vs low ranked families of various traits (left) and Venn-diagram showing shared DE microRNAs between different traits (right). In heat map, dark green and yellow colors indicate downregulation and upregulation respectively in high ranked family compared to low ranked family. Dark color indicates no differential expression of microRNAs. Note that the value of color limit (−5 to 5) does not reflect true fold change as values of fold change were transformed (log2) and color scale was adjusted to make the heat map more visible. True fold change is given in Supplementary Dataset 1B.
Interestingly, the direction of change of the shared DE microRNAs were consistently either positively or negatively correlated between traits. For example, all 18 DE microRNA that were downregulated in association with increased muscle yield were also downregulated in families showing high shear force (R 2 = 0.54) (Fig. 2 [heat map] and Supplementary Dataset 1C). In contrast, most of these downregulated microRNAs in families with high muscle yield (16 out 18) were upregulated in the high-whiteness families (R 2 = −0.56). In addition, 47 out of 49 of the shared DE microRNAs between the shear force and muscle whiteness groups showed opposite pattern of differential expression (R 2 = −0.87). Similarly, all 5 shared DE microRNAs between crude fat and whiteness groups were downregulated in high ranked families of both traits (R 2 = 0.87). Out of 3 DE shared microRNAs between shear and fat group, 2 microRNAs showed opposite fold change pattern between the traits. In accordance with this observation, growth and muscle quality phenotypes of 500 fish population used in this study showed correlation between traits ( Table 2). WBW showed positive correlation with muscle yield (R 2 = 0.32, p < 0.0001) and crude fat content (R 2 = 0.33, p < 0.0001). Similarly, muscle yield showed very weak but significant positive correlation with crude fat (R 2 = 0.061, p < 0.0001) and shear force (R 2 = 0.03, p = 0.0003), and negative correlation with whiteness (R 2 = −0.023, p = 0.0009). Crude fat content and whiteness had weak but positive correlation ( Table 2). This finding suggests that correlation among phenotypic traits could be, at least partially, explained by variation in expression level of DE microRNAs. Consistent with this observation, a recent report in salmon has indicated that crude fat content is negatively correlated with muscle shear force and positively correlated with L* (lightness) and b* (yellowness) of raw salmon muscle 19 .

Targets of DE microRNAs and their functional annotation.
A total of 6,837 protein-coding genes were identified as potential targets with high confidence for 90 DE microRNAs (Supplementary Dataset 2). To investigate the functional significance of the predicted target genes, we performed gene enrichment analysis (GEA) using DAVID 20,21 . A list of all overrepresented pathways in different GO categories are provided in Supplementary Dataset 3. As shown in Fig. 3, in the biological process category, genes involved in multicellular organism development (4.2-fold), transcription (3.7-fold) and mitosis/cell division (4.0-fold) were highly overrepresented. Other significantly enriched pathways included RNA processing, DNA repair, gene silencing by RNA, protein ubiquitination, lipid metabolism, muscle organ development, regulation of growth, cell proliferation and apoptosis. Similarly, among the signal transduction pathways, Wnt signaling pathway was overrepresented, which is involved in skeletal muscle growth/myogenesis [22][23][24] and regulation of growth control genes 25 .
In the molecular function pathways category, most of the overrepresented gene sets had functions related to epigenetic gene regulation such as methyltransferase activity, histone binding and chromatin binding activity (Fig. 3). Methyltransferases are important in epigenetic gene-regulation and are known to regulate skeletal muscle growth by regulating expression of myoD transcription factor 26 . Another enriched gene set included various ligases that post translationally modify proteins by adding specific amino acid; e.g. protein-glycine ligases, protein-glutamic acid ligase and tubulin-glycine ligase. These findings suggest that DE microRNAs may control growth and muscle quality phenotypes primarily via post-transcriptional regulation of genes involved in development, muscle growth as well as the genes involved in epigenetic gene regulation and protein modification.
MicroRNAs associated with growth and muscle quality traits. To estimate the association between microRNA expression and variation in phenotypes, 12 highly and commonly (among traits) DE microRNAs were selected for 'microRNA expression-phenotype' regression analysis. Expression level of microRNA was qPCR-quantified and correlated to phenotypes in 90 randomly selected fish from the same population. Out of the 12 microRNAs, 10, 12, 10, 5 and 6 correlated with WBW, muscle yield, crude fat content, shear force and muscle whiteness respectively (cut off: R 2 > 0.05, p-value < 0.05) ( Table 3). When the effect of all 12 microRNA expression levels was combined, about 31%, 42%, 22%, 13% and 26% variation in WBW, muscle yield, muscle crude fat, muscle shear force and muscle whiteness, respectively, was explained by variation in microRNA expression.
MicroRNAs associated with each phenotypes and their relevant protein coding target genes are described in separate section below.
Whole body weight (WBW). No microRNA was DE between the high and low ranked families of WBW phenotype. However, microRNA expression-phenotype regression analysis performed on 90 randomly sampled fish showed that variation in expression of 10 DE microRNAs explained significant variation in WBW (Table 3). These microRNAs included muscle specific myogenic microRNAs as well as other microRNAs that were DE in response to muscle yield and/or muscle crude fat content. Expression of individual microRNA explained 5.6% to 16.5% variation in WBW depending on the microRNA (Table 3). We calculated the fold change between high vs low ranked WBW families for these 10 microRNAs based on qPCR data. Majority of these microRNAs showed small but significant fold change between high vs low ranked WBW families. It is possible that RNA-Seq failed to detect such a small change in fold change, which qPCR did. Future mechanistic study involving gene knockout or dose-dependent response of individual microRNA on WBW may help understand biology of WBW growth. In addition, all 9 DE microRNAs that were chosen for phenotype-genotype regression analysis, negatively correlated with muscle yield (p < 0.05). Variation in expression level of each microRNA explained 6.8% to 18.8% of the variation in muscle yield, with the most significant microRNA being mir-1a-3p (Table 3). Although expression of mir-10c-5p and mir-181a-3p was not statistically different in the RNA-Seq data between high and low ranked muscle yield groups, their expression was significantly and negatively correlated with muscle yield based on qPCR quantification in random 90 individuals (Table 3). To the best of our knowledge, 8 out of 18 DE microRNAs (mir-19d, mir-100-5p, mir-99a-5p, mir-148-3p, mir-199-5p, mir-148-3p*, mir-19b-3p and mir-140) were not reported before as associated with muscle growth. The remaining 10 DE microRNAs are known to directly or indirectly regulate skeletal muscle development in different species by regulating genes involved in myogenesis 12,27 , muscle-growth-related signal transduction pathways 9, 28, 29 and cell cycle 30 . As an example, mir-1 and mir-133 are responsible for more than half of the microRNA-mediated gene regulation in muscle of zebrafish 11 . Similarly, mir-143 regulates expression of myoD in skeletal muscle of fish 12 . These findings suggest that differential gene expression has identified several microRNAs previously known to regulate myogenesis as well as several additional microRNAs potentially implicated in muscle growth. A total of 1,743 protein-coding genes were predicted as potential targets of 18 downregulated microRNAs in high muscle yield group (Supplementary Dataset 2). Interestingly, 10 significantly enriched gene pathways in the target gene list were directly involved in muscle growth e. g. muscle organ development (4.4-fold), skeletal muscle tissue development (4.5-fold) and muscle cell differentiation (3.1-fold) (FDR-p < 0.05) (Supplementary Dataset 4). Other enriched gene pathways included chromatin modification (5.4-fold), transcription (3.2-fold), cell cycle (4.7-fold), and multicellular organism development (4.6fold). These findings suggest that DE microR-NAs may regulate muscle yield by regulating genes directly involved in skeletal muscle development, transcription, cell cycle, and/or protein degradation 12,27,29,30 .
Three of the nine microRNAs (mir-10c-5p, mir-10b-5p and mir-181a-3p) were qPCR analyzed for association with variation in muscle fat content in 90 individual fish. The regression analysis showed that these microRNAs explained 6.3% to 8.1% of the variation in fat content in muscle (Table 3). In addition to these 3 DE microRNAs, 7 non-DE microRNAs that were associated with muscle yield (mir-1a-3p, mir-19b-3p, mir-148-3p*, mir-20-5p, mir-133a-5p, mir-143-5p and mir-99b) were also significantly correlated with muscle crude fat content (FDR-p < 0.05), (Table 3). These findings suggest that DE microRNAs may play crucial role in post-transcriptional regulation of genes associated with crude fat content in muscle.  Muscle Shear force. A total of 56 microRNAs were DE between high and low shear force groups, of them 46 microRNA were downregulated and 10 microRNAs were upregulated in fish with high shear force ( Fig. 2 and Supplementary Dataset 1C). Correlation between expression level and variation in muscle shear force was studied for 10 DE microRNAs in 90 fish, of them, 5 microRNAs (mir-126-3p, mir-148-3p*, mir-20-5p, mir-206 and mir-133a-5p) explained 6.5% to 9.0% variation in muscle shear force depending on the microRNA (FDR-p < 0.05) ( Table 3). To the best of our knowledge, association of microRNAs with muscle firmness has not been previously reported in salmonids. The connective tissue proteins, especially collagen, play an important role in determining muscle shear force. Density of collagen fibers in connective tissue determines fish muscle firmness 36 . Interestingly, out of 56 DE microRNA, 31 microRNAs targeted 49 genes coding for collagen proteins or regulators of collagen biosynthesis, metabolism or structure (Supplementary Dataset 2, see name and GO terms of target genes). In addition, 31 DE microRNAs targeted at least 61 genes coding for extracellular matrix proteins and their regulators other than collagens and collagen regulators (Supplementary Dataset 2, see GO terms of target genes). It has been observed that abundance of collagen and connective tissue in extracellular matrix determines stiffness and gaping in fish muscle 37 . These DE microRNAs, which target collagens and other extracellular matrix connective tissues and their regulators, may be used to develop suitable genetic markers to improve muscle firmness in rainbow trout.
Muscle whiteness index. Seventy-seven microRNAs were DE, 58 upregulated and 19 downregulated, in fish families with high muscle whiteness index compared to low whiteness index families ( Fig. 2 and Supplementary Dataset 1C). From DE microRNA, correlation between microRNA expression level and variation in phenotype was performed for 11 microRNAs, of them, variation in expression level of 6 microRNAs significantly correlated with variation in muscle whiteness index (  this is among the first genome-wide studies aimed at exploring the genetic/molecular basis of muscle whiteness in salmonids.

MicroRNA-target gene co-expression and transcriptional regulation. MicroRNAs and their target
genes are usually co-expressed and co-regulated by common transcription factors 39 . To investigate co-expression of DE microRNAs and their target genes, first, we calculated Pearson correlation between each DE microRNA and their target genes based on their expression values in 30 RNA-Seq samples (see methods section). Next, we investigated whether this co-expression was regulated via common transcription factors (TF). For this purpose, we scanned promoter sequences of the strongly correlated (R 2 > 0.70) DE microRNAs and their target genes for TF binding cis regulatory motifs. Out of 90 DE microRNAs, 15 microRNAs had strong positive expression correlation (R 2 > 0.70) with 194 different target genes, and these correlated microRNA-target gene pairs shared common TF binding motifs in their promoters. Selected microRNA-target gene pairs and common TF binding motifs in their promoters are given in Table 4. Some of those TFs are known to be heavily involved in muscle development (e. g. myoD, c-Fos, c-Jun, MAZ, NF-AT1, Smad3, Elk1, PEA3, NFI/CTF and NFY), development (e.g. HOXD9, HOXD10 and COE1), metabolism (e. g. HNF-3) and adipogenesis (e. g. C/EBPbeta). These findings suggest that myogenic TFs may regulate expression of muscle related microRNAs and their target genes. Wang and coworkers suggested that both miRNAs and TFs must stay active to simultaneously regulate their target genes 39 .
Genetic polymorphism in microRNA target sites. In order to explore phenotype-associated genetic polymorphism on microRNA and microRNA binding sites in target genes, first we predicted target genes of all rainbow trout microRNAs that were reported previously 18 . We performed SNP prediction on the mature microRNA sequences. In addition, SNPs on microRNA target genes were predicted by mRNA sequencing of the same fish families used in this study (data under review for publication elsewhere). A total of 249 SNPs existing in microRNA recognition element seed site (MRESS) of target genes showed allelic imbalance (>2.0 as an amplification and <0.5 as loss of heterozygosity) between high and low ranked families of different muscle traits. Out of 249 SNPs, 240 SNPs either destroyed or created a novel illegitimate microRNA target, and only 9 SNPs had no potential effect on microRNA binding (Supplementary Dataset 5). This alteration in target recognition is due to presence of the SNPs in MRESS which plays crucial role in microRNA-mRNA binding 1 . Any SNP in the MRESS has an important impact on microRNA recognition 7 . Out of 240 SNPs capable of destroying or creating a novel illegitimate MRESS, 204 SNPs were found to be true polymorphic SNPs by genotyping fish, and showed allelic imbalance between high and low ranked families (>2.0 as an amplification and <0.5 as loss of heterozygosity). Interestingly, 78 unique SNPs showed significant association with growth and muscle quality phenotypes based on genotype-phenotype association analysis performed on a larger number of fish from 2 generations (n = 786) (FDR-p-value < 0.05). Table 5 lists selected SNPs with high correlation coefficients with growth traits while complete list of SNPs in microRNA binding sites and their association with phenotype is given in Supplementary Dataset 5. A total of 55, 28, 7 and 9 MRESS-destroying SNPs in 3′ UTR were significantly associated with WBW, muscle yield, crude fat content and muscle whiteness respectively (FDR-p-value < 0.05). Previous reports indicated that SNPs in MRESS of target genes have important impact in phenotype determination 7 . Phenotype associated MRESS-destroying SNPs were present in various classes of genes including metabolic enzymes, transporter proteins, transcription factors, signaling molecules and muscle related proteins. Among the WBW associated SNPs, those present in 3′ UTR of troponin, malate dehydrogenase, ranbp-type and c3hc4-type zinc finger-containing protein 1, phosphate carrier mitochondrial like protein, cytochrome b-c1 complex, novel protein vertebrate nebulin, ATP-dependent 6-phosphofructokinase, ankyrin repeat and socs box protein 5 and calsequestrin-1 each SNP explained more than 2% variation in WBW phenotype (Table 5   muscle yield. Presence of MRESS-destroying alleles may stabilize the target genes against microRNA-mediated downregulation, which may contribute to the difference in phenotype between high vs low ranked families. As an example, allele destroying MRESS in troponin fast skeletal muscle was 11 times more frequent in high WBW family compared to low WBW family, and explained 4.2% variation in WBW (Table 5). Previous study performed in this trout population has also identified 3 SNPs present in 3′ UTR of troponin C associated with growth traits in trout 40 . Similarly, allele destroying MRESS in ATP-dependent 6-phosphofructokinase, a glycolytic gene, was 5.5 times more frequent in high WBW family compared to low WBW family, and explained 2.6% variation in WBW. Consistent with our finding, SNPs associated with growth traits in trout from our previous study were mainly present in genes involved in energy metabolism and muscle structure 40 . Above findings suggest that 3′ UTR SNPs capable of destroying or creating an illegitimate MRESS may have important functional consequences. In contrast to target genes, no SNPs with allelic imbalances were detected in the microRNA mature sequences. This may be due to high degree of negative selective pressure as hundreds of genes are regulated by the same microRNA 3 . Previous studies from mammals suggest that SNPs in microRNA sequence, especially seed region, are rare due to high degree of selective pressure 41 . These findings suggest that due to high selective pressure to conserve microRNA sequence, cell may introduce genetic variation in microRNA binding site of appropriate target gene to regulate the phenotype.

Conclusion
Improvement of muscle growth and quality in salmonids has long been sought by aquaculture industries. So far, little progress has been made toward genetic improvement of muscle growth and quality in salmonids. Muscle quality traits are lethally-measured traits which are hard to include in a breeding objective. In addition, complex interrelated relationship among muscle quality traits 19,42 , and influence of several genetic and non-genetic factors complicates genomic selection. Previous studies have found genetic variation for muscle traits in trout and other salmonids [43][44][45] . Genomic technologies could be used to exploit within-family genetic variation for these lethally-measured muscle quality traits to identify suitable genetic markers. In this study, we investigated microRNA expression and genetic polymorphism in microRNA-binding sites in target genes associated with phenotypic variation in muscle growth and quality in a pedigreed rainbow trout population undergoing selection for improved growth performance.
Muscle specific myogenic microRNAs (e. g. mir-1, mir-206 and mir-133) as well as several non-muscle specific microRNAs showed regulated expression in response to variation in muscle yield and other muscle quality phenotypes. Most of the DE microRNAs between high vs low muscle yield families were also DE between high vs low families of muscle whiteness and muscle shear force phenotype. This observation may be, in part, due to interrelated nature of these muscle growth phenotypes 19,42 . Biological pathways such as development, cell cycle, muscle growth, transcription and muscle proteolysis were significantly overrepresented in the list of DE microRNA targets suggesting that DE microRNAs may regulate growth and muscle quality traits by post-transcriptionally regulating the genes involved in these pathways. Presence of common cis regulatory motifs for myogenic TFs in DE microRNAs and their respective target genes may suggest that myogenic TFs may regulate expression of both myogenic microRNAs and their target genes.
Due to crucial role of MRESS in microRNA-mediated gene regulation, SNPs creating or destroying MRESS in target genes have important functional consequences 7,46 . In this study, we identified 72 true SNPs in 3′ UTR capable of abrogating or creating MRESS on several metabolic and growth-related genes, which explained significant variation in growth, muscle yield and other muscle phenotypes. However, functional significance of these SNPs in microRNA target recognition needs to be experimentally validated. To make the present study more robust in identifying potential genetic markers for muscle and growth phenotypes, we applied dual approach by analyzing differential microRNA expression as well as genetic variation analysis in their target genes. We believe that this approach is more productive as high negative selective pressure imposes constraints on expression/genetic variation in microRNA 41 . In support of this argument, we found several SNPs on microRNA binding site of target genes that showed allelic-imbalance between high and low WBW groups though microRNA itself did not show variation in sequence or expression between the two groups. In this study, we performed a genome wide approach to investigate variation in expression and variation in target recognition of growth and muscle-important microR-NAs, and the study may help identify suitable genetic marker for genetic selection of these traits.
In a separate study performed in the same set of fish, an unexpectedly low number of protein coding genes were observed as DE between high and low ranked families of these phenotypes (data will be published elsewhere). But, in this study, we observed significant number of DE microRNAs and their target protein coding genes that are associated with these phenotypes. These findings suggest that variation in these muscle quality phenotypes may be explained largely by variation in microRNA expression and genetic variation affecting recognition of microRNA targets rather than by direct regulation of mRNA expression. Consequently, the current study pinpoints to a greater role of microRNA-mediated gene regulation in determination of muscle traits in rainbow trout.
Using genomics approach, we have identified significant correlation of microRNA gene expression and SNPs in microRNA binding sites with growth/muscle quality traits of economic importance. Some of the genotype x phenotype associations reported in this study although significant, present low-to-medium values, not allowing to make definitive inferences about the behavior of some characteristics of the muscle growth and quality traits. The highly polygenic nature of traits demands identification of more genetic markers. Incorporation of genomics selection, based on genetic markers, can greatly help traditional breeding programs, based on quantitative evaluation techniques, in improving accuracy of breeding value prediction.

Materials and Methods
Ethics statement. Fish were maintained at the NCCCWA and all experimental protocols and animal procedures were approved and carried out in accordance with the guidelines of NCCCWA Institutional Animal Care and Use Committee Protocols #053 and #076.
Fish population and muscle sampling. Phenotypic data and muscle samples were collected from 471 fish representing 98 families (~5 fish/family) from USDA/NCCCWA growth selected trout line from each harvest year 2010 as previously described 40 . Briefly, fish families were produced and reared till ~13 month post-hatch as previously described 17 . Single-sire x single-dam mating was used to produce full sib families and eggs were reared in spring water. Water temperature was manipulated from 7 °C to 13 °C to synchronize the hatching time. Each family was reared in a separate 200-L tank at ~600 alevins/tank density. Random culling of fish was performed every month to maintain stock density of <50 kg/m 3 . At the age of 5 month, each fish was given unique identification PIT (passive integrated transponder) tag, and tagged fish were combined and reared in a big 1,000-L commercial tanks. Commercial fishmeal-based diet (16% fat, 42% protein; Ziegler Bros Inc., Gardners, PA) was fed using automatic feeder. Feeding rate was gradually reduced from 2.5% of body weight to 0.5% of body weight as fish grew older. WBW of all fish belonging to 98 families was measured, and families were ranked based on their WBW measurements. Second or third ranked fish from each family was selected for muscle sampling so that WBW of sampled fish is adjusted around median of each family. Selected fish were randomly assigned to one of the 5 harvest group (~100 fish/harvest group) and each harvest groups were sampled in 5 consecutive weeks. Fish were anesthetized in 100 mg/L tricaine methanesulfonate and weighted, slaughtered and eviscerated. Muscle samples were separated from dorsal musculature and were stored in liquid nitrogen until processing and phenotype measurement. Muscle quality phenotypes were measured at West Virginia University (WVU), meat processing laboratory. Muscle yield was calculated as percentage of WBW, and proximate analysis of muscle fillet was performed according to previously described protocol 47 . Briefly, fresh fillet surface color was measured with a chromameter (Minolta, Model CR-300; Minolta Camera Co., Osaka, Japan) calibrated using a standard white plate No. 21333180 (CIE Y 93.1; ×0.3161; y 0.3326). L* (lightness), a* (redness), and b* (yellowness) values were recorded at three locations above the lateral line along the long axis of the right fillet, and these values were used to calculate a fillet whiteness index according to the following equation; whiteness index = 100 − [(100 − L) 2 + a 2 + b 2 ] 1/248 . Muscle crude fat content was measured using Soxhlet solvent extractor with petroleum ether. For muscle peak shear force measurement, texture was analyzed using a five blade Allo-Kramer shear attached to the texture analyzer; and Texture Expert Exceed software (version 2.60; Stable Micro Systems Ltd., Surrey, U.K.) was used to record the peak shear force.
Library construction and sequencing. For differential gene expression analysis, 4 highest ranked and four lowest ranked families for each trait were sequenced. White muscle sample was isolated from five individuals belonging to each family and total RNA was isolated using Trizol protocol (Invitrogen, Carlsbad, CA) as described previously 49,50 . To allow enough transcriptome sequence total RNA from 5 individual fish from each family was pooled and sequenced on Illumina's HiSeq platform (Illumina Inc, CA, USA).

D a t a p r o c e s s i n g a n d p r e d i c t i o n o f t r o u t m i c r o R N A .
S e q u e n c i n g a d ap t e r 5′GCCTTGGCACCCGAGAATTCCA-3′ was trimmed and reads were annotated using miRBase microRNA reference (release 21) in CLC Bio small RNA analysis tool. Read alignment was run at default settings (i.e. mismatch ≤ 2, additional/missing upstream/downstream bases ≤ 2). MicroRNAs with mismatches and/or additional/missing upstream/downstream nucleotides were considered as variants of the same microRNA. Read count from all variants of the same microRNAs were summed and were used as expression value for that particular microRNA (default method of merging expression values in CLC Bio's small RNA analysis toolkit).
Identification of DE microRNAs. DE microRNAs were identified using EDGE test in CLC genomics workbench using expression values from the above step. The fold change in gene expression between two groups was considered significant if FDR-p < 0.05 and fold change <−2 or >2 fold.

Real time PCR validation of DE microRNAs and 'microRNA expression-phenotype variation'
correlation. Same RNA samples from high and low ranked families used for sequencing were used to validate RNA-Seq differential expression. cDNA was synthesized using miScript II RT kit (Qiagen, Valencia, CA, USA) and microRNA was quantified in Bio-Rad CFX96 ™ Real Time System (Bio-Rad, Hercules, CA) using miScript R SYBR R green (Qiagen, Valencia, CA, USA). A non-coding RNA U6 was used as a endogenous control for normalization, and fold change was calculated by ΔΔCt method 51 , as described previously 52,53 .
Correlation between microRNA expression and phenotype was studied in 90 random individual fish from 3 rd -generation families (2010 hatch year) from the growth-selected line. ΔCt value of each microRNA in all 90 fish was calculated. Pearson correlation and simple linear regression was run between ΔCt and quantitative value of muscle phenotype.
Bioinformatics prediction of DE microRNA target. For prediction of microRNA targets, 3′-UTR of trout mRNA were retrieved from the genome reference 54 . Due to difference in sensitivity and specificity of different target prediction algorithms, targets were predicted using 3 tools: miRanda, PITA and TargetSpy in small RNA analysis server sRNAtoolbox 55 . If the same target site is predicted by all 3 tools it was considered a potential microRNA target. For PITA, prediction parameters chosen were: seed length 6-8 nts, no G:U wobble allowed in seed of size 6 nts, one G:U wobble allowed in seed of size 7-8 nts, no mismatches allowed in seed of size 6 and 7 nts, one mismatch allowed in seed of size 8 nts, and no loop is allowed in microRNA or target for any seed size. Miranda parameters chosen for target prediction were: score threshold 150, gap-open penalty −4.0 and gap-extend penalty 9.0. For TargetSpy score a minimum threshold of 0.99 was used. For all tools, minimum energy threshold was chosen as −15Kcal/mole. Gene enrichment analysis of microRNA targets. Gene enrichment analysis (GEA) was performed by using DAVID 20, 21 (FDR-p < 0.05) and overrepresented pathways were visualized by EnrichmentMap app 56 in cytoscape 57 (p value < 0.005 and FDR-q < 0.05). Overlap between gene-sets was computed according to overlap coefficient, which was set to the default recommended value of 0.5.

MicroRNA-target gene co-expression and identification of cis regulatory promoter motifs.
Small RNA and mRNA sequencing reads from 22 families were used to calculate microRNA-target gene correlation. For target genes, mRNA sequencing reads were mapped to trout mRNA reference and transcript per million (TPM) was calculated for each mRNA. For microRNA, small RNA sequencing reads were mapped to mature microRNA sequence from miRBase release 21 (June 2014) and total count was calculated for each microRNA. TPM (for mRNA) and total count (for microRNA) were normalized and used to calculate expression correlation.
Transcription factor (TF) binding motifs were searched in the 500 upstream promoter sequences of DE microRNAs and correlated target genes using.Alggen Promo TF motif search tool 58,59 . Search parameters used were 'only teleost transcription factors' and 'only teleost motif sites' . Maximum dissimilarity rate between putative and consensus TF binding site was set to 5%, and RE equality/query (expectation of finding motif in random sequence) was set to <0.05.

Prediction of single nucleotide polymorphism (SNP). The same RNA samples used to sequence small
RNAs were used to sequence protein coding genes. Details of SNP prediction are described in a publication under review. Briefly, sequencing libraries were prepared using Illumina's Truseq RNA library preparation kit and sequenced on Illumina's HiSeq platform (Illumina Inc, CA, USA). GATK 60 and SAMTool 61 were both used to call SNPs. In the SAMtool approach, STAR alignment tool 62 was used to align sequencing reads from each family to the trout genome. SAMtools view/sort and mpileup functions were used to determine variant genotype. Popoolation2 package (version 1.201) was used to calculate allele frequencies 61,63 . Initial SNPs were considered at minimum reads >10 and minor allele count >4 and minimum allelic frequency (MAF) >0.05. SNPs were considered as putative trait-associated SNPs if allele frequency (allele A/allele B) ratio between high and low ranked families was ≥2.0 (as an amplification) or ≤0.5 (as loss of heterozygosity). In GATK approach, reads were aligned to the genome using STAR alignment tool and then Picard tool was used to sort the SAM files and to mark duplicates. Split and trim function was performed to reassign mapping quality, Indel realignment and local realignment around the Indel to clean up any mapping artifacts. Finally, base quality score recalibration was performed on the data resulting from Picard tools. Haplotype Caller was used to determine variants, followed by filtration of SNPs using strict thresholds: Qual By Depth (QD) 2. SNP genotyping and SNP-phenotype association analysis. SNP genotyping was performed as a part of development of a 50 K SNP chip for rainbow trout (full SNP chip study will be published elsewhere). Putative SNPs were genotyped in ~1,900 individuals from USDA, NCCCWA growth selected lines harvest year 2010 and 2012 using Affymetrix SNP array protocol (Geneseek Inc., Lincoln, NE, USA). Briefly, genomic DNA from each fish was extracted from a fin clip and was amplified by PCR. DNA samples were chemically fragmented and then biotinylated. Biotinylated DNA samples were hybridized to DNA probe in the SNP array and genotype call was determined based on 'probe DNA: sample DNA' hybridization.
For SSNP-phenotype association study, only 249 SNPs present in 3′ UTR microRNA binding sites and only 786 individual fish from 3 rd -and 4 th -generation families (2010 and 2012 hatch years) 17 with available phenotype measurements of interest were considered. Phenotypic data were rank transformed 64 , then data normality was checked by Kolmogorov-Smirnov and Shapiro-Wilk test to make sure that the traits are normally distributed to meet the assumption of quantitative trait association analyses. SNP-phenotype association was performed by quantitative trait analysis, linear and logistic methods in PLINK tool 65 and all three methods gave consistent results, and results from quantitative trait analysis were reported. Quantitative association analysis was performed between SNP genotype and quantitative value of phenotype measurements of 786 individuals. Linear association between SNP genotypes and phenotype measurement was performed to get correlation coefficient (R 2 ) between the SNP allele and the phenotype. Logistic regression was performed to check if the effect of one SNP on the phenotype was independent of other SNPs. Associations of all reported SNPs were significant in linear and quantitative analysis after entering other SNPs as covariates in logistic method, suggesting that the effect of each SNP for association to the phenotype is independent from other SNPs. Population stratification analysis was performed using genome-wide IBS that is integrated in PLINK. Multidimensional scaling plot on N x N matrix of genome-wide IBS pairwise distances showed fairly homogeneous sample with no obvious population stratification or significan clustering (Supplementary Dataset 1E).