Maternal 3’UTRs: from egg to onset of zygotic transcription in Atlantic cod

Methodological considerations; number of sequences threshold and external spiking

Transcriptome sequencing offers the possibility to discover novel expression patterns, high abundance transcripts, sequence motifs and splice variants [6]. In addition transcriptome sequencing measures abundance at a wider scale, as this technique does not get oversaturated. However challenges regarding quality do exist, it is for example crucial for the correct annotation to have a high-quality reference genome. Another aspect of quality is how reliable the expression patterns discovered in libraries really are. Highly expressed genes represent processes that dominate in the given developmental stage, and discovery and measurements of such gene expressions are more likely to be accurate than that of genes that are weakly expressed. Singletons (n = 1) are more likely to represent random sequencing or sequnecing errors than transcripts that are sequenced numerous of times for one gene. Likewise, in microarray studies, low number of transcripts causes a large variation in gene expression results [23]. This was also confirmed in this study as we revealed that for genes with few sequences (total n = 1-5) detected in the libraries, the expression profile from oocyte until gastrula could not be reflected in expression profiles measured with quantitative PCR (qPCR) for the same genes (Figure 2). For genes with a higher number (n = 12-61) of sequences, their expression profiles corresponded with the ones measured by qPCR (Figure 2). We therefore suggest for further transcriptomic studies to filter data based on analysis of number of transcripts to obtain reliable expression profiles for downstream analyses. A threshold of 9 or more total sequences for each gene was selected to be applied for our dataset, to be able to include only the most reliable data for analysis. Threshold 9 was chosen because it is in the middle of 5 and 12, and using this threshold did not considerably affect number of genes to be excluded for downstream analyses compared to when applying threshold 5. Given that only 1 EST was sequenced for approximately 50% of the genes in our dataset (data not shown), even a very low cutoff would exclude a high number of genes from the analysis. However some important information may be hidden within genes with low expression, and will therefore be lost by applying an EST number threshold. For example, the zebrafish oocyte possesses a considerable richness of maternally expressed genes of which many may be of importance for embryogenesis [6]. When the current experiment was carried out the cod genome was not yet sequenced. However, using the current genome [24] together with new sequencing technology could significantly increase number of reads analyzed in this study. For example, Aanes recently studied the zebrafish embryonic transcriptome by applying RNA sequencing [6], a method which provides a considerably more precise measurement of transcript abundance due to high number of sequence reads [20–22]. If we were to perform this study again with RNA sequencing, more sequence reads would be obtained and more genes could be included in the analysis as they would be represented by more than only 1 or a few sequenced transcripts.

Also if candidate gene expressions from cDNA libraries are to be verified by qPCR, stable reference genes are essential. When it comes to early embryo tissue like in this study, such stable references are difficult to obtain. For example, neither β-actin elongation factor 1- α (ef1α) 18S, actin related protein (arp) nor ubiquitin (ubi) showed stable reference gene expression in stages assayed for gene expression [Additional file 1, [25, 26]. For this reason we chose not to use any of the reference genes tested, but instead apply a technique (external spiking) where external RNA (rabbit alpha hemoglobin (hba)) is added in equal amounts to each RNA sample [27]. The considerably more stable expression pattern of hba can be seen in [Additional file 1. We therefore suggest that external spiking may serve as a good alternative to the traditional reference genes previously applied in qPCR, especially when measuring gene expressions in early embryo tissue.

Overall composition of cDNA libraries from oocyte, 1-2 cell, blastula and gastrula stage in Atlantic cod: amount of transcripts and genes, and mitochondrial contribution

In order to get an impression of the approximate timing of MZT (and therefore the shift from degrading to activating transcripts) we measured gene expression profiles of candidate genes known from other species to be maternal, zygotic or involved in events taking place at the initiation of MZT [Additional file 2. Our results together with the findings of Drivenes [19] indicate that the MZT takes place between blastula and gastrula in Atlantic cod, which is also true for zebrafish [12]. Together the cDNA libraries from oocyte, 1-2 cell, blastula and gastrula stage produced 26162 3’EST sequences, which have been added to GenBank; 6511, 6607, 6506 and 6538 sequences were obtained from oocyte, 1-2 cell, blastula and gastrula libraries, respectively. 22229 sequences (5538, 5875, 5238 and 5578 in oocyte, 1-2 cell, blastula and gastrula, respectively) were mapped to the Atlantic cod genome [24]. 17025 of these sequences hit cod genes with the rules given in material and methods section; they were predicted to encode 2611 genes, of which 2051 (78.6%) were Ensembl-annotated. 23 genes were pseudogenes, and were therefore removed from most of the analysis. Furthermore as described in methodological considerations, sequences with n ≤ 8 were removed from most of the analysis, leaving us with 12764 sequences predicted to be encoded by 257 genes. Compared to similar analyses [20–22], 257 genes is a low number; however these studies have included most transcripts in contrast to the current study. As mentioned, we showed no correlation between qPCR and low copy number of ESTs (N < 9, see Figure 2 and discussion first section; methodological considerations), and highly expressed genes (N > 9) are possible to classify as maternal, zygotic, degrading and activating. EST copy number filtering was a necessity to be able to study stability and complexity of dominating transcripts in early embryos of cod.

Strikingly, the overall composition of abundant transcripts obtained from maternal stages (oocyte, 1-2 cell and blastula) contained 45.9-50.6% mitochondrial transcripts, while at the zygotic stage (gastrula) only 12.2% of the transcripts were mitochondrial (Figure 3). Less than 5.1% of all the expressed genes in each developmental stage encoded these numerous mitochondrial transcripts (data not shown). This trend was also evident by the highly expressed mitochondrial genes of Table 1, showing the 20 most abundantly expressed genes when all transcripts from all libraries were included. Likewise it has been reported a high level of mitochondrial transcripts at the 2-cell stage in halibut [22] and in the oocyte of Senegalese sole [28]. High level of mitochondrial transcripts in mature oocytes may be linked to high metabolic need to complete synchronous continuous cell-divisions at early developmental stages, whereas at later stages the energy metabolism is more localized to target areas in specific development. It has even been suggested that a major part of maternal RNAs in egg-laying organisms lack specificity and only provides supplemental nutrition [29]. However a proportion of the maternal RNAs are essential for normal development in fish, as several studies report abnormal embryo development when mRNA translation is inhibited [30]. Also several processes in zebrafish morphogenesis, most notably brain formation, fail to occur normally when normal degradation of maternal RNA is disturbed [31].

Nuclear transcripts and expressed genes: distribution and gene ontologies

To facilitate analysis of expressed genes from the genome, mitochondrial transcripts (4786, encoded by 10 genes) were excluded. What remained were the nuclear transcripts (7978, encoded by 247 genes). A significant feature of the gastrula nuclear transcriptome was a high presence of ribosomal transcripts (50.6%) compared to the earlier stages (2.5–4.6%, Figure 4a). This increase in transcript abundance was not a consequence from many additional genes being expressed in the gastrula stage, as the number of genes encoding the ribosomal transcripts only increased from 6.7–13.1% in the maternal libraries to 19.7% in the gastrula library (Figure 4b). However, some of the ribosomal transcripts were uniquely found in the gastrula library (Table 2). A similar pattern of ribosomal transcripts has previously been found in zebrafish [17] and halibut [22]. The high presence of ribosomal transcripts in gastrula is most likely due to major onset of transcription since ribosomal proteins are required for the production of ribosomes at the site for transcription.

We also checked how the nuclear transcripts were distributed when including all nuclear genes (2575) without any threshold for number of transcripts. In this dataset a large proportion (56.5%) of the genes had only 1 transcript. Consequently, the genes expressed exclusively in maternal and gastrula stages showed an increase in their proportions (12.1–31.9% and 6.9–13%, respectively). The proportion of genes expressed in all developmental stages decreased markedly from 51.8 to 6.5% [Additional file 3]. It is therefore clear that the use of a threshold for number of transcripts makes a difference for describing gene expression patterns.

Functional profiles of the most common nuclear transcripts in oocyte, 1-2 cell, blastula and gastrula (obtained from the use of gene ontology (GO)-terms) revealed that both within cellular component, molecular function and biological process many of the transcripts did not change much in their abundance over time (data not shown). However some differences could be observed between maternal and zygotic stages [Additional file 4. An increase in proportion of transcripts that were part of macromolecular complexes from maternal to zygotic developmental stage possibly reflects specification of cell types which requires more complex molecules as the embryo enters the zygotic stage. At the molecular function level, a pronounced increase of transcript amount from maternal to zygotic stage was evident for structural molecule activity, reflecting a shift towards specification of cell function. Several biological processes including cellular component organization or biogenesis, multicellular organismal process and developmental process showed a marked shift between maternal and zygotic stage, with increase in transcript proportion in gastrula. Taken together, these changes reflect more specialized and complex processes developing after MZT in cod, as would be expected since gastrula is a developmental stage where intricate processes like germ layer formation takes place (reviewed in [33]).

3’UTR structure: isoforms, length and motifs

Based on manual annotation, 3’UTR length was predicted for all sequences with stop codon and poly-A tail (in total 191 genes). To be able to study the gene expression dynamics over time from a maternal to a zygotic embryo, genes (when possible) were strictly classified as having degrading, stable or activating transcripts. In order to be termed degrading, minimum 95% of the transcripts had to be located in maternal stages (oocyte until blastula). In order to be termed activating, minimum 95% of the transcripts had to be present in blastula and gastrula. Transcripts present in all 4 stages, with a distribution of minimum 10% and maximum 40% in one stage were termed stable. Applying this classification gave 35, 26 and 39 genes with degrading, stable and activating expression profiles, respectively. 91 genes could not be classified according to these specifications. To ensure reliable results, the following analyses on 3’UTR structure were also performed when applying less strict specifications [Additional file 5].

Totally 73 nuclear genes (29.6%) expressed more than one 3’UTR isoform. Among genes with degrading and stable transcripts a relatively high amount (45.7 and 34.6%, respecively) displayed 3’UTR isoforms (Figure 5). Alternatively spliced transcripts are common in early embryos in other species, for example in zebrafish where it has been shown 50–60% alternative splicing of transcripts [6]. Also in mouse embryos a significant proportion of alternatively spliced transcripts have been found [34], and in this context it was suggested that alternative splicing is of importance in the regulation of development. Activating genes of cod had a far less degree of 3’UTR isoforms (15.4%, Figure 5). This may be an indication that alternative splicing is a less important regulatory mechanism for genes that are more active in early zygotic stages, than for genes that are more active in maternal stages.

Mean 3’UTR length (when including all isoforms) in degrading (n = 48), stable (n = 35) and activating (n = 45) transcripts was 187.1, 208.8 and 76 bp, respectively (Figure 6). Activating transcripts had significantly (p < 0.0001) shorter 3’UTR length than degrading and stable ones. A higher proportion of long maternal 3’UTRs have previously been found across the animal kingdom [29] associated with lower turnover-rate of RNAs [35]. Our findings suggest that also in cod, transcripts from maternally active genes have longer 3’ UTRs. The shorter 3’ UTRs of the activating transcripts may also be explained by the fact that many of them were ribosomal (data not shown); such transcripts have been shown to contain short 3’UTRs [36]. Excluding ribosomal genes would lead to a small dataset. Consequently, we applied less strict classification criteria for degrading, stable and activating transcripts and found that also non-ribosomal activating transcripts have significantly shorter 3’UTRs [Additional file 5.

Although not statistically proved in this study, stable transcripts appeared to have the longest 3’UTRs (Figure 6). This corresponds to the suggestion that more stable RNAs have longer UTRs [37]. A strong exponential correlation has been found between 3’UTR length and morphological complexity (number of cell types in the organism) [38]. Increased regulatory complexity in long 3’UTRs may be a result of more miRNA binding sites. Furthermore, additional polyadenylation signals may produce transcript isoforms with several combinations of miRNA binding sites, which adds even more complexity [39]. It was observed an enrichment of 3’UTR isoforms in degrading and stable transcripts (abundant in the maternal stages) compared to activating ones (abundant in zygotic stage). This result suggests a more complex post-transcriptional regulation in the maternal stages since transcriptional regulation cannot be applied to control developmental processes prior to MZT.

The apparent complex features of transcripts that were abundant in early developmental stages of cod were further investigated by searching for motifs in the 3’UTRs of degrading, stable and activating transcripts. Using Multiple Em for Motif Elicitation (MEME) [40] the presence of the polyadenylation motif AATAAA was indicated. Manual inspection of the sequences revealed that AATAAA (the exact motif) was present at an increasing manner from degrading (68.8%), through stable (82.9%) and to activating (95.6%) transcripts (Figure 7a). AATAAA is essential for cytoplasmic polyadenylation to take place [41, 42]; that is, cleavage followed by adding of a poly-A tail to the 3’ end of the mRNA which then activates translation. Therefore polyadenylation is a key process in regulation of translation (for recent reviews see [43, 44]), and thus serves as an important step in post-transcriptional regulation in maternal stage embryos. Also, it was recently found that cytoplasmic polyadenylation of maternal RNAs prior to MBT is essential for normal morphogenesis post-MBT in zebrafish [6]. It is known that AATAAA is found in nearly all known mRNAs, and that it is highly conserved. Our finding that 31.2 and 17.1% of degrading and stable transcripts (with known 3’UTR length) did not contain AATAAA may be explained by the possibility of transcripts to contain alternative polyadenylation motifs, where AATAAA have been slightly modified (reviewed in [43]). Therefore, we investigated the presence of 11 other putative alternative polyadenylation motifs [45]; for each sequence we observed which motifs were present, and if there was more than one motif (Figure 7). 3 motifs (ATTAAA, AATATA and AATAGA) were found exclusively in degrading and stable transcripts, and AGTAAA was present only in stable transcripts. None of the motifs were unique for activating transcripts, and most of the motifs were present at the lowest rate in this group of transcripts (Figure 7a). Also, degrading and especially stable transcripts had a higher prevalence (41.7 and 54.3%, respectively) of multiple motifs (more than one of the motifs studied) than activating ones (31.1%, Figure 7b). A higher degree (diversity and quantity) of putative alternative polyadenylation motifs in maternal degrading and stable transcripts suggests that these transcripts may apply alternative polyadenylation more frequently than activating transcripts. Polyadenylation at different sites produces different variants of the 3’UTR, which contributes to a more complex transcriptome [46].

Egg and embryo collection

Eggs and sperm from Atlantic cod were obtained from Parisvatnet Research station, Institute of Marine Research, Øygarden, Norway in April 2009. Eggs were manually stripped from one female and sperm was collected from one male. Before freezing of unfertilized eggs, egg fluid was carefully aspirated and subsequently eggs were flash frozen liquid nitrogen. Remaining eggs were fertilized with the collected sperm and incubated at 7°C until they reached 1-2 cells, blastula and gastrula stages (Figure 1), embryos were staged according to Hall [47]. At sampling, staged embryos were collected and excess water was aspirated off before embryos were flash-frozen on liquid nitrogen. All collected samples were stored at −80°C until RNA isolation was initiated. All egg and embryo samplings were done according to guidelines approved by the Norwegian Animal Research Authority (NARA).

cDNA library construction

From eggs and each developmental stage, total RNA was isolated using iPrep™ Trizol® Plus RNA Kit (Invitrogen, Carlsbad, CA, USA). The quality of RNA showed RIN values above 9.5 for all samples used (Bioanalyzer 2100, Agilent, Santa Clara CA, United States). cDNA libraries were constructed from 1.0 μg of RNA from each sample using Creator SMART cDNA Library Construction Kit (Clontech Laboratories, Inc., Mountain View, CA, USA). The cDNA was size-selected for 500bp according to Evsikov [21]. Transformation of the ligation mix was performed using the DH5αMAX competent cells (Invitrogen, Carlsbad, CA, USA). The titers of libraries used were for oocytes 7,7x10⁵, 1-2 cell stage 4,7x10⁵, blastula 1,5x10⁶ and gastrula 1,4x10⁶ clones per μl of ligation mix. The clones were grown on robot plates (Nalgene Nunc International, Rochester, NY, USA) filled with Luria-Bertani (LB) agar and chloramphenicol (30mg/L). Positive clones were randomly selected by automated picking robots and transfered to 384-well plates containing HMFM freezing media. Culturing was performed for 16h at 37°C without agitation. For isolation and purification of plasmid DNA 384-deepwell plates containing 190μl 2YT media and 30 mg/L chloramphenicol were inoculated from the freezing stocks. Clones were grown for 16-18h at 1100 rpm and 37°C. Isolation and purification of plasmids was done by an automated system, applying alkaline lysis followed by size selective DNA precipitation in PEG/2-Propanol [48]. Quantified plasmids (UV photometry) were sequenced using T7 primers following the Big-Dye version 3.1 protocol and capillary gelelectrophoresis by an ABI 3730xl system (Applied Biosystems Inc., Foster City, Ca, USA) at Max Planck Institute for Molecular Genetics, Berlin, Germany.

Bioinformatical analysis of raw data

Sequence trace files from each library were processed employing the Phred base-calling algorithm [49, 50]. Quality trimming and removal of contaminating vector sequences were performed using cross match and RBR [51] with the Uinivec filtering database from the National Center for Biotechnology Information (NCBI). Clustering and assembly were carried out using the TIGR Gene Indices Clustering tools (TGICL) [52]. Sequences for the libraries were added to GenBank; oocyte library [GenBank: JK777860-JK784370], 1-2 cell library [GenBank: JK771253-JK777859], blastula library [GenBank: JK758209-JK764714] and gastrula library [GenBank: JK764715-JK771252], respectively.

Clustered sequences were assembled into contigs (Fasta sequences are given in [Additional file 6) which were further annotated with Ensembl predicted genes [53] by mapping the contigs to the Atlantic cod genome [24] with a threshold above 75% identity BLAT hit against a genome scaffold in a gene or within 2000 bp downstream of a gene. Sequences that passed the following rules were included in further analysis: 1) above 85% identity BLAT hit against a genome mapped EST-contig, or 2) above 75% identity BLAT hit against a genome scaffold in a gene or within 2000 bp downstream of a gene. The mapping of the contigs can be seen in [Additional file 7.

Functional profiles of transcripts were constructed by applying GO-terms to group the transcripts [54]. The Ensembl genes were GO-annotated using BLAST against SwissProt and UniRef90 and annotation files from [55]. BLAST and GO files for SwissProt and UniRef90 that are constructed for Ensembl-annotated genome assembly can be seen in [Additional file 8 and [Additional file 9, respectively.

Quantitative PCR

Statistical analyses

All data to be analyzed for significant differences were initially tested for Gaussian distribution by applying a D’Agostino-Pearson normality (omnibus K2) test. All mRNA levels quantified by qPCR had too few replicates (n = 3) for each developmental stage in order for the test to reveal normality or not. Therefore a non-parametric Kruskal-Wallis test was applied to reveal if the genes were differentially expressed in one or more of the developmental stages. Dunn’s multiple comparison post-test was used to detect at which developmental stages there were differences. Furthermore there was no evidence for Gaussian distribution in the datasets for 3’ UTR length analysis, including when they were converted to log values. Therefore a non-parametric Kruskal-Wallis test was also applied here to calculate if there were any differences in UTR length between different groups of transcripts. Dunn’s post test was applied to reveal which groups differed from each other. The statistical tests were performed using GraphPad Prism 5.04 (GraphPad Software Inc., La Jolla, CA92037, USA). A p-value of ≤0.05 indicates significant difference.