Early embryogenesis is the sensitive time window for vangl2 mutants
To determine the critical time window of ethanol sensitivity for vangl2 mutants and heterozygotes, we first initiated ethanol treatment at various stages comprising late blastula to early gastrula for 24 h. The inner lens-to-lens width was used as a measure of synophthalmia, and we calculated the occurrence of cyclopia to characterize ethanol-induced teratogenesis. In control conditions, the spacing of the eyes in vangl2 mutants resembled that of wild-type zebrafish. However, we note that the facial phenotypes in vangl2 mutants are variable and midline defects can occur even under control conditions [15]. Ethanol-exposed vangl2 mutants exhibited midline defects ranging in severity from synophthalmia to cyclopia across all time points examined, but these mutants were fully penetrant for cyclopia (100% fused; n=5/5) when ethanol was applied at shield stage (6 h post-fertilization, hpf) at the onset of gastrulation (Fig. 1a-b). Interestingly, heterozygotes only displayed cyclopia when ethanol was applied at a high stage (3.3 hpf) (22% fused; n=4/18), a time when treating wild-type embryos with higher concentrations of ethanol causes similar defects (Fig. 1b) [16, 17]. Thus, heterozygotes and homozygotes are divergent in their time window of greatest ethanol sensitivity. This may be due to a compensatory genetic mechanism in vangl2 heterozygotes, as zygotic gene expression initiates after high stage (4 hpf) in zebrafish.
The effect of ethanol on the early zebrafish transcriptome is subtle relative to the developmental age
To determine if ethanol caused transcriptional changes that could underlie the interaction with vangl2, we designed two RNA-seq experiments that largely overlapped in design (Fig. 2a-b). One possibility for the variability in the phenotypic outcome of ethanol exposure is variable transcriptional response across individuals. Therefore, we performed RNA-seq on single wild-type AB strain embryos to determine if ethanol exposure variably destabilized the transcriptome. Embryos were exposed to 1% ethanol in embryo media (171 mM), which equilibrates to approximately 50 mM tissue concentration [18]. For our RNA-seq experiments, ethanol exposure initiated at 6 hpf, when mutants were the most sensitive. Each sample consisted of an individual zebrafish embryo with five biological replicates per timepoint and treatment. For the first experiment, embryos were collected at mid-gastrulation (8 hpf) and the end of gastrulation (10 hpf) (Fig. 2a). A second experiment was performed to increase power and to include a time when the eye fields have completely separated, 14 hpf (Fig. 2b). The 6 hpf control samples were omitted in the second experiment since they lacked ethanol-treated samples for comparison. The data from both experiments were combined for subsequent analyses, controlling for batch effects.
To assess the effect of age, batch (experiments 1 and 2), and ethanol, on the early zebrafish transcriptome, we performed principal component analysis (PCA). Individuals of similar age clustered tightly along PC1, which accounted for 39% of the transcriptional variation observed across the datasets (Fig. 2c). This strong effect of age on the zebrafish transcriptome is in agreement with previous studies [19]. Clustering of samples by age, regardless of ethanol treatment suggested that control and ethanol-treated samples were accurately staged, indicating that ethanol did not delay developmentally regulated transcriptome patterns. There was greater discrimination of 14 hpf samples relative to earlier timepoints (6, 8, and 10 hpf), which is consistent with the greater distinction of this timepoint in terms of developmental time and morphology (Fig. 2c). PC2 largely captured batch effects between experiments 1 and 2 (Fig. 2d).
The majority of variation between samples did not appear to be due to treatment, with control and ethanol-treated samples randomly interspersed along PC1 and PC2 (Fig. 2e). Separation by treatment was observed along PC8 and PC9, which accounts for 3% of the variation in the data (Fig. 2f). Hierarchical clustering of samples based on correlation further corroborated this finding. The 6 and 14 hpf samples showed the greatest dissimilarity, whereas the 8 and 10 hpf samples showed the greatest similarity, irrespective of treatment (Additional file 1: Fig. S1). In summary, the transcriptional effect of exposure to this dose of ethanol on the early zebrafish transcriptome was subtle, while age provided the strongest transcriptional fingerprint.
Ethanol has effects on transcription that are largely distinct between different developmental timepoints
Although developmental age was a stronger source of transcriptional variation, we still detected substantial variation in gene expression following ethanol exposure. There were 1414 differentially expressed genes (DEGs), with a false-discovery rate (FDR) less than 0.1 (Fig. 3a; Benjamini–Hochberg procedure). There were more upregulated than downregulated DEGs among ethanol-treated individuals across timepoints (Additional file 2: Fig. S2A-D), and these DEGs were distinct between developmental timepoints (Fig. 3b). In summary, while some genes are generally affected by ethanol across several critical periods, ethanol largely elicits unique responses at distinct developmental stages spanning early embryogenesis.
Ethanol does not affect the Wnt/PCP pathway at the transcriptional level
While vangl2 mutants do not typically exhibit cyclopia, cyclopia is observed in double mutants between vangl2 and other Wnt/PCP pathway members [15, 20]. One potential mechanism for the interaction between ethanol and vangl2 would be the transcriptional misregulation of other Wnt/PCP pathway members. However, transcription of Wnt/PCP pathway members was largely unaffected by ethanol exposure. No Wnt/PCP pathway members were among the 14 shared DEGs across all timepoints (Fig. 3c). KEGG pathway enrichment analysis further confirmed that ethanol had little effect on the Wnt/PCP pathway at the level of transcription (Fig. 3d). Only two Wnt/PCP pathway members were moderately affected: ethanol exposure moderately decreased expression of the cofactor, glypican 4 (gpc4), (log2 fold= −0.237; p value = 0.036) across all timepoints [15, 21] and increased expression of rac3a (log2 fold = 0.737; p value = 8.89E−06), a member of the Rho family of small GTPases [22]. We plotted the normalized read counts across each time point from the RNA-seq experiment to further investigate the dynamics of the potential alteration in gpc4 levels. The expression of gpc4 was modestly reduced in ethanol-treated embryos across the RNA-seq dataset (Fig. 3e). Ethanol exposure consistently downregulated gpc4 at 8 and 10 hpf, but the magnitude of the downregulation was relatively modest. To statistically compare ethanol-treated and control embryos, expression of gpc4 at 10 hpf was investigated using real-time quantitative reverse transcription PCR (RT-qPCR). This result demonstrated that gpc4 was not significantly affected by ethanol exposure (p = 0.4631) (Fig. 3f). Thus, direct transcriptional alteration to the Wnt/PCP pathway is unlikely to explain the ethanol-induced phenotypes in vangl2 mutants and heterozygotes.
Modules of co-regulated genes related to ethanol exposure
To determine if ethanol disrupted networks of genes could explain the vangl2-ethanol interaction, we next performed Weighted Gene Co-expression Network Analysis (WGCNA) [23]. This unsupervised network analysis identifies groups of genes, termed “modules,” based on correlated expression patterns across the samples. Modules are summarized by the first principal component for the expression estimates of the included genes, termed the “module eigengene,” which can be correlated with sample traits to identify biological significance. The cluster dendrogram generated in this analysis illustrates the presence of highly distinct and clustered modules (Additional file 3: Fig. S3A). The merging of similar modules produced eleven total modules (Additional file 3: Fig. S3B). Consistent with our PCA analysis, the module eigengenes were primarily correlated with time. We used these modules to further validate our ethanol exposure regimen.
Previous work has shown that ethanol delays development in a dose-dependent manner at concentrations equal to or greater than 1.5% ethanol [24]. For the experiments herein, ethanol-treated samples were morphologically stage-matched to control samples to exclude differences due to developmental age or delay. To confirm that the ethanol samples were indeed age-matched to the control samples at the level of transcription, we compared expression patterns of developmentally regulated genes between the ethanol and control samples from each timepoint. For developmentally regulated genes, we used the gene with the highest module membership, the “hub gene,” from each of the WGCNA modules that was associated with time (p < 0.05) (Additional file 4: Fig. S4). Consistent with the results from the PCA, the slope of the expression levels for each gene were similar in control and ethanol-treated samples across age, most clearly demonstrated in the magenta4, thistle1, and blueviolet modules (Additional file 4: Fig. S4B-C, and G). Interestingly, we find a time-specific difference in the expression of hub genes in the honeydew1, greenyellow, and saddlebrown modules at 10 hpf (Additional file 4: Fig. S4F, H, and K). Together, these data indicate samples were accurately age-matched and the observed changes were biologically relevant.
Two modules (mediumpurple4 and darkolivegreen4) correlated with ethanol treatment (Additional file 3: Fig. S3B-C). The purple and green modules positively and negatively correlated with ethanol treatment, respectively (Additional file 3: Fig. S3D). A significant differentially expressed gene from the green (Additional file 3: Fig. S3E) and purple module (Additional file 3: Fig. S3F) were selected for independent validation on biological replicate samples derived from the same wild-type zebrafish line using RT-qPCR. These results indicate that our RNA-seq faithfully represents transcript levels. The green module is only weakly downregulated and is enriched for genes encoding zinc finger (ZnF) proteins, which included znf1015 (Additional file 3: Fig. S3G). This module is almost entirely composed of genes on chromosome 4 that have been shown to be co-regulated during these stages of development [25, 26]. The purple module revealed GO enrichment of transmembrane transporters, which included slc16a9a (Additional file 3: Fig. S3H). As an upregulated module enriched in transporters, these genes are likely to represent the physiological response to ethanol. Thus, while we did identify two modules associated with ethanol (Additional file 5: Table S1), the module membership did not provide clues to the interaction between ethanol and vangl2.
Cyclopamine is predicted to mimic the effects of ethanol
One challenge in RNA-seq analyses is inferring the mechanism underlying a diseased or environmentally perturbed state from a large set of differentially expressed genes. Individual functional analyses of significant gene-ethanol interactions are time-consuming and thus inefficient. To circumvent this problem, we adopted a bioinformatics approach, utilizing the Library of Integrated Network-Based Cellular Signatures (LINCS L1000) toolkit [27]. This transcriptomic dataset includes gene expression data from nine human cancer cell lines exposed to thousands of small molecule drugs [27]. We queried the top 100 and 150 upregulated and downregulated genes induced by ethanol exposure against the LINCS L1000 dataset using the clue.io platform (https://clue.io). The query generated a list of small molecules predicted to have a positive or negative correlation to the input signature (i.e., small molecules predicted to either mimic or antagonize the transcriptional effects of ethanol). We were particularly interested in those chemicals with a positive correlation to ethanol as they would give insight into the mechanism of ethanol teratogenicity and highlight potential co-factors that exacerbate ethanol teratogenicity.
Interestingly, cyclopamine, a hallmark Shh pathway inhibitor that inhibits the core Shh pathway protein Smoothened (Smo), was predicted to positively correlate with the ethanol signature. The Shh pathway is critical for midfacial development and mice deficient in Sonic Hedgehog (Shh) exhibit severe brain and face malformations, including holoprosencephaly and a single-medial eye (cyclopia) [28]. In zebrafish, reduction of shh or null mutations in smo similarly results in severe loss of craniofacial midline structures of the anterior neurocranium [29]. Animal models have demonstrated ethanol is an environmental risk factor for holoprosencephaly, resulting in a characteristic set of midfacial defects, a hypomorphic forebrain, and in severe cases, cyclopia [16, 17, 30,31,32,33]. Thus, attenuation of the Shh pathway could mechanistically explain the vangl2-ethanol interaction.
Ethanol indirectly attenuates Shh signaling
To investigate the interaction of cyclopamine and ethanol, we first exposed wild-type zebrafish embryos to cyclopamine (50 μM) at shield stage (6 hpf) for 24 h, mimicking the ethanol exposure window for the vangl2 mutants. Embryos were fixed at 4 dpf and the cartilage and bone were stained with Alcian blue and Alizarin red, respectively. We observed a range of midfacial defects with the most severe phenotype being a complete loss of the anterior neurocranium and reduced spacing between the eyes (Fig. 4a). Since cyclopamine was predicted to mimic the effects of ethanol, we next combined this low dose of cyclopamine with ethanol. Strikingly, all embryos presented with synophthalmia or cyclopia and significant reductions and defects of the cartilages of the neuro- and viscerocranium. To quantify this combinatorial effect, we measured the inner lens-to-lens width and observed a significant reduction in co-exposed embryos relative to those exposed to either ethanol or cyclopamine alone (p < 0.0001), suggesting a strong synergistic interaction (Fig. 4b).
In chick, ethanol exposure during somitogenesis has been proposed to suppress Shh signaling and induce apoptosis in cranial neural crest cells that make up the craniofacial skeleton [34]. However, work in zebrafish shows only a modest increase in cell death within the eye field at a much higher dose of 2% ethanol [35]. To ensure that cells in the eye field are not simply undergoing apoptosis, we performed a TUNEL cell death assay in ethanol-treated vangl2 mutants at 11 hpf, prior to optic vesicle evagination. As expected, we failed to detect an increase in apoptotic cells within the eye field in ethanol-treated vangl2 mutants or their siblings (Additional file 6: Fig. S5). Thus, Shh-mediated pro-survival signals do not appear to be significantly disrupted by ethanol.
To test for a direct effect of ethanol on Shh signaling, we quantified the relative expression of patched 2 (ptch2), a canonical read-out of Shh pathway activity, at bud stage (10 hpf). Ethanol exposure had no effect on the expression levels of ptch2 (p = 0.9966), consistent with RNA-seq results (Fig. 4c). While cyclopamine significantly reduced Shh signaling (p = 0.0006), ethanol did not further reduce ptch2 levels significantly (p = 0.1115). KEGG analysis from the RNA-seq confirmed that ethanol does not affect the Shh pathway as a whole, nor does it reduce the levels of direct targets of Shh signaling (i.e., gli, ptch) (Fig. 4d). Thus, at this dose, ethanol does not appear to directly attenuate Shh signaling.
Ethanol disrupts convergent extension
The ethanol-induced vangl2 mutant phenotype closely mirrors those in compound mutants between vangl2 and other Wnt/PCP pathway members [15, 20]. These double mutants display a further reduction in convergent extension, as evidenced by a shorter and broader body axis [15]. Based on these data, we hypothesized that ethanol disrupts convergent extension, which would mislocalize the Shh signal and result in eye defects.
We performed in situ hybridization on control and ethanol-treated vangl2 embryos to examine convergent extension and to test the hypothesis that the Shh expression domain is mislocalized in ethanol-exposed vangl2 mutants. We analyzed the expression of shha in the axial mesoderm and paired box 2a (pax2a) in the midbrain-hindbrain boundary at bud stage (10 hpf) (Fig. 5a) [36]. Ethanol was initiated at a high stage (3.3 hpf). We chose this timepoint because it is when heterozygotes are most sensitive to ethanol-induced synophthalmia. Because vangl2 mutants have convergent extension defects and display synophthalmia, we reasoned that effects may be more apparent in heterozygotes. We observed a gene and ethanol-dose-dependent reduction in the length of the shha expression domain (extension) and an increase in the width of the pax2a expression domain (convergence).
To quantify the effect of ethanol on convergent extension, we plotted the normalized expression values of shha/pax2a (Fig. 5b). Post hoc analyses (Tukey’s) revealed significant differences between control embryos across vangl2 genotypes, indicating that vangl2 gene dosage affects convergent extension. Similarly, we found significant differences between ethanol-treated mutants and their heterozygous (p = 0.0037) and wild-type siblings (p = 0.0052). This observation confirms that loss of vangl2 results in reduced convergent extension movements, as evidenced by their short body axis. This phenotype was exacerbated with ethanol treatment, resulting in reduced convergent extension for vangl2 heterozygotes and homozygotes, compared to their untreated counterparts. While there was a trend, we did not observe a difference in convergent extension between ethanol-treated vangl2 heterozygotes and their wild-type siblings in their shha/pax2a ratio (p = 0.9554). However, synophthalmia occurred in ethanol-treated heterozygotes but not wild-type embryos, demonstrating that the nonsignificant reduction in convergent extension in the heterozygotes can have a phenotypic consequence.
Separation of the eye field into two eye primordia is accomplished via Shh signaling to the anterior neural plate. We measured the distance between shha and dlx3, a marker of the anterior neural plate (Fig. 5a and c), to determine if ethanol alters the source and target of Shh. Here again, we find a combinatorial effect of ethanol and vangl2 gene dosage, with ethanol-exposed mutants having the greatest distance between the shha and dlx3 expression domains (Fig. 5c). Ethanol exposure has a significant effect on this distance in wild-type embryos (p = 0.003). The distance is slightly, albeit non-significantly, elevated in ethanol-treated mutants, compared to unexposed mutants (p = 0.7133). Given that vangl2 mutants are predisposed to cyclopia, this slight increase may again, have phenotypic relevance. These results are consistent with a model in which the interaction between ethanol and vangl2 is due to an indirect effect of ethanol on the Shh pathway via combined disruption of convergent extension cell movements.
Ethanol alters six3 and rx3 expression in the eye field
Transcription factors involved in the specification of the eye field have also been implicated in the mechanism of eye field separation. The expression of three of these transcription factors, six3, rx3, and rx1 is altered by ethanol exposure [35]. We examined the effect of ethanol on six3a and rx3 in ethanol-treated vangl2 mutants at the initiation of optic vesicle evagination. In situ hybridization of six3a in 11 hpf embryos shows a heart-shaped expression pattern in the prospective forebrain. The caudal indentation marks the splitting of the eye field into bilateral domains (Fig. 5d). This expression pattern becomes more diffuse with loss of vangl2. In untreated homozygous mutants, we observed a shortening along the anterior-posterior (AP) axis and a broadening along the mediolateral axis at 11 hpf, consistent with the convergent extension defect. This expression pattern was further exacerbated in ethanol-treated mutants with complete loss of the caudal indentation. We observed a similar effect of genotype and ethanol on rx3 expression at mid-evagination (12 hpf) (Additional file 7: Fig. S6). At this stage, rx3 is localized to the prospective forebrain and retina [37]. Ethanol-treated vangl2 homozygous mutants exhibit a compressed expression domain, clearly displaying a reduction in convergent extension. Our expression analyses suggest the eye field is specified but mislocalized and fails to separate into bilateral domains, likely due to defects in mesodermal migration. However, we cannot rule out a specific transcriptional effect of ethanol on a small cell population that receives Shh signaling, and scRNA-seq would prove useful to identify such interactions.
Mutation in gpc4 enhances cyclopia in a dose-dependent manner
Our data support a model in which ethanol interacts with vangl2 via a combinatorial disruption of convergent extension. If ethanol disrupts convergent extension, which leads to interactions with vangl2, then further genetic disruption to convergent extension should exacerbate the ethanol-vangl2 phenotype. Embryos deficient in gpc4 similarly have a defect in convergent extension as evidenced by their shortened, broadened body axis [15]. Previous work in zebrafish has shown a functional interaction between vangl2 and gpc4, where vangl2;gpc4 double mutants were invariably cyclopic [15]. We conducted additional functional analyses to further examine the relationship between these two genes in the context of ethanol exposure. Consistent with Marlow et al., double mutants were fully penetrant for cyclopia with or without ethanol. While 1% ethanol exposure altered the facial morphology of gpc4 mutants, it did not cause cyclopia (Additional file 8: Fig. S7A). However, ethanol concentrations greater than 1%, which did not cause cyclopia in wild-type siblings, resulted in cyclopia in gpc4 mutants (Additional file 8: Fig. S7B) [15]. Additionally, embryos carrying 3 mutant alleles (either vangl2m209/m209;gpc4fr6/+ or vangl2m209/+;gpc4fr6/fr6) were greatly sensitized to cyclopia when exposed to 1% ethanol (Additional file 8: Fig. S7C). Collectively, these data indicate that ethanol disrupts convergent extension, and in sensitized genotypes, mispositions a source of Shh that is essential for the separation of the eye field into bilateral primordia.
Blebbistatin phenocopies the ethanol-induced defect in vangl2 mutants
If ethanol is having an effect on convergent extension, a known inhibitor of convergent extension should similarly interact with vangl2. We extended our analyses to test the effect of blebbistatin, a myosin II inhibitor that disrupts cell movements and reduces axis elongation, in vangl2 mutants [38]. Blebbistatin treatment during gastrulation (6–10 hpf) induced synophthalmia in vangl2 heterozygotes and homozygotes (Fig. 6a). We then quantified this interaction by measuring the inner lens-to-lens width and found a significant reduction in blebbistatin-treated vangl2 homozygotes (Fig. 6b). Furthermore, 26.92% and 100% of blebbistatin-treated heterozygotes and homozygotes, respectively, were cyclopic (Fig. 6c). These data suggest ethanol further inhibits convergent extension in the sensitized background, vangl2, which not only disrupts their ability to elongate their body axis, but also inhibits the splitting of the eye field.
Ethanol affects vangl2 filopodia organization on migrating cells
Based on our data, cell dynamics underlying convergent extension are likely involved in the vangl2-ethanol interaction. During gastrulation, mesodermal cells adopt a bipolar shape and orient their actin-based cytoskeletal projections on their leading and trailing edge, with respect to the notochord, allowing for proper convergent extension [39]. In vangl2 mutant embryos, these gastrula cells fail to elongate and align with the primary axis [40]. Furthermore, vangl2 mutant ectodermal cells possess a greater number of filopodial protrusions than their wild-type siblings and display projections in an unpolarized manner [41].
To determine whether ethanol may be exacerbating this phenotype, we analyzed filopodia protrusion length and orientation around the circumference of migrating live mesodermal cells. We injected memGFP mRNA into a single blastomere of vangl2 embryos at the 8–32 cell stage, to obtain mosaic expression. We then treated a subset of the injected embryos with 1% ethanol from 6 to 10 hpf. Homozygous mutants were identified based on their reduced body axis elongation at 10 hpf. Cells were imaged at 40–60° from the primary axis. Spike-like filopodia were highly unpolarized in ethanol-treated vangl2 mutants compared to their wild-type and heterozygous siblings (Fig. 7a). We quantified the polarity and number of filopodia for each cell, with respect to the head (90o), the primary axis/notochord (180o), and the tail (270°) (Fig. 7b). Rose diagrams suggest ethanol-treated vangl2 mutants have more un-polarized projections on the anterior/posterior edge (red) as opposed to the leading (green) or trailing (blue) edge. We next reduced the number of bins in the rose diagrams to 90° and noticed a clear over-representation of filopodia in vangl2 mutants on their anterior/posterior edge (Fig. 7c). Tukey’s post hoc analyses confirmed ethanol-treated vangl2 mutants had significantly more filopodia on their anterior/posterior edge compared to other axes and their wild-type or heterozygous siblings (Fig. 7d). These data confirm that ethanol exposure is associated with an increase in the number of filopodial projections oriented in the wrong direction, or perpendicular to the path of dorsal migration, in ethanol-treated vangl2 mutants.