Distinct gene expression signature sets of ERK1 and ERK2 knockdown embryos
A morpholino knockdown approach was used to block translation of either ERK1 or ERK2 by injection of 0.4 mM (= 3.4 ng/embryo) morpholinos (MO) targeting ERK1 (ERK1MO) or ERK2 (ERK2MO). The knockdown embryos, also referred to as morphants, showed severe phenotypes after depletion of ERK2. These embryos did not enter epiboly at 4.5 hpf and the blastula cells remained on top of the yolk, preventing further development of the embryo (Fig. 1C). In addition, ERK2 knockdown induces disorganization of the margin. (Fig. 1F). Wild type embryos reached 30% epiboly at this time (Fig. 1A,D) [13]. In contrast, ERK1 morphants did not show any obvious phenotypes at this point yet and had entered epiboly (Fig. 1B,E). However, ERK1 morphants did show strong phenotypes at later stages in embryogenesis (manuscript in preparation and additional file 2).
The specificity of ERK1 and ERK2 knockdown phenotypes was rescued by co-injection of synthetic mRNA (data not shown, manuscript in preparation), western blot analysis (Fig. 1J) and immuno-localization in wildtype, ERK1MO and ERK2MO injected embryos at 4.5 hpf (Fig. 1G–I) using a phospho-specific ERK antibody (dpERK). ERK1 morphants (Fig. 1H) show similar levels of dpERK staining at the dorsal margin compared to wild type embryos (Fig. 1G), but the dpERK signal in ERK1 morphants is reduced at the ventral half of the margin. ERK2MO injected embryos hardly show any dpERK staining and the active ERK signal is depleted from the marginal ring in these embryos (Fig. 1I). Quantification of a western blot analysis of zebrafish wild type, ERK1MO and ERK2MO injected embryos, probes with a global ERK antibody (Santa Cruz Biotechnology), recognizes both zebrafish ERK1 (45.7 kD) and ERK2 (43.3 kD) protein clearly shows the specific knockdown of either ERK1 or ERK2 by the corresponding morpholino (Fig. 1J).
Addition of different MEK specific inhibitors (U0126 or PD98059, Cell Signaling technologies), did not result in the same phenotypes as obtained by the ERK2MO mediated knockdown. The inhibiting effects of these drugs were confirmed by Western blot analysis, but apparently these effects were not efficient enough to block epiboly (data not shown). Because it is not possible to address the specific functions of either ERK1 or ERK2 using these chemical inhibitors, we did not proceed with these experiments.
These data prove the functionality of the morpholino-mediated knockdown of either ERK1 or ERK2. In addition, the severe phenotype of ERK2 morphants indicate defects in crucial early developmental processes and most likely affects the expression levels of a larger number of genes than knockdown of ERK1.
Distinct ERK-knockdown gene expression profiles in time
To identify specific gene pools affected by the knockdown of ERK1 or ERK2, and to identify possible downstream targets, microarray based transcriptome analysis was performed using Agilent zebrafish microarrays. As a control for aspecific morpholino effects, a standard control morpholino (GeneTools Philomath, OR, USA) was injected in the same concentration. This did not result in any phenotypes during zebrafish development. The RNA from these standard control MO injected embryos was used as a reference to compare the transcriptomes of both ERK1MO and ERK2MO injected embryos. We annotated the Agilent 22K-zebrafish microarray chip by BLAST searches with all oligonucleotide sequences in the zebrafish genome. From the complete number of 21495 oligonucleotides from the Agilent 22K zebrafish chip, 16675 oligonucleotides were assigned a Unigene ID (build #105). The phenotypic effect of ERK2 depletion was observed at 30% epiboly indicating an altered gene expression profile at earlier stages. Therefore we analyzed the gene-expression profile of ERK2 morphants at more time points than ERK1 morphants (Fig. 2A,B). We obtained gene expression profiles for ERK1 morphants at 4.5 hpf and 8 hpf, and for ERK2 morphants at 3.5 hpf, 4.5 hpf, 6 hpf and 8 hpf (equivalent to oblong-stage, 30% epiboly, shield-stage and 80% epiboly time points), as shown in figure 2C and 2D.
Comparison of the gene expression profiles of ERK1 and ERK2 morphants at various stages showed a larger number of Unigene identifiers with significant changes (p < 10-5) in ERK2 compare to ERK1 morphants in time, as illustrated in a Venn-diagram (Fig. 2A and 2B). These Venn diagrams also show that 207 genes are affected in expression at both 30% and 80% epiboly in ERK1 morphants, whereas in the ERK2 morphants time-series, we find 186 genes to be significantly changed in expression in all time points. At 30% and 80% epiboly the numbers of genes with an altered expression was larger and with a higher fold of change for ERK2 than for ERK1 morphants (Fig. 2C and 2D). This is in agreement with the phenotype of ERK2 knockdown embryos that indicates a more prominent role for ERK2 in early development (Fig. 1). The effect of ERK1 knockdown becomes more pronounced at 80% epiboly (Fig. 2A,C). This indicates that ERK1 may become relatively more important at later developmental stages.
Comparing the effect of ERK1 and ERK2 knockdown, we found distinct gene expression signature sets during embryonic development. (Fig. 2E,F and 2G). In addition to commonly affected genes (Fig. 2G, target gene pool C), we found distinct genes that were specifically regulated by either knockdown of ERK1 (Fig. 2G, target gene pool A); 198 vs. 281 up-regulated, 109 vs. 345 down-regulated at 30 or 80% epiboly respectively) or knockdown of ERK2 (Fig. 2G, target gene pool B); 1311 vs. 1228 up-regulated, 934 vs. 786 down-regulated), or genes which were regulated in an anti-correlated manner: 32 genes (30% epiboly) or 106 genes (8 hpf, equivalent to 80% epiboly time point) were up-regulated by knockdown of ERK1 whereas they were down-regulated by knockdown of ERK2 (Fig. 2G, target gene pool D1; anti-correlated gene-pool 1) and 16 genes (30% epiboly) or 204 genes (8 hpf, equivalent to 80% epiboly time point) were down-regulated by knockdown of ERK1 whereas they were up-regulated by knockdown of ERK2 (Fig. 2G, target gene pool D2; anti-correlated gene-pool 2). These results confirm that ERK1 and ERK2 MAPK are key regulators of distinct gene signature sets during embryonic development. This is supported even when comparing multiple gene expression profiles from different developmental time-points (Fig. 2).
Because we observed a strong activated ERK signal in the margin at the onset of epiboly, we compared the expression levels of a selection of genes that are described to be expressed in the margin at the onset of epiboly in time. To do so, a gene-expression trend-line of the selected margin genes for ERK1 and ERK2 morphants was constructed (Fig. 3A,B). Most of the selected genes did not give a significant difference in time upon ERK1 knockdown, suggesting different developmental functions for ERK1, whereas in ERK2 morphants the expression levels of most of the selected 'margin'-genes was affected. A common trend in the expression-levels of the selected genes was observed upon ERK2 depletion, as most genes showed stabilization in their expression levels between 30% epiboly and shield stage, and even recovery between 6 to 8 hpf. Despite this, the presumptive blastula cells remained on top of ERK2 morphants. This indicates that the obtained gene expression profiles of later stages of ERK2 morphants (6hpf and 8hpf) are the results of a prolonged epiboly arrest, most likely due to multiple secondary developmental defects (Fig. 2 and 3). To analyze possible apoptotic effects in the ERK2 morphants, we also made similar trend lines with a selection of genes that are associated with apoptotic responses (Fig. 3C,D). The apoptotic responses by the ERK1MO treatment are minimal and also ERK2MO injected embryos do not show obvious responses in the earlier stages (3.5–6hpf). However, the apoptosis responsive genes casp8 and casp3 revealed an increased expression at 8hpf in ERK2 morphants. These combined results were the rational for limiting the further comparisons of the effects of ERK1 and ERK2 knockdown at 30% epiboly.
The identified gene-sets of correlated and anti-correlated regulated genes by knockdown of either ERK1 or ERK2 at 30% epiboly are listed and annotated [see Additional file 1, tables 1 and 4]. To identify the ERK1MO and ERK2MO specific genes, we focused on the genes that were most significantly affected. Therefore we used the following criteria: the absolute fold change must be at least 1.5 in each independent replicate and the common p-value provided by the error-model taking into account all hybridizations must be smaller than 10-5. The genes that were only found in either ERK1MO or ERK2MO gene-pools were manually annotated and assigned gene designations [see Additional file 1, tables 5 and 6].
Quantitative real time PCR analyses confirm the different ERK1- and ERK2- knockdown gene expression profiles
To confirm the results of the microarrays experiments, quantitative reverse transcriptase PCR analysis was performed on seven regulated genes at 4.5 hpf (30% epiboly) that were chosen as hallmarks of the differences between the ERK1 and ERK2 morphant expression profiles. The expression levels were tested on the same RNA samples as used for the microarray analysis for cdh2 (cadherin 2, neuronal, NM_131081), mycn (v-myc, myelocytomatosis viral related oncogene, neuroblastoma derived, NM_212614), erm (ets related protein erm, NM_131205), cfos (FBJ murine osteosarcoma viral oncogene homolog, NM_205569), mos (moloney murine sarcoma viral oncogene homolog, NM_205580), snai1a (snail homolog 1a Drosophila, NM_131066) and vegf (vascular endothelial growth factor A, NM_131408) (Fig. 4). β-actin was taken as reference to determine the relative expression levels of the selected genes in ERK1MO, ERK2MO and standard control MO injected embryos. The obtained qPCR data show that the expression levels of cdh2, mycn, erm and snai1a are down-regulated in both ERK1MO and ERK2MO conditions, compared to standard control MO (Fig. 4A,B,C and 4F), whereas the expression level of mos is up-regulated in both ERK1MO and ERK2MO. The analysis of cfos and vegf confirmed the anti-correlated regulation comparing ERK1 and ERK2 knockdown to standard control MO conditions. Both genes are down-regulated by knockdown of ERK1 and up-regulated by knockdown of ERK2, compared to the expression-level of fos and vegf in standard control MO treated embryos (Fig. 4D,G).
In summary, the qPCR data confirmed the change in expression levels of the selected genes as observed by microarray analysis for all genes tested, thereby confirming the unique gene expression profiles for ERK1MO and ERK2MO mediated knockdown in early zebrafish development at 4.5 hpf (30% epiboly).
Gene Ontology (GO) analysis
The gene expression signatures of the ERK1 and ERK2 morphants were used to perform gene ontology (GO) analysis. This provides an unbiased biological gene enrichment analysis based on biological properties (GO-terms) assigned per gene. Gene ontologies describe gene products in terms of their associated biological processes (GO:0008150), cellular components (GO:0005575) and molecular functions (GO:0003674) in a species-independent manner. The results of this analysis showed a significant relative over- or under-representation of the number of Unigene IDs in ERK1 versus ERK2 morphants within the GO categories (Fig. 5). For ERK1 and ERK2 knockdown signature sets we obtained remarkable differences in the significantly enriched categories in the highest analyzed GO-level (level 4): for instance 5 vs. 14 enriched GO-terms are associated with Biological processes (Fig. 5A), 3 vs. 15 enriched GO terms are associated with cellular components (Fig. 5B) and 3 vs. 7 enriched GO terms are associated with Molecular functions (Fig. 5C), respectively.
Comparing the ERK1 and ERK2 knockdown signature sets various particular differences in over- or under-represented GO-terms were found. For example, both the GO-terms cell cycle (GO:0007049) and apoptosis (GO:0006915) are significantly enriched upon ERK2 knockdown. However, looking at the gene-lists in more detail inhibitory factors of apoptosis are mostly down-regulated, whereas positive regulators of cell cycle were up-regulated, indicating that apoptosis was not induced by the depletion of ERK2 at 30% epiboly (also see Fig. 2D) confirming our earlier conclusion from the time series. Cell adhesion (GO:0007155) and the cellular component GO terms 'tight junction' and 'cell junctions' are only significantly under-represented in the signature set of ERK2 morphants. Regulation of cell adhesion and the organization of tight- and cell-junctions are crucial for cell migration processes. Specifically for ERK1 knockdown a significantly enrichment of the 'translator regulator activity' (GO:0030528) GO-cluster was found. In contrast, the relative enrichment of this GO term in ERK2 morphants showed an under-representation. A significant overrepresentation of the GO term biosynthesis in ERK1 morphants correlates with these observations.
The GO-enrichment analysis showed that the number of genes within the GO-cluster 'development' (GO:0007275) were significantly under-represented for both ERK1 (19 genes) and ERK2 (136 genes) morphants. From the 19 development-related genes whose expression was affected by ERK1 knockdown, 12 genes (63%) were not found in the ERK2 knockdown signature set. This supports the notion that ERK1 and ERK2 may have distinct functions during embryogenesis by affecting the gene-expression of common and distinct genes sets during vertebrate development.
GenMAPP Pathways for zebrafish
To further analyze putative down stream targets of ERK1 and ERK2 involved in early development, we focused on essential signaling pathways that are involved in early embryonic differentiation and patterning; Nodal, FGF, Wnt and BMP-signaling pathways (Fig. 11) [12]. For our study, we used the signaling pathway analyzing software program, GenMAPP (Gene Microarray Pathway Profiler) [10, 11]. This program is designed for viewing and analyzing gene expression data in the context of biological pathways and allows microarray-mediated gene expression signature sets to be displayed on biological (signaling) pathways [10]. In contrast for human and mouse gene expression data-sets, where most signaling pathways are available for this program, there are no GenMAPP pathways based on zebrafish literature available yet for analyzing our gene expression datasets. Therefore, we first generated the in silico GenMAPP pathways for the zebrafish Nodal, FGF, (canonical) Wnt and BMP signaling pathways (Fig. 6, 7, 8, 9). This provides a valuable tool for the research community that makes use of zebrafish. The construction of these GenMAPP signaling pathways is based on what is specifically described in literature for zebrafish development, supported by the described knowledge for other vertebrate signaling processes and canonical signaling models, found on the Science's STKE Connections Map Database [14]. Although it is clear that the Nodal, FGF, Wnt and BMP pathways are all interconnected, resulting in a complex signaling network, we performed a pathway-based analysis focusing on separate signaling pathways since the ways these signaling pathways exactly interconnect on a molecular scale is hardly understood yet.
Pathway Analysis of ERK1MO and ERK2MO mediated knockdown expression profiles
The Unigene ID linked ERK1MO and ERK2MO signature sets that were used for GeneMAPP analysis were not limited by fold change but instead we used all genes that had a combined p-value for changed expression, compared to the standard control morpholino treated embryos, smaller than 10-5. As previously mentioned, the number of genes that showed a changed expression in ERK2MO compared to ERK1MO injected embryos was far larger. Therefore, as expected, more genes with changed expression levels were found in the in silico GenMAPPs signaling pathways for ERK2MO, than for ERK1MO.
Knockdown of ERK1 did show only one gene (smurf1) with a significantly changed expression level within our BMP signaling GenMAPP (Fig. 9). However, more genes were affected in FGF signaling: fgf17b (-1.37 fold) the MAPKKK mos, (+3.448 fold), transcription factor cmyc (-1.71 fold) and srf (serum response factor, -1.39 fold) showed significant changes in expression. In the nodal pathway, the Nodal antagonist lft1/antivin1 (+2.55 fold) and the EGF-CFC co-receptor oep (one eyed pinhead, -1.53 fold) were the only components found to be affected in ERK1 morphants. Furthermore, the ventrally expressed Wnt8-mediated organizer inhibitory gene vent [15] was down-regulated (-1.46 fold, Fig. 8). Other genes involved in Wnt-signaling affected by ERK1 knockdown were dab2 (disabled homolog 2, +1.47 fold), ck2b (casein kinase II beta subunit, -1.24 fold) and ppp2r5e1 (Protein phosphatase 2A, regulatory B subunit, B56, +1.30 fold). These genes are also considered to be involved in early embryonic pattering pathways. Two genes involved in regulating gastrulation cell migration, one-eyed pinhead (oep) and quattro [16, 17], were altered in expression.
The effect of depletion of ERK2 was far more severe in most of the analyzed signaling processes (Fig. 6,7,8,9). Key components of the FGF-pathway (fgf8, fgfr4 frs2, bRaf, aRaf and mek1l) and downstream target genes (erm, eve1, pea, mkp3, spry2, ntl, spt/tbx16 and tbx6) were down-regulated, indicating a block of the FGF-ERK pathway by ERK2 knockdown (Fig. 7). Expression of some of these (mesoderm) target genes is initiated by Nodal. The Nodal-genes like boz/dharma, squint/ndr1 and smad2 are up-regulated, whereas inhibiting genes lefty1 (lft1, -6 fold) and the ventral genes vox (-1.9 fold) and ved (-4 fold) are down-regulated in ERK2 morphants (Fig. 6). Other nodal signaling mediator genes that are down-regulated are oep (-4 fold), p300 (-2.03 fold), foxh1/sur (schmalspur, -2 fold) and the negative regulator of TGFβ signaling TGIF (-2 fold). The nodal-mediated endoderm gene sox32/casanova, expressed in the margin, was down-regulated (-6 fold), and also the downstream target-gene axial/foxA2 (-2 and -4 fold). Interestingly, squint/ndr1 also functions as a positive regulation of fibroblast growth factor receptor signaling pathway [18].
The Wnt ligand Wnt11 and receptors (frz7a, 7b, 8a, 9 and 10) and the central mediator β-catenin1 were down-regulated in ERK2 morphants, suggesting a severe inhibitory effect or even complete block of these pathways at this level (Fig. 8). This inhibition of the Wnt pathway is also supported by the up-regulation of axin2/conductin, a scaffold protein from the β-catenin destruction complex, responsible for the degradation of beta-catenin [19]. Down-regulation of the putative Wnt-target genes vox, vent, but also otx2, sp5, and lim1 further support impaired Wnt-signaling. However, ERK2 knockdown also led to the down-regulation of the inhibitors dkk1 and sfrp1, and up-regulation of the intracellular Wnt-signaling components fxd8c, dab2, β-catenin2 and tcf1.
The effect of ERK2 knockdown on BMP signaling is also complex, as bmp4 is up-regulated whereas bmp1a/tolloid and bmp6 are down-regulated (Fig. 9). This opposing effect is also found in the BMP antagonists, as chordin (chd) and the ventrally expressed membrane bound bmp-inhibitor bambi were down-regulated, whereas a different BMP antagonist gremlin is up-regulated. Adding to this complexity is the fact that the agonist twisted gastrulation (twsg1a) is up-regulated. The results clearly show that that dorsal-ventral patterning and also mesoderm patterning is severely affected but it is difficult to speculate about the downstream effects of all these changes of expression in the BMP pathway.
Biological confirmation of Pathway Analysis based prediction
To confirm predicted effects of the GenMAPP pathway analysis experimentally and to add information on the localization of expression, we performed whole mount in situ hybridization on ERK1 and ERK2 morphants at 30% epiboly with marker genes regulated by Nodal, BMP, Wnt and FGF (Fig. 6, 7, 8, 9 and Fig. 11). Different components of the Wnt-β-catenin pathway showed lower expression levels in ERK2 morphants. We showed that goosecoid (gsc) [20], a downstream marker gene for the Wnt pathway at early developmental stages (Fig. 10A–C) is not expressed in the ERK2 morphants. Knockdown of ERK1 did lead to a significant effect on the expression of gsc, but after knockdown of ERK2 no expression of gsc was detected by whole mount in situ hybridization. This confirms that canonical Wnt signaling was severely affected in ERK2 morphants, preventing subsequent expression of the Wnt-target gene gsc.
The lefty 1 (lft1, antivin1) gene is a member of the TGF-beta super-family that regulates left-right axis formation during embryogenesis via antagonistic activity against nodal, another TGF-beta super-family member. Expression starts at blastula stage, immediately after initiation of zygotic transcription, and is localized in the whole blastula margin at late blastula – 30% epiboly stage [21]. Whole mount in situ hybridization with lefty1 probe (Fig. 10D–F) at 30% epiboly shows a possible increase of lefty1 expression in ERK1 morphants (Fig. 10E), but the decrease of expression in ERK2 morphants (Fig. 10F) was clearly visible. As lefty1 is both an antagonists of Nodal signaling as well as a Nodal responsive gene, an increase of lft1 expression could mean that the signaling is sufficient and must be inhibited (ERK1MO), like in a wild type situation. A decrease in expression would mean that Nodal signaling not yet sufficient. Expression of mesoderm-genes in the margin indicates that (nodal mediated) mesoderm initiation took place in ERK2 morphants, however at a much reduced level (Fig. 10O,R).
In zebrafish, vox and vent interact with bozozok (boz), which is the earliest expressed dorsal-specific gene, and studies of boz embryos and the effects of ectopic boz expression indicate that it functions at the top of a hierarchy. Vox and vent are proposed to be repressors of boz expression since ectopic vox and vent eliminated the appearance of boz to establish the dorsal organize [15]. The expression signatures from the ERK1 and ERK2 morphants revealed that vox expression was not significantly changed in ERK1 morphants, but was down-regulated in ERK2 morphants., whereas for vent -expression this seemed to be opposite, as its expression was down-regulated in ERK1 morphants, but not significantly changed in ERK2 morphants (Fig. 6). The expression patterns of these genes revealed a possible reduction of vox expression in ERK1 morphants, which was more obvious on the putative dorsal side of the embryo where a clear cap was observed. The expression of vent was also reduced in ERK1 morphants and did not extend as far dorsally (K) compared to wild type embryos (J), indicating a mild dorsalization of ERK1 morphants. In ERK2 morphants, vox expression seemed to be reduced to a greater extend at the ventral side, but in the rest of the blastula the reduction of vox expression was not as significant and expression of vent was only detected at the ventral side of the blastula margin (L). In support of these finding, the expression of boz was found also to be up-regulated (+1.4 fold) in ERK2 morphants. Combined, these findings confirm that knockdown of ERK2 leads to impaired Wnt-mediated vox and vent expression which is reported to be involved in mesoderm patterning and maintenance.
The zebrafish ntl gene is, like the tbx6 gene, a member of the Brachyury-related T-box family of genes. Notail (ntl/brachyury) is involved in mesoderm development, as described in the legend to figure 11. At 30% epiboly ntl is expressed in the blastula margin [22]. This expression is synergistically regulated by FGF and Nodal signaling pathways [23, 24]. Both of these pathways show a negative regulation in the ERK2 morphants, as shown by the GenMAPP analysis (Fig. 6, 7, 8, 9). The negative effect on these pathways and the array-data itself suggested a down-regulation of the ntl-gene upon ERK2 knockdown, and was confirmed by whole mount in situ experiments (Fig. 10M–O). The ntl gene expression in ERK1 morphants was comparable to expression in wild type embryos, but ntl expression was decreased in ERK2 morphants. Strikingly, expression of ntl was not constant in the marginal ring, as stronger expression was detected in the putative dorsal side of the ERK2 morphants.
Tbx6 is exclusively expressed in the ventral mesendoderm and its expression is linked to ventral mesoderm specification [25]. In ERK1 morphants the in situ hybridization experiment showed that tbx6 expression was not extended as far dorsally as in wild type embryos, as tbx6 expression at the putative dorsal side of these embryos was severely reduced (Fig. 10P,Q). In ERK2 morphants, tbx6 expression was greatly reduced and was only detected at the ventral margin (Fig. 10R). In ERK2 morphants tbx6-expression an even more severe reduction of tbx6 expression was down-regulated (-2.3 fold).
The obtained results by whole mount in situ hybridization using gsc, lft, vox, vent, ntl and tbx6, confirm or support the predictions made by the GenMAPP analysis, as the changes in their expression levels are in agreement with the predictions obtained by the signaling pathway analysis of the microarray data.