Comparisons among rainbow trout, Oncorhynchus mykiss, populations of maternal transcript profile associated with egg viability

Abstract

Background

Transcription is arrested in the late stage oocyte and therefore the maternal transcriptome stored in the oocyte provides nearly all the mRNA required for oocyte maturation, fertilization, and early cleavage of the embryo. The transcriptome of the unfertilized egg, therefore, has potential to provide markers for predictors of egg quality and diagnosing problems with embryo production encountered by fish hatcheries. Although levels of specific transcripts have been shown to associate with measures of egg quality, these differentially expressed genes (DEGs) have not been consistent among studies. The present study compares differences in select transcripts among unfertilized rainbow trout eggs of different quality based on eyeing rate, among 2 year classes of the same line (A1, A2) and a population from a different hatchery (B). The study compared 65 transcripts previously reported to be differentially expressed with egg quality in rainbow trout.

Results

There were 32 transcripts identified as DEGs among the three groups by regression analysis. Group A1 had the most DEGs, 26; A2 had 15, 14 of which were shared with A1; and B had 12, 7 of which overlapped with A1 or A2. Six transcripts were found in all three groups, dcaf11, impa2, mrpl39_like, senp7, tfip11 and uchl1.

Conclusions

Our results confirmed maternal transcripts found to be differentially expressed between low- and high-quality eggs in one population of rainbow trout can often be found to overlap with DEGs in other populations. The transcripts differentially expressed with egg quality remain consistent among year classes of the same line. Greater similarity in dysregulated transcripts within year classes of the same line than among lines suggests patterns of transcriptome dysregulation may provide insight into causes of decreased viability within a hatchery population. Although many DEGs were identified, for each of the genes there is considerable variability in transcript abundance among eggs of similar quality and low correlations between transcript abundance and eyeing rate, making it highly improbable to predict the quality of a single batch of eggs based on transcript abundance of just a few genes.

Background

Egg quality is fundamental to reliable seed stock production in aquaculture and yet what makes an egg developmentally competent to be fertilized and subsequently develop into a normal embryo is poorly understood [1,2,3]. Fertilization rates are often high in the rainbow trout industry but the quality of eggs in fishes can be affected by intrinsic factors such as the genetics and age of the brood fish [1, 4,5,6,7,8,9,10] and extrinsic factors that can vary with hatchery environments and practices [11,12,13,14,15]. Female rainbow trout broodstock do not volitionally oviposit in captivity and therefore must be stripped of their eggs following ovulation. The female gamete obtained by this stripping process or when spawning naturally is an oocyte arrested in metaphase of the second meiotic division that should be competent for fertilization. The oocyte is largely transcriptionally silent from the end of oocyte growth until the zygote genome is activated, referred to as zygotic genome activation (ZGA), which begins at about the mid-blastula transition (MBT) in most vertebrates. The oocyte therefore serves as a reservoir for RNAs as well as other biomolecules including proteins and lipids accumulated during oogenesis, for utilization from oocyte maturation through early embryonic development [16, 17]. Levels of biomolecules in the egg including proteins, lipids, and RNAs have been linked to egg viability in many fishes including rainbow trout [1,2,3, 18].

The almost total reliance of the late stage oocyte and early embryo on maternally derived RNAs has led to investigations of associations between the maternal transcriptome and measures of developmental competence in several species of fish and has been reviewed [3, 19, 20]. Most investigations identified mRNAs that reflect differences in egg quality by simply comparing transcript expression profiles among eggs or early embryos exhibiting variation in measures of developmental competence, usually including progression to a specific developmental stage or a developmental abnormality [21,22,23,24,25,26,27]. A number of studies, primarily in rainbow trout, have identified mRNAs differentially expressed among eggs of different quality in response to treatments used to alter time of spawning through photoperiod manipulation or hormone treatment [15, 28] and in response to being overripe due to post-ovulatory aging [29, 30]. In addition to mRNAs, profiles of microRNAs and mitochondrial genome-encoded small RNAs were related to egg deterioration caused by post-ovulatory aging in rainbow trout [31, 32]. Recently, we identified over 1000 differentially expressed transcripts or genes (DEGs) in unfertilized rainbow trout eggs that are associated with eyeing rate [27]. However, these differences were only found when the libraries used for sequencing were prepared following polyadenylation capture and not rRNA-removal, suggesting differences in egg quality may derive in part from differences in maternal transcript activation and cytoplasmic polyadenylation before ovulation.

Much has been learned about the contribution of maternal mRNAs to egg quality in fish. As expected, many of the transcripts that appear dysregulated in poor quality eggs are in pathways known to be involved in critical processes taking place at the developmental stages investigated [3, 19, 20]. Nevertheless, there is considerable disparity in DEGs identified among the studies. This may be due to differences in species, stages investigated, measures of egg quality, intrinsic and extrinsic causes of the differences in quality, and molecular and statistical approaches employed. Furthermore, studies thus far have focused on identifying possible DEGs for dysregulation but compared transcriptomes of few individuals. The aim of the present study is to further evaluate the robustness of genes identified as possible markers of egg quality using a commercially important species, rainbow trout. To meet this aim we designed an assay based on the nCounter analysis data system (Nanostrings Technologies; Seattle, WA) to compare expression of 65 mRNAs previously identified as being differentially expressed with egg quality (Additional file 1: Table S1). The nCounter analysis data system was chosen in part because it is relatively easy to customize the multiplex CodeSet to update or meet specific needs of the user. Most of the genes incorporated in the assay are DEGs from our previous transcriptome analysis of egg viability in rainbow trout using RNA-Seq, [27], but also includes 10 additional transcripts reported as dysregulated in poor quality eggs in rainbow trout [28,29,30], and also igf-3 since many IGF-system genes were already in the assay. The genes from [27] were selected for the assay primarily based on magnitude of statistical differences and fold-change. Three populations of broodstock were compared including two different year classes from a commercial line, referred to as Group A1 and A2 respectively, and females from the 2015 year-class at the National Center for Cool and Cold Water Aquaculture (NCCCWA) referred to as Group B. One of the year classes from the commercial line, Group A1, included eggs from the same females used in our RNA-Seq study [27]. In all 152 families were included in the study. The present study had four aims. The first aim (i) was to determine if DEGs identified in a limited number of fish were DEGs in a broader sample; the second aim (ii) was to determine if the identified DEGs were consistent year to year within the same line; the third aim (iii) was to determine if they varied from line to line, and the fourth aim (iv) was to determine if a small set of genes can act as a reliable universal marker for egg quality in rainbow trout.

Results

Eyeing rate and early embryo viability

Eyeing rate was assessed at ~ 250 accumulated temperature units (ATUs) post fertilization. This timepoint is slightly after retinal pigmentation but often used by hatcheries because embryos are resistant to handling or mechanical shock [33, 34], most of embryonic mortality has already occurred [35], and it still allows time for dead and subviable egg removal and shipment to hatching facilities. Eyeing rates were collected for all families in each of the broodstocks that made up that year’s cohort for genetic selection for that line. A total of 192, 143, and 325 families were evaluated for Groups A1, A2, and B respectively, with mean eyeing rates of 78.3% + 0.015, 79.1% + 0.015, and 49.7 + 0.017 (Fig. 1abc). The data on eyeing rate for Group A1 were previously presented in Ma et al. [27]. Historical eyeing rates are higher for the commercial hatchery lines from which groups A1 and A2 were collected, than for the NCCCWA line from which group B was collected. Nevertheless, there were fewer egg lots with survival less than 30% than has usually been observed (Kyle Martin, personal communication) with only 6 and 2 families yielding eyeing rates below 30% in Groups A1 and A2 respectively, and all these were below 10%. Transcript abundance analysis was determined for 48, 44, and 60 families for Groups A1, A2 and B respectively including all families with less than 30% eyeing in Groups A1 and A2 (Fig. 1def).

Fig. 1

Eyeing rates of all the surveyed rainbow trout families in the breeding groups (a-c) and those selected for mRNA analysis (d-f)

Families from Group A1 with eyeing rates under 80% were generated with sperm that also generated families with eyeing rates over 78%, supporting subfertility was due to the eggs and not the sperm. Sperm used in Group A2 to fertilize each of the 27 families with eyeing rates between 20 and 80% also produced families with eggs from a different female that yielded 22 families with eyeing rates over 70% and 18 over 80% support eyeing rates were mainly due to egg quality. Although the sperm lot used to produce the family with an eyeing rate of 0% used in the present study also yielded a family with an eyeing rate of 83.1%, sperm from the family with 1.4% eyeing in the present study was not used to produce a second family making it unclear if the low eyeing rate was due to the quality of the egg or sperm, although normalized read values are consistent with reduced egg quality (Additional file 1: Table S5B). Sperm quality could not be ruled out as contributing to eyeing rates in Group B since sires were only used once. Egg lots used in the study showed no obvious visible signs of poor quality including overripening when examined before fertilization.

Mortality before eyeing has been previously investigated in line A including for 20 of the families used in Group A1, and found to predominantly take place before the 32-cell stage [27, 36]. In the present study embryo cleavage was assessed at about 19–20 h post fertilization at ~ 10 °C and early embryo development or streak rate was estimated at about 10 days post fertilization for the 60 families in Group B (Table 1; Additional file 1: Table S2). Fertilization rate was high with families averaging 89.6% of zygotes completing first cleavage. The majority of the zygotes of families with eyeing rates greater than 80%, which we consider families with high quality eggs, reached at least the 16-cell stage, 91.6%, with some reaching the 32-cell stage, 43.2%. Those zygotes not reaching the 8-cell stage were therefore considered subviable and on average 76.7% of zygotes reached this stage. This is well above the mean eyeing rate of 35.9%. We prefer assessing early stage mortality after most of the embryos in the families with greater than 80% eyeing rates reach the 32-cell stage, which we failed to meet in Group B samples. We therefore included a measure of streak rate which as evaluated is only a rough estimate of development to an elongating embryo. The average streak rate among the families was 63.7% which is still well above the eyeing rate supporting mortality was taking place throughout development to eyeing in Group B.

Table 1 Assessment of early embryo development in Group B selected families

Transcriptome abundance analysis

Overall, there were 32 transcripts identified as DEGs among the three populations or groups by regression analysis (Tables 2, 3, 4; Fig. 2a-c; Additional file 1: Tables S3AB). More DEGs were shared between Groups A1 and A2 which were within the same line, than between these groups and Group B which is from a different line. Group A1 (Table 2) had the most, 26; A2 (Table 3) had 15, 14 of which were shared with A1; and B (Table 4) had 12, 7 of which overlapped with A1 or A2. Six transcripts, all from nuclear genes, were found to be differentially expressed in all groups (Fig. 2a-c; Table 5). Low raw read counts limited the detection of differences in the same 10 genes in each of the three groups and two additional genes among the groups (Additional file 1: Tables S4A-D).

Table 2 Group A1 Normalized reads

Table 3 Group A2 Normalized reads

Table 4 Group B Normalized reads

Fig. 2

Normalized reads for nuclear genes that were differentially expressed with egg quality. a) Group A1. b) Group A2. c) Group B. Values are mean + SEM for families designated as having low-quality, medium-quality, or high-quality eggs based on eyeing rates

Table 5 Summary of gene transcripts with significant correlations between transcript abundance and eyeing rate

In Group A1 regression analysis identified 25 nuclear and one mitochondrial gene with transcript levels correlated with eyeing rate (Table 2; Fig. 2a). Twenty-two of the nuclear genes and the mitochondrial gene, mt-cyb, had increased transcript abundance with increased survival and three decreased. The coefficient of determination ( R[2]) values for normalized untransformed data were at or below 0.2269 for all genes. Three genes, impa2, linb7, and mrpl39-like had over three times more transcripts in the high-quality eggs (80–100% eyeing) than in the low-quality eggs (0–20% eyeing). There were five genes in which the medium-quality eggs (20–80% eyeing) had the highest and one the lowest number of reads, and apoc1 had about three times more abundant reads than either the low- and high-quality eggs which had read amounts similar to each other. The numerical means for all the mitochondrial genes in the high-quality eggs were 46–105% above the low-quality eggs.

In Group A2 there were 14 nuclear genes and one mitochondrial gene with correlated transcript abundance and eyeing rates (Table 3; Fig. 2b). All but fbxo5 were also significant for A1 (Table 5). Transcript abundance and eyeing rates were positively correlated for all DEGs and R² values were at or below 0.1878. As in A1, impa2, linb7, and mrpl39-like had over three times more transcripts in the high-quality eggs than in the low-quality eggs, as did samm50 in A2. There were no DEGs in which the medium-quality eggs had the highest or lowest number of reads. The numerical means for all the mitochondrial genes in the high-quality eggs were 58–143% above the low-quality eggs.

In Group B Regression analysis identified 11 nuclear and one mitochondrial gene with transcript levels correlated with eyeing rate (Table 4; Fig. 2c). Transcript abundance of seven of the nuclear genes increased with eyeing rate whereas the remaining four along with the mitochondrial gene mt-dlp, decreased. Six of the nuclear genes were also significant for both A1 and A2, nasp was also significant for A2, whereas the remaining 4 and the mitochondrial gene mt-dlp were only significant for B (Table 5). Furthermore, transcript abundance of all the DEGs shared with A1 or A2 were positively correlated with eyeing rate whereas the remaining transcripts including mt-dlp, were all negatively correlated. The R² values were at or below 0.2075 for the DEGs and differences among low-, medium- and high-egg quality family means were less than 2-fold for all DEGs. The numerical means for all the mitochondrial genes other than mt-dlp were 25–65% greater in the high-quality eggs than the low-quality eggs.

Discussion

Assay performance

The nCounter analysis data system was selected with possible use by hatchery managers in mind. Nanostring Technologies designs the custom CodeSets for nCounter® analysis based on submitted target sequences and conducts the genomic analyses required to avoid non-specific hybridization, and provides free software, nSolver, to quality check, normalize, and analyze the data. In addition, the system can include over 800 genes and new probes can easily be exchanged in the CodeSet. In our assay the mean raw reads for a given gene were consistent across the three groups or populations with an average CV of 18.6% (Additional file 1: Table S4D) comparing among the three groups. The high abundance of reads for the five mitochondrial genes limited read values for many nuclear genes with 14 below acceptable limits including two of the reference genes ef1a and ppia (Additional file 1: Table S4D). The CodeSet can be designed to attenuate high abundance genes to alleviate this problem. The cause for apparent instability of b-actin in A1 is not known but it was also found to be unstable among the 20 samples from A1 previously analyzed by RNA-Seq [27]. Nevertheless, the reliability of reference genes can be questionable with egg quality in which the proportion of mitochondrial and nuclear transcripts can vary with the quality of the eggs [27] and the efficiency of methods to capture polyadenylated transcripts can vary for shorter poly(A) tail lengths as with stored maternal mRNAs. Whereas the majority of cytosolic nuclear transcripts in most cells are polyadenylated with a poly(A) tail greater than 80 nucleotides [37, 38], stored maternal nuclear transcripts possess a short poly(A) tail of around 15–40 nucleotides that are elongated to over 80 nucleotides through cytoplasmic polyadenylation during activation [37, 39,40,41,42]. Oligo(dT) capture approaches in general are not very efficient at capturing mRNAs with shorter poly(A) tails [39, 43, 44]. Assay performance and normalization may be improved with the use of an RNA spike-in as described by Bettegowda et al. [45] for use in bovine oocytes. An RNA with a poly(A) tail with a minimum of 25 nucleotides that would likely be efficiently capture by most oligo(dt) methods, and can be used starting at homogenization, should be considered.

Eyeing rate and early embryonic viability

Previous studies have characterized early mortality in line A. We characterized mortality in the 20 families from Group A1 when used in our previous analysis of transcript abundance [27]. In these 20 families almost all the mortality observed by eyeing, took place between the 8- and 32-cell stages. Stoddard et al. [36] studying the same line determined most of the mortality in subfertile embryos took place by the second cleavage interval. Although the timing was slightly different, in both studies most of the mortality took place before the MBT and therefore before ZGA. Although the timing of the ZGA has not been determined for rainbow trout, the major wave of ZGA in most fish investigated takes place after the 64-cell stage [46,47,48,49,50]. Furthermore, the fish species investigated develop more rapidly than rainbow trout and in general the number of cleavage divisions completed before ZGA is greater in animals that develop more slowly [51]. We normally see a more prolonged period of mortality in embryos from our NCCCWA broodstock as observed with the families selected for the present study. Unfortunately, we sampled the embryos in Group B earlier than when we sampled A1, and therefore were not able to determine what percentage of mortality took place before the 32-cell stage which was the end of most pre-eyeing mortality in A1. Nevertheless, there was considerable mortality between fertilization and the 8-cell stage; the 8-cell stage and the 16–32-cell stages at about 20 h; ~ 20 h and streak; and between streak and eyeing (Table 1). As mentioned, the streak rate is a very rough estimate, but the values support at least about 28% of mortality by eyeing took place after the ZGA, suggesting the causes for mortality are likely, at least in part, different for the A and B populations.

Gene expression

Mitochondrial genes

The rainbow trout mitochondrial genome encodes 13 polypeptides, two rRNAs, 22 tRNAs and a non-coding region [52]. All 13 polypeptide genes and the mt-dlp region transcripts were found to be reduced to a similar degree in eggs with low eyeing rates by Ma et al. [27]. The present assay included five mitochondrial genes, one for a polypeptide from each of the four complexes of the electron transport chain, and the non-coding mt-dlp region. In the zebrafish, mitochondria are required for oxidative phosphorylation to generate ATP for early cleavage events as the maternal ATP pool is insufficient to sustain the proteasomal pathway required for protein degradation needed to advance beyond the 32-cell stage [53]. Only two of the mitochondrial genes were found to be differentially expressed among the three studies, with mt-cyb being increased with eyeing rate in A1 and A2 and mt-dlp transcripts being negatively correlated with eyeing rate in B (Tablse 2, 3, 4, 5). The P-values reported in Ma et al. [27] for the differences in expression for the mitochondrial genes were generally much higher than for the nuclear genes selected to be in the present assay, which may also contribute to why a greater proportion of the selected nuclear genes were identified as DEGs in the present groups. Nevertheless, all the mitochondrial genes in all groups trended towards increasing with eyeing rate except for mt-dlp transcripts in group B, consistent with decreased mitochondrial gene expression being a common feature of eggs with reduced developmental competence. Monitoring a suite of mitochondrial genes for small but similar differences in expression may be a reliable approach to identifying eggs with compromised mitochondrial function.

Nuclear genes

A total of 48 nuclear genes were at detectable levels for all three groups. Transcript abundance was correlated with eyeing rate in 30 of these nuclear genes in at least one of the three groups (Table 5; Fig. 2a-c). All the DEGs except ctsz, cycB and mr-1 were among the 49 nuclear genes included in the assay because they were identified as DEGs in our previous RNA-seq study on Group A1. Of the 17 genes included in the assay that had previously been reported as DEGs in studies of rainbow trout other than just Ma et al. [27] [28,29,30], eight were identified as a DEG in at least one of the present groups, seven including igf-2 were not significantly altered, and two (igf-1 and npm2) were below the detection limit in all groups. Therefore, among the genes that were represented at detectable levels, about half of the transcripts selected from our previous study on Group A1, and half of the transcripts identified from studies by different investigators, served as markers for egg quality in at least one of our study groups supporting a fair degree of overlap in genes altered with egg quality among rainbow trout populations.

The highest number of nuclear DEGs, 25, was in Group A1 (Table 2; Fig. 2a) as would be expected since most of the genes in the assay were identified as DEGs using samples from this group [27]. About half as many genes were DEGs in A2 (Table 3; Fig. 2b), with 13 of the 14 being DEGs in both A1 and A2 and only fbxo5 significant in A2 and not A1 (Table 5; Fig. 2ab). All 26 genes identified as DEGs for either A1 or A2 in the present study, including fbxo5, were DEGs in Ma et al. [27]. Transcript levels of all 14 genes that were DEGs in both A1 and A2 were positively correlated with eyeing rate and only myo1b, galnt3 and tob1 were negatively correlated in A1. The mitochondrial gene mt-cyb was also increased with eyeing rate in both A1 and A2. This uniformity in the profiles of maternal transcript dysregulation with poor egg quality for both year classes for population A supports a consistent cause for reduced egg quality in this line. A consistent pattern of dysregulation increases the chances that further investigation can elucidate physiological reasons and an underlying condition for variation in egg quality among females of a broodstock. Again, these transcripts derive from a global transcriptome analysis of egg quality in this line [27].

Despite being from a different line than the one from which most of the genes in the assay were selected, Group B (Table 4; Fig. 2c) had a similar number of nuclear genes with expression levels that correlated with eyeing rate as A2, 11; of which seven were shared with A1 and A2 (Table 5; Fig. 2a-c). The four genes that were unique DEGs to Group B were also the only four that decreased in abundance with increased eyeing rate. One of the four genes that were unique DEGs to Group B, rpl30, was identified as a DEG in A1 following RNA-Seq analysis [27]. Therefore, whereas all the genes identified as DEGs in A1 and A2 were previously reported in [27], Group B had three genes that were not. These genes included the two highest expressed nuclear genes in the assay, cycB and ctsz. Transcript abundance of cycB increased with post-ovulatory ageing and malformation rate at yolk sac resorption, but not survival at eyeing [29]; and ctsz increased with post-ovulatory ageing and decreased with eyeing rate [30]. The abundance of the remaining unique DEG, mr-1, decreased in eggs from fish induced to ovulate with hormone injections; a treatment that also results in more deformed embryos at yolk sac resorption [28]. Whereas ctsz and cycB generally increased with measures or treatments associated with reduced egg quality among the studies, mr-1 which increased with eyeing rate in our studies, decreased in response to hormone induced spawning.

Two genes, tubb which increased with eyeing rate in A1 and A2, and krt18, which increased with eyeing rate in A1, were also identified as DEGs with eyeing rate by Aegerter et al. [30]. Bonnet et al. [28] found rpl24, myo1b, and apoc1, which were DEGs in A1, to have altered expression in response to treatments that changed spawning time, which were in turn shown to increase malformation rates at yolk sac resorption. However, the direction of the effects did not agree among the studies. Transcripts for tubb were increased with eyeing rate in both studies. On the other hand, krt18 decreased with increasing eyeing rate in Aegerter et al. [30] but was increased with increased egg eyeing rate in Group A1. Nevertheless, the R² was very low for A1, 0.0453 (Table 2), and the families with medium fertility had the highest mean transcript levels in Groups A1, A2, and B (Tables 2, 3, 4; Fig. 2a-c). Transcript abundance for myo1b was negatively correlated with eyeing rates in the present study and increased in eggs from fish induced to spawn with photoperiod shifting [28], whereas apoc1 and rpl24 were positively correlated with eyeing rates in the present study but were increased in eggs of fish injected to spawn with a gonadotropin-releasing hormone analog [28].

Many genes that were not significantly correlated with eyeing rate or were below detection in the present study were identified as DEGs in the previous studies of egg quality in rainbow trout. Genes for which no correlation between transcript abundance and eyeing rates were found in the present study included genes of the IGF system, igf-2 and igfr1b. Aegerter et al. [29] found igf-1, which was below detection in the present assay, and igf-2, to be associated with eyeing rate, and although igfr1b was not associated with eyeing rate, it was decreased with increased malformation rates at yolk sac resorption. Transcript abundance of krt8 and ptgs2 were not altered in our groups although they were negatively correlated with eyeing rate by Aegerter et al. [30]. Transcript abundance of npm2 and igf-1 were increased with eyeing rate in the same study but were below detection in our assay. Detection limits for several of these genes in the present study impaired a comparison among studies. Transcript abundance of pyc and ntan1 were not altered in our groups but were increased in fish exposed to a shifted photoperiod and hormone injection to affect time of spawning, respectively [28].

In all, differences in transcript profiles can be found among all studies and populations in which a relation between egg quality and maternal transcriptome have been investigated in rainbow trout. The bases for these differences among studies are not known but there were differences in when mortality took place and treatments associated with decreased developmental competence among the studies that likely influenced or were influenced by the transcript profiles. In the present study almost all the mortality that would take place by eyeing likely took place by 24 h in line A based on studies by Ma et al. [27] on Group A1 and Stoddard et al. (2015), looking at an earlier year class of the same line. Later mortalities in Group B suggests at least some different factors associated with mortality before eyeing. Further analyses including time of mortality within population B families may suggest if there are transcript profiles associated with early and later embryonic mortality. Aegerter et al. [29, 30] used post ovulatory ageing to decrease egg quality whereas it was purposely avoided as in the present study. Genes identified as DEGs in the studies on post ovulatory ageing that were not DEGs in the present study or had opposite correlations might suggest those genes that are more specific or responsive to post ovulatory ageing. This includes krt8, krt18, ptgs2, and possibly igf-2 or the IGF system in general. Different expression profiles between the present study and the study by [28] may not be surprising considering the transcript levels in [28] were in response to treatments that were in turn associated with increased malformations at yolk sac resorption as a measure of egg quality. Thus, the transcript profiles in this previous study were not only not in association with the same egg quality metric as in the present study, they were not even in direct association to that egg quality metric, malformations at yolk sac resorption. Nevertheless, the differences in responses of many genes among the studies shows these transcripts cannot be used as general markers of egg quality but does not negate important roles for these transcripts in embryo development.

The identification of six DEGs shared among our three groups, dcaf11, impa2, mrpl39_like, senp7, tfip11 and uchl1, supports some consistency among populations in transcripts that are altered with the egg quality metric eyeing rate. Considering the six genes were identified as DEGs in our previous study that included in-depth genomic analyses within the context of all the identified DEGs in Group A1 [27], we will not again discuss possible functions of these genes in detail. The genes do however belong to both disparate and overlapping functional pathways. Several are involved in proteolysis with dcaf11 and uchl1 functioning in ubiquitin pathways, and senp7 and impa2 having hydrolase activity. On the other hand, tfip11 is involved in RNA processing and mrpl39 at least, is involved in mitochondrial function. Unfortunately, they were not included as target genes in previous studies on rainbow trout in which eyeing rate was a metric [29, 30]. Nevertheless, they were dysregulated in all three of our groups and therefore are the most promising as markers for poor egg quality. Unfortunately, even for these genes the regression coefficients and R² were low and therefore the expression of any one gene would not be effective at predicting quality of a single batch of eggs. As previously suggested, a suite of genes would likely be required to identify eggs of different developmental competence [20, 23]. Still more data would be required if the bases or causes of gene dysregulation, or more specific quality outcomes are to be predicted.

Conclusions

The present study confirmed DEGs for eyeing rate identified through a comparison of a small number of individuals by RNA-Seq can be extended to the broader population. However, for each of the DEGs identified there is considerable variability in transcript abundance among eggs of similar quality and low correlations between transcript abundance and eyeing rate, making it highly improbable to predict the quality of a single batch of eggs based on transcript abundance of just a few genes. The DEGs were more consistent among these two year classes of the same line than between either of these two groups and Group B which is from a different line. Greater similarity in dysregulated transcripts across year classes within the same line than among lines suggests patterns of transcriptome dysregulation may provide insight into causes of decreased viability specific to a hatchery population.

Although not as similar as among year classes within the same line, there appears to be commonality in genes that are dysregulated even among diverse populations and metrics for evaluating egg quality. Transcript abundance of over half of the 17 genes in the assay that were identified as DEGs among rainbow trout eggs of disparate quality based on a range of egg quality metrics in populations other than those in the present study, were found to be correlated with eyeing rate in at least one of our study groups, although not always in the same direction as in the previous studies. There were six genes with transcript abundance correlated with eyeing rate in all three groups, dcaf11, impa2, mrpl39_like, senp7, tfip11 and uchl1, and therefore have the greatest potential to serve as general markers for egg quality among those in our current study.

Other than the mitochondrial genes, the genes selected for the present study were primarily based on magnitude of response and statistical differences reported in previous studies. Despite this lack of focus in the genes selected, many of these genes were found to be differentially expressed in response to differences in egg quality in our other populations; and similarities and differences in expression profiles among the studies and our groups were identified and useful inferences were deduced. Future assays need to be designed for investigating specific pathways involved with egg quality parameters and designed to address dysregulation in specific hatchery populations.

Methods

Sample collection

Eggs were collected from 192 individual two-year-old broodstock rainbow trout from the May even-year selective breeding program at Troutlodge Inc. Sumner, WA, USA (Group A1); and 142 individuals from the May odd-year line (Group A2); and 325 individuals from the 2015 NCCCWA line (Group B). Samples were collected for Groups A1 and A2 one year apart, following the same procedures as described in Ma et al. [27] for Group A1. Briefly, approximately 90 unfertilized eggs from each female were frozen in liquid nitrogen for mRNA analysis and 50 unfertilized eggs were fixed in modified Davidson’s fixative [54] to be examined for overripening or other abnormalities. Sperm harvested from a single neomale was used to fertilize the remaining eggs from each female and the collected semen from each male was used for two or three crosses. Eyeing rate was evaluated at about 250 ATUs (sum of mean daily water temperature in degrees Celsius) (Fig. 1ab). About 25–60 embryos per family in Group A1 were fixed in Stockard’s solution [55] at about the 32-cell stage. Enumeration of embryos reaching each cleavage stage was used to assess viability. The data for Group A1 are reported in Ma et al. [27], and samples to evaluate early embryonic survival were not collected from Group A2. Sample collection for Group B followed similar procedures as for Groups A1 and A2 with minor exceptions. Batches of unfertilized eggs or spawns with evidence of overripe eggs or other abnormalities were not collected, but a 5-ml sample of unfertilized eggs from retained lots were also collected in neutral-buffered formalin to be examined again for elimination of batches of eggs with such traits. Two 5-ml samples of about 50 unfertilized eggs each were collected from each female and immediately frozen in liquid nitrogen for mRNA analysis. Semen derived from a single neomale sire was used to fertilize eggs from a single female, so no half-sibling crosses were generated to evaluate sperm quality. Eyeing rate was determined at about 250 ATUs (Fig. 1c). About 40 embryos from each family were collected after 19–20 h of incubation at 10 °C to be evaluated for early embryonic viability as described with Group A1 (Table 1; Additional file 1: Table S2). In addition, at approximately 10 days post fertilization 10 embryos with normal coloration from each spawn were fixed in Davidson’s fixative to provide a rough estimate of survival later in development (Table 1; Additional file 1: Table S2). Only unfertilized eggs and early cleavage stage embryos were sampled for this study and therefore no animals were euthanized.

Selection of rainbow trout females for mRNA analysis of unfertilized eggs

Egg samples were selected for mRNA analysis based on eyeing rate (Fig. 1d-f) with 130–218 individuals from each family evaluated at eyeing for Group A1, 54–207 for Group A2, and all eggs for Group B (Fig. 1a-c). Unfertilized eggs, eggs with precipitated yolk in response to egg shocking, and those with poorly developed eyes were classified as not eyed for Groups A1 and A2, whereas poorly developed eyes were not a criteria in Group B. Low-quality eggs were those with less than 20% viability at eyeing, medium-quality eggs were those with 20–80% viability, and high-quality eggs were those with above 80% viability. Eggs from the 20 females from Group A1 used in our previous RNA-seq study [27] were included in the present study. Only six females in A1 and two females in A2 yielded eying rates below 30% so eggs from all these females were included in the study. The samples from the medium and high viability families for all groups, and the low viability families for group B, were selected to provide a range of eying rates.

Assessment of early embryo development

Embryo development for females from Group B selected for mRNA analysis were examined for early embryonic development at 19–20 h post fertilization and about 10 days post fertilization. The embryos fixed in Stockard’s solution at about 19–20 h post fertilization were immersed in 0.5% methylene blue overnight before evaluation. The cell number of each of 25 embryos per family was counted or confirmed to be greater than 32 cells using a stereo microscope (Nikon SMZ660).

The ten embryos from each family with normal coloration collected at ten days post fertilization and fixed in Davidson’s fixative were examined by eye to enumerate the embryos that had reached development of the neural keel [56]. The percentage of embryos reaching this stage, was then adjusted by the estimated percent of eggs that did not have normal coloration or were white in appearance indicating yolk precipitation, based on a rough visual estimate of the eggs in the batch, and this was recorded as the streak rate (Table 1; Additional file 1: Table S2). The streak rate was meant as a hatchery tool and is a very rough estimate of early embryonic survival.

RNA isolation

Total RNA was isolated from a pool of 50 eggs per fish as described in Ma et al. [27]. Briefly, Eggs were homogenized while frozen in Tri Reagent (Sigma, St. Louis, MO) using a Qiagen Retsch MM300 TissueLyser (Retsch Inc., Haan, Germany). Total RNA was isolated using the manufacturer’s suggested protocol. The aqueous phase was separated from the organic phase using Phase Lock Gel tubes (5 PRIME, Inc., Gaithersburg, MD) and Phase Separation Reagent (Molecular Research Center, Cincinnati,OH). Isolated RNAs were purified further by a lithium chloride precipitation and then treated with DNase. Polyadenylated RNA was then collected using Oligotex mRNA Mini-Kits (Qiagen, Germantown, MD).

Transcript abundance analysis

The study used an nCounter analysis data systems assay comprised of 65 transcripts previously identified as being associated with measures of egg quality in rainbow trout, along with four reference genes (Additional file 1: Table S1). Annotated sequence data for each gene was submitted to Nanostrings Technologies (Seattle, WA) for CodeSet design. Nanostring Technologies designed custom CodeSets for nCounter® analysis using the submitted target sequences and rainbow trout transcriptome sequence data to avoid non-specific hybridization. CodeSet details are provided in Additional file 1: Table S1. Transcript measurement was conducted following the nCounter analysis system workflow protocols. Ten ng of mRNA was used for each sample for A1, 7.6 ng for A2, and 6.15 ng for B. The raw data before normalization are provided in Additional file 1: Tables S4A-D and normalized data used for transcript abundance analysis are presented in Additional file 1: Tables S5A-C.

Normalization of nCounter transcript abundance data

The nSolver Analysis Software (V2.0) was used for normalization of transcript abundance data. The average geometric mean of positive spike-in RNA controls was used across all samples to normalize for technical aspects of assay performance. Although four reference genes were included in the assay, both ef1a and ppia were below detection limits and b-actin was not stable in Group A1 samples. Therefore, the geometric mean of g6pd and b-actin were used for normalization in Groups A2 and B, and g6pd alone for Group A1. The assay included both mitochondrial and nuclear genes. The mean raw reads per gene was 6107, however, the mean raw reads for the five mitochondrial genes was 76,186 whereas the mean raw reads for the nuclear genes was 632 (Additional file 1: Table S4D). The high reads for the mitochondrial genes overwhelmed the nuclear genes resulting in 12 genes in addition to two of the reference genes being less than the geometric mean plus two standard deviations of the negative controls in at least one of the groups (Additional file 1: Table S4D). The average CV among the three studies for the raw reads for the individual genes was 18.6% (Additional file 1: Table S4D).

Statistical analysis

Statistical analyses were conducted separately for the three studies and for nuclear and mitochondrial transcripts. The transcripts selected for the assay were based primarily on studies in which low egg quality was classified as having eying rates below 20% and high egg quality as having eyeing rates above 80% [27, 29, 30]. Families were therefore classified as having low-quality eggs, 0–20% eyeing, medium-quality eggs, 20–80% eyeing, and high-quality eggs, 80–100% eyeing in the present study and mean and standard errors of the mean (SEM) are presented (Tables 2, 3, 4; Fig. 2 a-c). Unfortunately, Group A1 had only six and Group A2 had only two samples with eyeing rates below 30%, hampering the ability to use multivariate analysis, and therefore regression and analysis of variance were used to test for correlations and measure the effect of gene transcript abundance on the egg quality phenotype eyeing rate.

Prior to conducting statistical data analysis, the datasets were transformed as needed to meet normality of distribution and equal variance requirements of the statistical data analyses. We performed log10, arcsin and square-root data transformation, and assessed the fit of the transformed data to normal distribution using statistical tests available in the Procedure Univariate from SAS software [57]. Then, to quantify the correlation or impact of gene transcript abundance on the egg quality phenotype eyeing rate, we performed two types of statistical data analysis: First, we estimated the regression coefficient of the continuous phenotype eyeing rate on gene transcript abundance using the Procedure Regression from the SAS software [57]. Second, we performed analysis of variance to determine whether each discretized eyeing rate phenotype (i.e., low, medium, high) was associated with a distinct level of gene transcript abundance using the Procedure GLM from the SAS software [57]. For both analyses, the threshold significance level was set at the nominal P-value of P < 0.05. Although data for all genes were analyzed, those with average raw read levels less than that of the geometric mean plus two standard deviations of the negative controls for that group, were considered unreliable and therefore below detection. The geometric means plus two standard deviations for the negative control raw reads were 40.3, 38.3 and 30.6 for Groups A1, A2, and B respectively (Additional file 1: Tables S4A-C).