Zebrafish mRNA sequencing deciphers novelties in transcriptome dynamics during maternal to zygotic transition

Abstract

Maternally deposited mRNAs direct early development before the initiation of zygotic transcription during mid-blastula transition (MBT). To study mechanisms regulating this developmental event in zebrafish, we applied mRNA deep sequencing technology and generated comprehensive information and valuable resources on transcriptome dynamics during early embryonic (egg to early gastrulation) stages. Genome-wide transcriptome analysis documented at least 8000 maternal genes and identified the earliest cohort of zygotic transcripts. We determined expression levels of maternal and zygotic transcripts with the highest resolution possible using mRNA-seq and clustered them based on their expression pattern. We unravel delayed polyadenylation in a large cohort of maternal transcripts prior to the MBT for the first time in zebrafish. Blocking polyadenylation of these transcripts confirms their role in regulating development from the MBT onward. Our study also identified a large number of novel transcribed regions in annotated and unannotated regions of the genome, which will facilitate reannotation of the zebrafish genome. We also identified splice variants with an estimated frequency of 50%-60%. Taken together, our data constitute a useful genomic information and valuable transcriptome resource for gene discovery and for understanding the mechanisms of early embryogenesis in zebrafish.

Figure 1.

Overview of mRNA-seq data and mapping to the zebrafish genome. (A) Zebrafish embryonic developmental stages analyzed in our study. Libraries were generated from eggs, 1-cell, 16-cell, 128-cell, 3.5-hpf, and 5.3-hpf embryos. The graph depicts general expression patterns during these developmental stages. (B) Total number of mRNA-seq reads mapped in each developmental-stage library. Differences in the number of generated reads were observed, but the proportions of mapped reads were similar between all stages. (Blue bars) The number of reads generated from sequencing; (red bars) the number of reads mapped to the genome; (green bars) the number of reads mapped to Ensembl-annotated genes. (C–F) Mapping of mRNA-seq reads as viewed in the UCSC Genome Browser. Horizontal red bars at the bottom of each panel represent annotated exons (tallest boxes), UTRs (half-sized boxes), and introns (lines with arrowheads). (Arrowheads on bars) Orientation of transcript from 5′ to 3′. Mapped reads are represented by red histograms. (C) Strand-specific mapping in a region of chromosome 5, showing reads mapping to mobkl1aa in the negative, and grsf1 in the positive strands. (D) Reads mapping to the htt gene, containing 67 exons. Note the sharp exon–intron boundary of the mapping. (E) Mapping of exon–exon junction spanning reads (red horizontal lines) to klf4 locus, allowing the derivation of splice junctions. (F) Mapping of reads to regions beyond annotated 3′ (zic3) or 5′ UTRs (foxa3), indicated by arrows. Di-tags (blue horizontal bars at top of panel) provide support for this mapping.

Figure 2.

Expression clusters of early developmental genes. (A) Heatmap showing distinct expression profiles of different clusters at each developmental stage. Values were scaled for each cluster; color intensity represents expression level relative to its own cluster's average. (B–D) Graphs showing mean expression of all genes in each cluster from the egg to the 5.3-hpf stage. (B) Two clusters of the maternal supercluster, showing two different degradation patterns. (C) Two clusters of the pre-MBT supercluster, both characterized by a significant increase at the 16-cell stage, with either decrease (pre-MBT1) or a continued increase after 3.5 hpf (pre-MBT2). (D) Zygotic supercluster consisting of genes activated at 3.5 hpf (MBT) and 5.3 hpf (post-MBT), or maternal–zygotic cluster genes that showed stable accumulation of transcripts during pre-MBT and with a sudden increase at MBT and post-MBT.

Figure 3.

Validation of mRNA-seq expression clusters. (A–C) RT-qPCR of the maternal gene cldng and the zygotic genes id1 and hspb1 show similar expression patterns to those detected by mRNA-seq. WISH confirmed the expression pattern of these genes as observed in mRNA-seq. (D) RT-qPCR of the pre-MBT gene xbp1 showed a different expression profile from that of mRNA-seq when random primers (r.p.) were used, while it showed good correlation when oligo d(T) primers were used. (E) WISH of xbp1 demonstrates the presence of maternal transcript in the egg and the 1-cell stage and confirmed r.p. RT-qPCR data. (F) Measurement of xbp1 poly(A) tail length by RT-PCR in the egg and 1-cell, 16-cell, and 5.3-hpf stages. Total RNA was G/I-tailed and amplified with the xbp1-specific forward primer (F) and a universal oligo d(C) reverse primer (U). A control RT-PCR with xbp1-specific F and R primers gave a 180-bp fragment (S). The increasing poly(A) tail appears as a smear that increases gradually from the 1-cell stage to the 5.3-hpf stage, indicating the presence of an extended poly(A) tail.

Figure 4.

Inhibition of cytoplasmic polyadenylation by 3′-dA. (A) Schematic of treatments at various developmental time frames. Horizontal arrows represent exposure to 3′-dA. (B–D) Control embryos showing normal development. (Red arrows) The level of marginal cells, signifying the progress of epiboly. Note the presence of embryonic shield (S) at the 6-hpf stage. (E–G) Embryos treated with 3′-dA from the 1-cell to 3.5-hpf stage, showing a severe delay in development. (H) U-rich motif found in genes of the pre-MBT supercluster using the MEME algorithm. (I) Distance between Hex and CPE in different groups of genes. The boxplot shows the distribution of the shortest distance between any Hex and CPE pair in the 3′ UTR of genes in maternal–zygotic and pre-MBT clusters, as well as in a randomly chosen group of genes (random). Notice the longer distance between the two elements in the pre-MBT group compared to the other two groups.

Figure 5.

Discovery of NTRs. (A) Genome-wide distribution of NTRs plotted at their mapped chromosomal positions. The color code depicts at which stage an NTR was detected. (B) Novel exons of ENSDARG00000035540, some of which correspond to exons of its human homolog, BRWD1. Splice junctions were also identified among these novel exons. RT-PCR performed using reverse primer at annotated exon 2 (R) and forward primers at each of the 12 NTRs (inset) showed that at least 11 of these NTRs are expressed and are part of the transcription unit. (C) Novel exon in foxa found 5′ of annotated transcription start site. Note the presence of the 5′ annotated exon at 3.5 hpf and 5.3 hpf (white arrowheads). Di-tag (blue horizontal lines) provides support for an isoform that includes the novel exon. Junction mapping at 5.3 hpf (green horizontal lines) indicated the presence of two distinct transcript isoforms. RT-PCR using a reverse primer (R) designed on the 3′ distal exon and forward primers on either 5′ annotated (*) or NTR (+) confirmed the existence of a new 5′ distal exon (inset). The 5′ annotated exon (*) is skipped in the maternal mRNAs (egg to 128-cell stage).

Figure 6.

Analysis of splice variants. (A) The number of alternatively spliced transcripts (blue bars) and alternatively spliced exons (red bars) at each stage. Both were stable through all stages. (B) Location of alternatively spliced exons in all transcripts. Most spliced exons were located between exons, but with a substantial fraction corresponding to distal 5′ exons. (C) Splice variants in msna at 5.3 hpf. Junction mapping shows skipping of exons 1 to 8, 2 to 8, and 3 to 11 (blue arrows). Di-tags (red bars with arrowheads) support the presence of several splice isoforms.