Jasmijn A Baaijens

Leen Stougie

Alexander Schönhuth
email: alexander.schoenhuth@cwi.nl

Strain-aware assembly of genomes from mixed samples using variation graphs

Keywords: viral quasispecies, haplotype reconstruction, abundance estimation, variation graph

The goal of haplotype-aware genome assembly is to reconstruct all individual haplotypes from a mixed sample and to provide corresponding abundance estimates. We provide a referencegenome-independent solution based on the construction of a variation graph, capturing all diversity present in the sample. We solve the contig abundance estimation problem and propose a greedy algorithm to efficiently build maximal-length haplotypes. Finally, we obtain accurate frequency estimates for the reconstructed haplotypes through linear programming techniques. Our method outperforms state-of-the-art approaches on viral quasispecies benchmarks and has the potential to assemble bacterial genomes in a strain-aware manner as well.

Background

Many sequencing data sets contain fragments of several closely related sequences, such as bacterial mixtures, environmental samples (metagenomics), or mixtures of viral strains, as extracted from infected hosts (viral quasispecies [START_REF] Domingo | Viral quasispecies evolution[END_REF]). High mutation rates affecting the species contained in such a sample usually yield that a particular species comes as an ensemble of closely related genomes each of which pertains to a particular strain of the species. Because relevant properties such as escape from medical treatment or host immune response happen at the strain, and not the species level [START_REF] Crotty | RNA virus error catastrophe: direct molecular test by using ribavirin[END_REF][START_REF] Vignuzzi | Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population[END_REF], one of the current most driving analytical challenges is to reconstruct the individual genomes present in a such a sample at the strain level. When emphasizing the genetic variation that characterizes them, these individual genomes are usually referred to as haplotypes. Therefore, this challenge is more generally also referred to as haplotype-aware genome assembly.

It is further characteristic that the relative abundance (or relative frequency) of the strains contained in such mixed samples varies across the strains. It is possible that some strains dominate the sample by their high relative frequency and, as a consequence, mask others that come at low relative frequency. Also because such low-frequency bacterial or viral strains are often relevant factors when considering escape from medical treatment or host immune response [START_REF] Crotty | RNA virus error catastrophe: direct molecular test by using ribavirin[END_REF][START_REF] Vignuzzi | Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population[END_REF], obtaining reliable, accurate estimates of the relative abundances of strains can be crucial when analyzing mixed samples.

In this paper, we present VG-flow, a method for reconstructing genomes at the strain level with integrated abundance estimation. Because genomes from mixed samples are usually affected by substantial mutation rates, which can decisively hamper the use of linear reference genomes because of the biases they induce when dealing with genetic variation, VG-flow avoids using standard linear reference genomes altogether. Instead, VG-flow is based on variation graphs as underlying, flexibleto-construct reference systems that account for haplotype-specific mutations in an unbiased manner. We construct these variation graphs in a de novo manner, that is, without resorting to any external means such as a reference genome.

In this, first and foremost, VG-flow presents a comprehensive solution to the de novo viral quasispecies assembly problem. Thereby, VG-flow comes with significant improvements over earlier related work presented in the literature. Note that in comparison to bacteria and eukaryotes, viruses have relatively short genomes, subject to mutation rates that exceed those for bacterial genomes; this is even more expressed for RNA viruses than for DNA viruses [START_REF] Duffy | Why are RNA virus mutation rates so damn high?[END_REF]. In particular the shortness of the genomes simplifies the analysis, and allows for solutions that are specific to reconstructing viral quasispecies.

Thanks to the computational efficiency of the underlying method, establishing a novelty from a theoretical point of view, VG-flow also establishes the first solution that can also be applied to bacterial sized genomes. This points out its potential for applications in metagenomics type settings, beyond its applicability for strain-aware assembly of virus genomes. Note that reconstructing all of the individual haplotypes present in a metagenome, that is de novo, strain-aware metagenome assembly, is challenging also because of various other reasons, all of which require specialized tools and further methodological progress [START_REF] Sczyrba | Critical assessment of metagenome interpretation -a benchmark of metagenomics software[END_REF]. In that respect, VG-flow offers an important, and so far missing building block for a full solution of this problem.

VG-flow takes as input a next-generation sequencing (NGS) data set and a collection of strainspecific contigs assembled from the data, and produces full-length haplotypes and corresponding abundance estimates. Our method is centered on estimating contig abundances in a contig-variation graph, a graph that captures all quasispecies diversity present in the contigs [START_REF] Paten | Genome graphs and the evolution of genome inference[END_REF][START_REF] Garrison | Variation graph toolkit improves read mapping by representing genetic variation in the reference[END_REF]. We build a flow network to accompany the variation graph and estimate contig abundances by solving a flow-like optimization problem: variables represent flow values on the edges of the flow network and we impose flow constraints, while the objective function evaluates the difference between estimated contig abundances and read coverage for every node in the variation graph. This objective function is convex, which renders the flow problem polynomial time solvable [START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF].

The flow solution presents abundance estimates for the input contigs, which are of value in its own right in various mixed sample applications [START_REF] Lindner | Metagenomic abundance estimation and diagnostic testing on species level[END_REF][START_REF] Fischer | Abundance estimation and differential testing on strain level in metagenomics data[END_REF][START_REF] Bray | Near-optimal probabilistic RNA-seq quantification[END_REF]. We use the contig abundance estimates in a combination of greedy algorithms to extract candidate haplotypes from the variation graph, where candidate haplotypes reflect concatenations of subpaths associated with the input contigs. Finally, we solve an optimization problem whose variables represent the haplotype abundances and the difference between read coverage and haplotype abundance is minimized over all nodes. Thus, we obtain a selection of candidate haplotypes that represents the quasispecies, along with haplotype abundance estimates.

Existing viral quasispecies assemblers include widely evaluated tools like [START_REF] Prabhakaran | HIV haplotype inference using a propagating dirichlet process mixture model[END_REF][START_REF] Zagordi | ShoRAH: estimating the genetic diversity of a mixed sample from next-generation sequencing data[END_REF][START_REF] Prosperi | QuRe: software for viral quasispecies reconstruction from next-generation sequencing data[END_REF], as well a variety of methods introduced more recently [START_REF] Ahn | aBayesQR: A bayesian method for reconstruction of viral populations characterized by low diversity[END_REF][START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF][START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF][START_REF] Barik | Qsdpr: Viral quasispecies reconstruction via correlation clustering[END_REF][START_REF] Chen | De novo haplotype reconstruction in viral quasispecies using paired-end read guided path finding[END_REF][START_REF] Knyazev | CliqueSNV: Scalable reconstruction of intra-host viral populations from NGS reads[END_REF]. These methods can be divided into two classes: reference-guided and reference-free (also referred to as de novo). De novo approaches do not require any prior information, such as a reference genome or knowledge of the quasispecies composition. This has been shown to have advantages over reference-guided reconstruction, since using a reference genome can induce significant biases. Especially at the time of a viral disease outbreak, an appropriate reference genome may not be available due to high mutation rates. However, most of the aforelisted tools are reference-guided; only [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF], [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF], and [START_REF] Chen | De novo haplotype reconstruction in viral quasispecies using paired-end read guided path finding[END_REF] present de novo approaches.

Moreover, many of these specialized viral quasispecies assemblers aim at single gene reconstruction, rather than whole genome assembly. In [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF] we took a first step towards full-length de 2 novo viral quasispecies assembly. There, we have shown that de novo haplotype reconstruction with integrated haplotype abundance estimation yields assemblies that are more complete, more accurate, and provide better abundance estimates. While this approach is guaranteed to find a selection of haplotypes that is optimal in terms of being compatible with the read coverages, its runtime is exponential in the number of contigs. Since the number of contigs generally increases on increasing genome length, we found [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF] unsuitable for genomes larger than ∼ 10 kb. With VG-flow, we provide a reference-free solution to the full-length viral quasispecies reconstruction problem that scales well to longer genomes. As another benefit of the theoretical rigorosity of our problem formulation and efficiency of the solution, we also experience considerable improvements in terms of accuracy compared to existing tools. Some of the challenges that have to be dealt with in viral quasispecies assembly can also be found in RNA transcipt assembly, where the goal is to reconstruct an unknown number of transcripts and predict the relative transcript abundances. Not surprisingly, many RNA transcipt assemblers define graph optimization problems similar to our flow formulation [START_REF] Tomescu | A novel min-cost flow method for estimating transcript expression with RNA-seq[END_REF][START_REF] Rizzi | On the complexity of minimum path cover with subpath constraints for multi-assembly[END_REF][START_REF] Bernard | Efficient RNA isoform identification and quantification from RNA-seq data with network flows[END_REF][START_REF] Pertea | StringTie enables improved reconstruction of a transcriptome from RNA-seq reads[END_REF][START_REF] Trapnell | Transcript assembly and quantification by RNA-seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF]. Although dealing with related problems, these methods cannot be applied in a viral quasispecies setting so easily: they require a collection of reference genomes representing all possible haplotypes as input, which is not available in our setting. Nevertheless, the theory behind these approaches is related to what we do. In [START_REF] Tomescu | A novel min-cost flow method for estimating transcript expression with RNA-seq[END_REF], node and edge abundance errors are used to define a min-cost flow problem; note that this formulation does not take subpath constraints into account. On the other hand, [START_REF] Rizzi | On the complexity of minimum path cover with subpath constraints for multi-assembly[END_REF] describes how subpath constraints can be incorporated into a minimum path cover formulation. This results in an optimization problem that is solvable in polynomial time, but does not minimize node abundance errors. We use the best of both worlds by defining an optimization problem that takes subpath constraints, minimizes node abundance errors, and is polynomially solvable; this establishes a theoretical novelty. Because this novelty gives way to different types of analyses in other settings, and immediately connects to extensively treated theoretical issues [START_REF] Rizzi | On the complexity of minimum path cover with subpath constraints for multi-assembly[END_REF][START_REF] Tomescu | A novel min-cost flow method for estimating transcript expression with RNA-seq[END_REF], we feel that it is of value also in its own right.

Among more generic assemblers, [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF] has been shown to be capable of reconstructing individual haplotypes from mixed samples, up to a certain degree. This method was designed for bacterial genomes and scales well to human genomes, but is unable to reconstruct low-frequent haplotypes [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF]. Haplotype-aware assembly of metagenomes is a big challenge, which tends to result in scattered genome fragments and missing strains [START_REF] Sczyrba | Critical assessment of metagenome interpretation -a benchmark of metagenomics software[END_REF]. Metagenomic assemblers such as [START_REF] Nurk | metaSPAdes: a new versatile metagenomic assembler[END_REF][START_REF] Boisvert | Ray meta: scalable de novo metagenome assembly and profiling[END_REF][START_REF] Peng | Meta-IDBA: a de novo assembler for metagenomic data[END_REF][START_REF] Li | MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph[END_REF] aim to reconstruct mixtures of viral and bacterial populations at strain level. The contigs obtained with these methods, or any other assembler, can also be used as input for VG-flow. Although we emphasize the reconstruction of viral quasispecies, the mathematical framework presented here is generic and could be applied in other scenarios as well; as such, VG-flow has the potential to make a big step ahead in haplotype-aware genome assembly in general.

Results

We present VG-flow, a new approach to haplotype aware genome assembly from mixed samples. This algorithm takes as input a data set of next-generation sequencing reads and a collection of strain-specific contigs; note that de novo assembly into strain-specific contigs can be performed using various tools (e.g. [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF][START_REF] Chen | De novo haplotype reconstruction in viral quasispecies using paired-end read guided path finding[END_REF][START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF][START_REF] Nurk | metaSPAdes: a new versatile metagenomic assembler[END_REF], depending on the application). The output of VG-flow consists of maximal length haplotypes along with relative abundance estimates for each of these sequences. In addition, the algorithm yields abundance estimates for each of the input contigs.

Our approach consists of five steps, as depicted in Figure 1:

(1) We construct a contig variation graph V G C by performing Multiple Sequence Alignment (MSA) on the input sequences. Node abundances are obtained by mapping sequencing reads (2) We build a flow network F G using V G C .

(3) We define and solve a flow-like optimization problem on F G to obtain contig abundance estimates.

(4) We generate a set of candidate haplotypes P cand based on the estimated contig abundances through multiple greedy heuristics.

(5) We obtain a selection of haplotypes H from P cand by solving another linear optimization problem, defined on V G C . The solution to this problem presents estimates for the relative abundances of all candidate haplotypes in P cand , thereby eliminating any false haplotypes.

The final output is presented as a genome variation graph V G H capturing the haplotypes in H, along with the estimated relative abundances.

Steps (1) and ( 5) are based on the Virus-VG algorithm [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF] and are used without further adjustment. Steps (2), (3) and (4) are entirely novel. They incorporate a new problem formulation, and based on the solution of this problem, provide a way to estimate contig abundance, as part of an overall efficient alternative to the exponential brute-force routines from [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF]. VG-flow easily scales to data sets of higher complexity and thus provides a clear view towards haplotype reconstruction for mixtures of bacterial strains or even metagenomic data. For further details on algorithm design we refer the reader to the Methods section.

Benchmarking preliminaries

We perform benchmarking experiments where we compare VG-flow to existing methods for fulllength viral quasispecies reconstruction. In these experiments, we make use of the specialized de novo viral quasipecies assembler SAVAGE [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF] for generating a set of strain-specific contigs. We compare performance of VG-flow to Virus-VG [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF], another de novo approach, and to referenceguided viral quasispecies reconstruction tools PredictHaplo [START_REF] Prabhakaran | HIV haplotype inference using a propagating dirichlet process mixture model[END_REF] and ShoRAH [START_REF] Zagordi | ShoRAH: estimating the genetic diversity of a mixed sample from next-generation sequencing data[END_REF]. More recent viral quasispecies assemblers aBayesQR [START_REF] Ahn | aBayesQR: A bayesian method for reconstruction of viral populations characterized by low diversity[END_REF], QSdpR [START_REF] Barik | Qsdpr: Viral quasispecies reconstruction via correlation clustering[END_REF], and PEHaplo [START_REF] Chen | De novo haplotype reconstruction in viral quasispecies using paired-end read guided path finding[END_REF] focus on reconstruction of relatively short genomic regions and were unable to process full-length quasispecies data sets at ultra-deep coverage. We provide the reference-guided methods with a consensus reference genome obtained by running VICUNA [START_REF] Yang | De novo assembly of highly diverse viral populations[END_REF] on the data set. This procedure simulates a de novo setting where the viral agent and its genome may be unknown. Moreover, the consensus reference sequence may be a more accurate representation of the data set under consideration than the standard reference genomes available.

We evaluate all assemblies by comparing the assembled contigs to the ground truth sequences using QUAST [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF]. This assembly evaluation tool aligns the assembled contigs to the true haplotypes, which are provided as a reference, and calculates several standard evaluation metrics. For each assembly, we report the number of contigs, percent target genomes covered, N50, NG50, and error rate. If an assembly consists of only full-length contigs, the number of contigs can be interpreted as the estimated number of strains. Target genome coverage is defined as the percentage of aligned bases in the true haplotypes, where a base is considered aligned if there is at least one contig with at least one alignment to this base. The N50 and NG50 measures reflect assembly contiguity. N50 is defined as the length for which all contigs in the assembly of at least this length together add up to at least half of the total assembly size. NG50 is calculated in a similar fashion, except that the sum of contig lengths is required to cover at least half of the total target length. Error rates are equal to the sum of mismatch rate, indel rate, and N-rate (ambiguous bases). We do not report unaligned bases or misassemblies as we did not encounter any of these.

In addition to the above QUAST assembly metrics, we evaluate strain abundance estimates by comparing estimated values to true strain abundances. For each assembly, let n be the number of true strains and let x i , x i denote the estimated and true abundance, respectively, of strain i. For each ground truth haplotype, the abundance estimates of sequences assigned to this haplotype were summed to obtain the strain abundance estimate x i . We only evaluate abundances for strains that are present in the assembly (i.e. x i > 0) since we cannot expect an assembler to estimate the abundance of a missing strain. Therefore, the true strain abundance values x i are also normalized, taking only the assembled sequences into account. Then, we calculate the absolute frequency error (AFE) and the relative frequency error (RFE) a follows:

AFE = i∈I |x i -x i | |I| , RFE = i∈I |x i -x i | |I| • x i
, where

I = {i ∈ [n] : x i > 0}
VG-flow scales well to bacterial size genomes

Our main goal of designing VG-flow was to enable generation of high-quality assemblies for data sets of higher complexity compared to what is possible with other de novo approaches. While Virus-VG [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF] performs well on the quasispecies benchmarking data sets, the path enumeration step will quickly become too expensive as data sets become more complex: the number of candidate haplotypes is exponential in the number of input contigs. In particular, the number of haplotypes grows exponentially in the genome size, making candidate path enumeration infeasible for larger haplotypes.

In order to explore the limits of VG-flow and to highlight its advantages over Virus-VG, we simulated 28 data sets of increasing complexity using SimSeq [START_REF] St | An illumina paired-end and mate-pair short simulator[END_REF] (2x250 bp paired-end reads, Illumina MiSeq error profile). We created data sets with genomes of increasing size (2500 bp, 5000 bp, 10.000 bp, 20.000 bp, 40.000 bp, 100.000 bp, 200.000 bp) and an increasing number of strains [START_REF] Crotty | RNA virus error catastrophe: direct molecular test by using ribavirin[END_REF][START_REF] Duffy | Why are RNA virus mutation rates so damn high?[END_REF][START_REF] Paten | Genome graphs and the evolution of genome inference[END_REF][START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF]. Each data set has a total coverage of 1000x. Each strain was created by randomly introducing mutations at a mutation rate of 0.5% into a randomly generated nucleotide sequence of the desired length; hence, haplotypes have a pairwise divergence of 1%. For data sets of 2, 4, 6, and 8 strains, the relative strain abundances were set to ratios of 1:2, 1:2:3:4, 1:2:3:4:5:6, 5 and 1:2:3:4:5:6:7:8, respectively. Although the genome sequences are artificial, allowing us to vary genome size and number of strains flexibly, the relative abundances and pairwise divergence reflect plausible real-world, and challenging scenarios in metagenomics [START_REF] Sczyrba | Critical assessment of metagenome interpretation -a benchmark of metagenomics software[END_REF].

In Figure 2 we show VG-flow runtimes as a function of the genome size for a fixed number of strains. As expected, runtime increases as the genome size and the number of strains increase. Even the data sets with a genome size of 200.000 bp are easy to process with VG-flow. Virus-VG, on the other hand, was unable to process any genomes larger than 20.000 bp (2 strains), 5000 bp (4 strains), or 2500 bp (>4 strains).

Remark. Currently, the limiting factor for processing genomes larger than 200.000 bp with VGflow is the pre-assembly step. VG-flow requires pre-assembled strain-specific contigs as input and we use SAVAGE [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF] for this. SAVAGE has proven to produce assemblies of very high quality, but this assembler does not scale well to large genomes. Inspired by results from [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF], we experimented with SPAdes [START_REF] Bankevich | SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing[END_REF] assemblies as input for VG-flow. Although SPAdes does not produce strain-specific contigs as well as SAVAGE, it performs reasonably well and VG-flow is able to build full-length haplotypes from these contigs. Results and further details are shown in the Supplementary Material.

VG-flow outperforms existing tools

We evaluate performance of VG-flow on three simulated viral quasispecies data sets from [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF] and one real HIV benchmark presented in [34], also referred to as the labmix. The simulated data sets are based on true genomic sequences from the NCBI nucleotide database; the characteristics of all data sets (virus type, genome size, number of strains, relative strain abundances, and pairwise divergence) are described in Table 1. All data sets consist of Illumina Miseq reads with an average sequencing depth of 20.000x. For each data set, including the labmix, the true haplotypes and their relative abundances are known.

Table 2 presents assembly statistics for all methods on the three simulated data sets (HCV, ZIKV, and Poliovirus) and Table 3 methods produce full-length haplotypes for all simulated data sets, with much higher assembly N50 and NG50 values than SAVAGE. The improved assembly contiguity comes with only slightly higher error rates compared to the SAVAGE contigs. Table 2 shows that VG-flow builds contigs with even lower error rates than Virus-VG (0.108% versus 0.115% on ZIKV data and 0.036% versus 0.064% on Poliovirus data for VG-flow and Virus-VG, respectively). On the Poliovirus data set we do not only observe a lower error rate for VG-flow compared to Virus-VG, but also a higher target coverage (90.2% for VG-flow versus 80.7% for Virus-VG). The frequency estimation errors (AFE and RFE) in Table 2 show that the increase in assembly accuracy also leads to lower frequency estimation errors.

On real data (Table 3) we observe that VG-flow produces the same number of contigs as Virus-VG, leading to identical target coverage and N50 values; the only difference between the assemblies is a slightly lower NG50 for VG-flow (4608 versus 4642) and a slightly higher error rate (0.535% versus 0.324%). These differences may be explained by the highly uneven coverage of this data set, which affects the contig abundance estimation and hence also the greedy path extraction (see Figure 1). However, the contigs produced by VG-flow are much longer than the input contigs, with the N50 value more than doubled.

Compared to the state-of-the-art for full-length viral quasispecies reconstruction, we notice a clear advantage for VG-flow in terms of target coverage and error rate. Table 2 shows that PredictHaplo and ShoRAH are unable to reconstruct all haplotypes in any of the simulated data sets. The strains that could be reconstructed have higher error rates than VG-flow, as well as much higher frequency estimation errors. On the labmix, PredictHaplo and ShoRAH both achieve 8 a target coverage of 100%. In other words, they assemble each of the five HIV strains at full-length. However, PredictHaplo does so at almost twice the error rate of VG-flow (1.066% for PredictHaplo versus 0.535% for VG-flow) and the ShoRAH assembly has an even higher error rate of 3.910%. Moreover, ShoRAH greatly overestimates the number of strains in all data sets considered. Figure 3 shows the relative frequency estimation errors per method as a function of true abundance per strain, for each of the simulated data sets. Results are divided into bins (binsize=0.05) and average errors are shown. Figure 3 highlights the advantage of de novo methods VG-flow and Virus-VG, which have much smaller relative errors than PredictHaplo and ShoRAH. On the HCV and ZIKV data sets, VG-flow and Virus-VG show nearly identical performance; on the Poliovirus data, we observe a small advantage for VG-flow.

De novo approaches achieve highest precision and recall

In addition to the standard assembly quality metrics presented in the previous section, we analyze relevance of the reported solutions and the amount of similarity between true and reconstructed haplotypes. In the following, we define true positives by their relative edit distance to the corresponding true haplotype (i.e., edit distance divided by alignment length). A contig is considered a true positive if it aligns to a true haplotype with relative edit distance ≤ α; a haplotype is considered correctly reconstructed if at least one contig aligns to it with relative edit distance ≤ α. Figure 4 presents precision, recall, and F-measure per data set. These measures evaluate the number of true positive contigs relative to the total number of contigs (precision), the number of correctly reconstructed haplotypes relative to the total number of haplotypes (recall), and the harmonic average of precision and recall (F-measure = 2*precision*recall/(precision+recall)). We consider various thresholds for the relative edit distance and plot precision, recall, and F-measure as a function of the threshold α. Each of these measures takes values between 0 and 1, with 0 the worst possible score and 1 the best possible score.

We observe that, in general, SAVAGE achieves high values for all three measures already at low relative edit distance. However, SAVAGE only assembles short contigs. VG-flow and Virus-VG show very similar performance, with slightly better values for VG-flow on the ZIKV and Polio data sets, and slightly better performance of Virus-VG on the labmix. All other methods are outperformed by these de novo approaches: PredictHaplo and ShoRAH do not achieve comparable F-measure scores on the simulated data. On the labmix these methods obtain similar scores only at an allowed relative edit distance of 4%, which is nearly as high as the maximal pairwise divergence between strains in this data set.

Runtime and memory usage

The haplotype reconstruction steps used in VG-flow are highly efficient: on the benchmarking data sets presented in Table 1 compared to Virus-VG. However, total runtime for VG-flow is mostly determined by the contigvariation graph construction step, which involves multiple sequence alignment and read mapping. This graph construction step is shared by VG-flow and Virus-VG. Hence, when considering the complete approach on the simulated quasispecies benchmarks, we observe identical runtime and memory usage for VG-flow and Virus-VG: runtime varies between 3.6-12.5 CPU hours and peak memory usage is between 0.6-0.9 GB. Both approaches require as input pre-assembled contigs, for which we used SAVAGE. This de novo assembler constructs an overlap graph from the sequencing reads, which is a very expensive procedure (30.6-276 CPU hours). In comparison, PredictHaplo is faster (2.0-7.4 CPU hours) and ShoRAH is slower (209-814 CPU hours, unable to process the poliovirus data set). However, it is important to realize that all methods except PredictHaplo are able to profit from multithreading, leading to competitive wall clock times on sufficiently large computing clusters. Moreover, constructing a consensus reference genome to be used for reference-guided methods ShoRAH and PredictHaplo also incurs some additional costs (0.07-0.44 CPU hours). More details are presented in the Supplementary Material.

Analysis of an HCV patient sample

In order to demonstrate utility of VG-flow on real data, we ran our method on a patient sample (plasma) of a Hepatitis C virus infection (subtype 1a). This sample is part of a deep sequencing initiative of HCV genomes [START_REF] Hedegaard | High resolution sequencing of hepatitis C virus reveals limited intra-hepatic compartmentalization in end-stage liver disease[END_REF]. It consists of 349268 reads (2x250 bp, Illumina MiSeq), covering the HCV reference genome (NC 004102.1) from position 2296 to 7328 with an average sequencing depth of 34704×. We performed de novo assembly with SAVAGE and obtained 133 contigs varying in length from 152 to 1238 bp, with an N50 value of 472. After running VG-flow on this set of contigs, we obtained 33 contigs with an N50 value of 2342. Among the contigs were 7 full-length haplotypes (>4500 bp), in agreement with the analysis performed by [START_REF] Hedegaard | High resolution sequencing of hepatitis C virus reveals limited intra-hepatic compartmentalization in end-stage liver disease[END_REF] using single genome 10 amplification. The estimated relative frequencies varied between 6.0% and 33.3%. Two of the fulllength haplotypes show a large insertion of 573 bp at position 3546 of the HCV reference genome; this insertion falls into the NS3 gene, which is involved in viral RNA replication through helicase activity. Overall, the 7 full-length haplotypes have 98.8-99.4% pairwise sequence identity, while they share only 93.8-94.7% of their sequences with the HCV reference genome. The overall assembly process took 92 minutes using 12 CPU's (51 minutes for SAVAGE, 41 minutes for VG-flow) and used 1.2 GB of RAM.

Discussion

Many genomic data sets contain mixtures of closely related sequences, such as viral quasispecies or bacterial populations, where the number of haplotypes is generally unknown and relative abundances may differ per haplotype. VG-flow addresses these challenges: we successfully reconstructed haplotypes from several mixed samples, both simulated and real, and obtained highly accurate frequency estimates for each haplotype. VG-flow performs full-length haplotype aware genome assembly, without using existing linear reference genomes, but by constructing variation graphs from pre-assembled contigs. This approach establishes a reference system that allows for analyses that do not suffer from any kind of haplotypespecific mutation-induced biases. We compute abundance estimates for the input contigs, which are of value in its own right. We also enable haplotype reconstruction in polynomial time, with runtimes depending linearly on genome size in practice. We have shown that VG-flow scales well to bacterial sized genomes, hence proving its to contribute also to metagenomic assembly. In benchmarking experiments on simulated viral samples, our method outperformed the state-ofthe-art in full-length viral quasispecies reconstruction in terms of assembly completeness, assembly accuracy, and abundance estimation quality. Finally, we also demonstrated the value of our method on real HIV and HCV data sets.

In view of general benefits of Virus-VG [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF] and the fact that its brute-force solution experiences severe limitations with respect to genome size and the number of contigs, our main goal was to find a polynomial time solution producing high quality assemblies like Virus-VG. Our algorithm did not only achieve decisive speed-ups, but also improved on Virus-VG in terms of assembly accuracy. Virus-VG and VG-flow both aim to reconstruct individual haplotypes from pre-assembled contigs and perform abundance estimation, but VG-flow uses a radically different approach to generating haplotypes. Our results show that the greedy path extraction step in our algorithm selects a subset of possible haplotypes that represents the quasispecies sufficiently well. By limiting the path abundance optimization to this subset of haplotypes, many false haplotypes are excluded from optimization and hence the algorithm gets less confused by false haplotypes.

Interestingly, for some data sets VG-flow was able to improve on the input contigs in terms of target coverage. A possible explanation is that for haplotypes which are not fully represented by pre-assembled contigs, VG-flow is able to reuse contigs from other strains in the same region. This hypothesis is supported by the fact that error rates for VG-flow are higher than for the input contigs constructed using SAVAGE.

The results in Figure 2 show that VG-flow has the potential to process larger genomes. Currently, there are two limiting factors that prevent us from processing bacterial or metagenomic data. First, our method requires as input a collection of pre-assembled, strain-specific contigs. We obtained best results using SAVAGE [START_REF] Baaijens | De novo assembly of viral quasispecies using overlap graphs[END_REF] to generate these input contigs; however, in its current state SAVAGE does not scale well to larger genomes. Experiments using SPAdes to generate input contigs have shown that VG-flow can also generate full-length haplotypes from these assemblies. In fact, any haplotype aware assembler could be used to generate the input contigs, but the quality of the input contigs has a significant impact on the quality of the output.

Another limiting factor is that VG-flow depends on multiple sequence alignment for constructing the contig variation graph. This step can become quite expensive as the number of contigs grows, which could lead to difficulties when processing metagenomic data sets. Nevertheless, we have moved to a much wider range of feasible genome sizes than what was possible before. Note finally, that variation graph construction is part of a current, very active area of research, such that improvements on that end are to be expected [START_REF] Marschall | Computational pan-genomics: status, promises and challenges[END_REF][START_REF] Paten | Genome graphs and the evolution of genome inference[END_REF][START_REF] Garrison | Variation graph toolkit improves read mapping by representing genetic variation in the reference[END_REF].

Note that each of these factors imposes the same limitations to the approach in [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF]. Addressing these challenges would make VG-flow truly capable of processing bacterial or metagenomic data. Therefore, each of the points discussed above provides an interesting starting point for future work.

Conclusions

While multiple approaches to reference-free viral quasispecies assembly have been introduced recently, efficient reconstruction of full-length haplotypes without using a reference genome is a major challenge. Although de novo methods have shown advantages over reference-guided tools, the resulting assemblies often consist of rather short contigs. In this paper, we have proposed VG-flow as an efficient solution to extend pre-assembled contigs into full-length haplotypes, based on variation graphs as reference systems that allow for a bias-free consideration of all haplotypes involved. Benchmarking experiments have shown that VG-flow outperforms the state-of-the-art in viral quasispecies reconstruction in terms of accuracy of haplotype sequences as well as abundance estimates. Moreover, we have shown that our method scales well to bacterial sized genomes, thus proving its potential for processing larger data sets like bacterial mixtures or metagenomic data.

Methods

Variation graphs

Variation graphs are mathematical structures that capture genetic variation between haplotypes in a population [START_REF] Paten | Genome graphs and the evolution of genome inference[END_REF][START_REF] Garrison | Variation graph toolkit improves read mapping by representing genetic variation in the reference[END_REF]. These graphs provide a compact representation of a collection of input sequences by collapsing all shared subsequences between haplotypes. Definition. Let S be a collection of sequences. We define the variation graph V G(S) as a tuple (V, E, P, a). The nodes v ∈ V store sequences seq(v) of nucleotides (of arbitrary length) which appear as a substring of some s ∈ S. The edges (v 1 , v 2 ) ∈ E indicate that the concatenation seq(v 1 )seq(v 2 ) also appears as a substring of some s ∈ S. In addition to nodes V and edges E, a variation graph stores a set of paths P representing the input sequences: for every s ∈ S there is a path p ∈ P (i.e. a list of nodes, linked by edges) such that the concatenation of node sequences equals s. Finally, we store path abundances using an abundance function a : P → R which assigns an absolute abundance value to each path in P . Approach. Following [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF], we distinguish between two types of variation graphs: contig-variation graphs and genome-variation graphs. Let C be a set of pre-assembled contigs and let H be the collection of haplotypes we aim to reconstruct. The contig-variation graph V G(C) = (V C , E C , P C , a C ) organizes the genetic variation that is present in the input contigs and the abundance function a C gives contig abundance values for every input contig. The genome-variation graph V G(H) = (V H , E H , P H , a H ) stores the haplotypes within a population and the abundance function computes haplotype abundances. Constructing a genome-variation graph is the goal of our method; the key idea is to use the contig-variation graph to get there.

Contig-variation graph construction

We construct a contig-variation graph from pre-assembled contigs C using existing techniques for variation graph construction, similar to [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF]. This entails three steps:

(1) Multiple sequence alignment (MSA). We run vg msga [START_REF] Garrison | Variation graph toolkit improves read mapping by representing genetic variation in the reference[END_REF] on the input contigs; the resulting MSA is represented as a graph (V, E, P ).

(2) Compactification. We compactify the graph by contracting any non-branching path into a single node. For every contig, we update the corresponding path p ∈ P such that it stores the path through the compacted variation graph. Thus, we obtain a graph (V C , E C , P C ).

(3) Node abundance computation. We use vg map [START_REF] Garrison | Variation graph toolkit improves read mapping by representing genetic variation in the reference[END_REF] to align the sequencing reads to the compacted variation graph. From these alignments we compute the average base coverage for every node in the graph, also referred to as the node abundance.

Note that the computed node abundances do not yet give us the abundance function a C : P C → R.

To complete the construction of V G(C), we construct a flow network and solve a minimum-cost flow problem as described below.

Flow network construction

We construct a flow network F G = (V, E, c, d), which allows us to compute contig abundances by solving a variant of the minimum-cost flow problem. Network flows are defined on directed graphs, where every edge has a given capacity and receives a certain amount of flow [START_REF] Ravindra | Network Flows: Theory, Algorithms, and Applications[END_REF]. A flow network has a source node and a sink node, which have only incoming and outgoing flow, respectively. For all other nodes, the amount of incoming flow must always equal the amount of outgoing flow, so-called flow conservation. We define our graph as follows.

Nodes. We start by creating a source s and a sink t. Then, we introduce two vertices for every contig c i ∈ C, thus obtaining the vertex set V = {s, t} ∪ {v - i , v + i | c i ∈ C}. Edges. We introduce directed edges (arcs) of three types: contig-arcs, overlap-arcs, and auxiliaryarcs. For each contig c i we add a contig-arc e i : v - i → v + i . For each pair of contigs i and c j there is an overlap-arc e ij from vertex v + i to vertex v - j if a suffix of c i has a non-conflicting overlap with a prefix of c j . In other words, the sequences of i and c j are identical on their overlap. Finally, we add auxiliary-arcs s → v - i for any v - i which has no incoming overlap-arcs, and auxiliary-arcs v + i → t for any v + i which has no outgoing overlap-arcs. Capacities. All edges have infinite capacity. Costs. To every edge e ∈ E, we assign a cost d e where d e =    1, for contig-arcs; -1, for overlap-arcs; 0, for auxiliary-arcs.

The intuition behind this construction is that haplotypes can be found as s -t paths in F G and flow along the edges reflects accumulated haplotype abundances. The edge costs in the flow network allow for the definition of a minimum-cost flow problem that computes contig abundances that are optimal in terms of being compatible with the node abundances in the contig-variation graph, as described in the next section. The construction of the flow network is illustrated in Figure 5.

Contig abundance computation

The problem of estimating contig abundances has applications in metagenomics [START_REF] Li | Microbiome, metagenomics, and high-dimensional compositional data analysis[END_REF] and RNA transcript assembly [START_REF] Conesa | A survey of best practices for RNA-seq data analysis[END_REF]. Existing methods make use of read mapping, either to a reference genome [START_REF] Lindner | Metagenomic abundance estimation and diagnostic testing on species level[END_REF] or to the contigs [START_REF] Fischer | Abundance estimation and differential testing on strain level in metagenomics data[END_REF][START_REF] Bray | Near-optimal probabilistic RNA-seq quantification[END_REF]. Such techniques may cause ambiguous alignments when contigs overlap or share identical sequence. Here, we avoid these issues by mapping reads to the contig-variation graph and solving a flow-like optimization problem. Problem formulation. Candidate haplotypes in the contig-variation graph V G C can be obtained by concatenating overlapping contig subpaths. Therefore, any maximal length path in the variation graph corresponds to an s-t path in F G. We denote by δ + (v) and δ -(v) the set of arcs, respectively, entering and leaving v ∈ V . Recall that V C denotes the set of nodes in V G C and let a u denote the abundance of node u ∈ V C , as computed from the read alignments. For a node u ∈ V C and edge e ∈ E, we write u ∈ e if the contig (or overlap) associated with the contig-arc (or overlap-arc) passes through node u in the contig-variation graph. We define the following flow problem, in which the variables x e decide the amount of flow going through arc e ∈ E:

min u∈V C |a u -{e∈E|u∈e} d e x e | (1) s
.t. e∈δ + (v) x e = e∈δ -(v) x e , ∀v ∈ V \ {s, t} x e ≥ 0, ∀e ∈ E.
Motivation. The objective function evaluates the node abundance errors, defined as the absolute difference of the node abundance and the sum of contig abundance estimates of all contigs whose path passes through the node under consideration. However, contigs belonging to the same haplotype may have overlaps due to conserved regions between haplotypes; we need to avoid doublecounting the contig abundances for nodes corresponding to such overlaps. The edge costs d e ensure that for any pair of overlapping contigs, for any node u ∈ V C in the overlap, the estimated abundance is only added to the sum once.

Solution. The objective function is convex in the flow-variables x e (Supplementary Material, Lemma 3.1), so we have the problem of minimizing a convex function over a set of linear constraints. Such problems can be solved in polynomial time [START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF]. Given a solution to this optimization problem, the flow values on the contig-arcs reflect contig abundance estimates. We use these values to define the abundance function a C on the contig-variation graph. An evaluation of abundance estimates for all simulated data sets is presented in the Supplementary Material. Below, we explain how to use these contig abundances to extract candidate haplotypes for further optimization.

Greedy path extraction

The outcome of the above algorithm gives us a flow value for each edge in the flow network. In the biological context of the problem, we are interested in a decomposition of this flow into a set of s-t paths representing the reconstructed haplotypes. Finding such a flow decomposition can be done in polynomial time, as follows from any constructive proof of the Flow Decomposition Theorem [START_REF] Ravindra | Network Flows: Theory, Algorithms, and Applications[END_REF]. In general, we are interested in a parsimonious solution; that is, a solution with a small number of paths. Finding a decomposition with a minimal number of paths, however, is NP-complete. Many approximation algorithms have been developed for finding a minimum path flow decomposition, e.g. [START_REF] Shao | Theory and a heuristic for the minimum path flow decomposition problem[END_REF][START_REF] Kloster | A practical fpt algorithm for flow decomposition and transcript assembly[END_REF], but these algorithms could not even handle our smallest data set (a mixture of 2 haplotypes of length 2500 bp). Therefore, we resort to other, more efficient means for obtaining a set of haplotypes from the given flow solution [START_REF] Vatinlen | Simple bounds and greedy algorithms for decomposing a flow into a minimal set of paths[END_REF].

We consider a generic greedy heuristic to obtain a selection of candidate paths (Algorithm 1). This approach iteratively selects an s -t path p from the flow network, then updates the flow solution by subtracting the largest possible flow on contig-arcs in p. It terminates when all flow on contig-arcs is below a user-defined threshold for the minimal haplotype abundance. The order in which paths are selected depends on the optimality criterion: we consider maximum capacity paths, minimum capacity paths, and shortest paths-this essentially gives rise to three heuristics. Note that we do not the flow values on overlap-arcs or auxiliary-arcs into account, because we want to avoid any preliminary restrictions on the contig overlaps used.

P cand ← ∅ 4: F G R ← F G 5:
while F G R has at least one s -t path do 6:

Find an s -t path p in F G R that is optimal w.r.t. opt return P cand It varies per data set which optimality criterion gives best results: the maximal capacity criterion extracts paths in order of decreasing abundance, hence leads to paths which are most reliable. However, if a sample contains low-frequency strains, it can be beneficial to select haplotypes in order of increasing abundance (minimum capacity). Since we do not know the composition of the quasispecies beforehand, we combine the results of all three heuristics into one set of candidate haplotypes for further optimization. Earlier work has shown that merging a pool of high quality approximations allows for efficient solutions to well-known optimization problems [START_REF] Cook | Tour merging via branch-decomposition[END_REF][START_REF] Bosman | A solution merging heuristic for the steiner problem in graphs using tree decompositions[END_REF]. We compare performance of our combined approach and the individual greedy heuristics in the Supplementary Material.

Path abundance optimization

Given a collection of candidate haplotypes P cand in the form of paths through the contig-variation graph, the only task remaining is to compute relative abundances for these haplotypes. Although the greedy path extraction algorithm produces preliminary path abundance estimates, these can be improved by the following linear programming approach, also described in [START_REF] Baaijens | Full-length de novo viral quasispecies assembly through variation graph construction[END_REF].

Problem formulation. Let a v denote the abundance of node v ∈ V C , which was computed from the read alignments to V G C . We define variables x p ∈ R ≥0 for p ∈ P cand , representing the estimated abundance for haplotype p, and consider the following optimization problem:

min v∈V C a v - p v x p s.t. x p ≥ 0 ∀ p ∈ P cand . (2)
The objective function is similar in spirit to the objective in Equation ( 1), where abundance estimation errors are evaluated per node in the contig-variation graph. Only now, we compute the absolute difference between the node abundance value and the sum of abundance estimates for all haplotypes passing through this node. This is a convex programming formulation, which can be linearized and solved using an LP solver.

Genome-variation graph construction

Given the candidate haplotypes P cand and the abundance estimates x p for p ∈ P cand , we obtain a final selection of haplotypes H = {p ∈ P cand : x p ≥ m}. Here, m is a user-defined minimal path abundance, by default set to 1% of total sequencing depth. H, we can transform the contig-variation graph V G C into the genome-variation graph V G H , a complete representation of the viral quasispecies.

Figure 1 :

 1 Figure 1: Algorithm overview

Figure 3 :

 3 Figure 3: Abundance estimation results per data set. Abundances were only evaluated for assemblies containing at least 2 full-length haplotypes.

Figure 4 :

 4 Figure 4: Precision, recall, and F-measure per data set.

Figure 5 :

 5 Figure 5: Flow network construction: source (s), sink (t), vertices (v - i , v + i ), contig-arcs, overlap-arcs, and auxiliary-arcs.

Algorithm 1

 1 Greedy path extraction given a flow solution Input: flow network F G, contig-arcs E , flow solution x, min abundance m, optimality criterion Output: a selection of candidate paths P cand 1: function GreedyPaths(F G, E , x, m, opt)

7 :w

 7 ← min e∈p∩E {R e } 8: R ← R -wp 9: F G R ← F G \ {e ∈ E : R e < m} 10: P cand ← P cand ∪ {p} 11:

Table 1 :

 1 presents results on the labmix. The de novo approaches VGflow and Virus-VG both use the contigs obtained with SAVAGE as input. We observe that both Quasispecies characteristics of benchmarking data sets. All data sets consist of Illumina Miseq reads with an average sequencing depth of 20.000x.

	6

Table 2 :

 2 Assembly results on simulated data (Illumina MiSeq, 20.000x coverage). ER = Error Rate (N's + mismatches + indels), AFE = Absolute Frequency Error, RFE = Relative Frequency Error. Frequency errors were only computed for assemblies containing at least 2 full-length haplotypes. * If contigs are full-length, this number reflects the estimated number of strains in the quasispecies.

Table 3 :

 3

	we measured a decrease in haplotype reconstruction time of 9.2-92%
		# contigs * target (%) N50 NG50 ER(%)
	SAVAGE	68	97.9 1026	1450	0.066
	Virus-VG	23	90.6 2130	4642	0.324
	VG-flow	23	90.6 2130	4608	0.535
	PredictHaplo	6	100.0 8825	8825	1.066
	ShoRAH	250	100.0 8775	8775	3.910

Assembly results on the labmix (5-strain HIV mixture, real Illumina MiSeq, 20.000x coverage).

Acknowledgements

Not applicable.

Funding

This work was supported by the Netherlands Organisation for Scientific Research (NWO) through Vidi grant 679.072.309 and Gravitation Programme Networks 024.002.003.

Data simulation

All synthetic data sets were generated using the software SimSeq [START_REF] St | An illumina paired-end and mate-pair short simulator[END_REF] to simulate Illumina MiSeq reads from the genome of interest. In order to obtain realistic sequencing error profiles, we used the MiSeq error profile provided with the software. The genomes used for each data set are listed in the Supplementary Material.

Availability of data and material

Software and analysis scripts are publicly available at https://bitbucket.org/jbaaijens/vg-flow. The synthetic benchmarking datasets analysed during the current study are available at https: //bitbucket.org/jbaaijens/savage-benchmarks. The real HIV data (labmix) is available at https://github.com/cbg-ethz/5-virus-mix. The real HCV data are available in the Sequencing Read Archive under accession number SRR3951347.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

JAB and LS designed the algorithmic components. JAB and AS designed the experiments. JAB implemented the software and performed experiments. JAB wrote the manuscript, with the help of LS and AS. All authors read and approved the manuscript.


De-Li Shi
email: de-li.shi@upmc.fr

Decoding Dishevelled-mediated Wnt signaling in vertebrate early development

Keywords: Wnt signaling, Dishevelled, Wnt/ß-catenin, Wnt/PCP, embryonic axis, gastrulation, neurulation, convergence and extension, morphogenetic movements, mouse, Xenopus, zebrafish

Dishevelled proteins are key players of Wnt signaling pathways. They transduce Wnt signals and perform cellular functions through distinct conserved domains. Due to the presence of multiple paralogs, the abundant accumulation of maternal transcripts, and the activation of distinct Wnt pathways, their regulatory roles during vertebrate early development and the mechanism by which they dictate the pathway specificity have been enigmatic and attracted much attention in the past decades. Extensive studies in different animal models have provided significant insights into the structure-function relationship of conserved Dishevelled domains in Wnt signaling and the implications of Dishevelled isoforms in early developmental processes. Notably, intra-and intermolecular interactions and Dishevelled dosage may be important in modulating the specificity of Wnt signaling. There are also distinct and redundant functions among Dishevelled isoforms in development and disease, which may result from differential spatiotemporal expression patterns and biochemical properties and post-translational modifications. This review presents the advances and perspectives in understanding Dishevelled-mediated Wnt signaling during gastrulation and neurulation in vertebrate early embryos.

Introduction

Wnt signaling plays critical roles in a wide variety of biological process, including embryonic axis formation, cell proliferation, differentiation and migration, polarity establishment, and stem cell self-renewal [START_REF] Steinhart | Wnt signaling in development and tissue homeostasis[END_REF][START_REF] Wiese | Wnt signalling: conquering complexity[END_REF]. Upon stimulation by Wnt ligands, membrane receptors (Frizzled) and co-receptors (LRP5/6, glypican-3/4, ROR, and RYK) assemble into complexes to activate divergent pathways [START_REF] Niehrs | The complex world of WNT receptor signalling[END_REF][START_REF] Green | The role of Ryk and Ror receptor tyrosine kinases in Wnt signal transduction[END_REF][START_REF] Stricker | ROR-family receptor tyrosine kinases[END_REF]. The activity of Frizzled receptors is further controlled by an auxiliary regulatory system involving RSPO1-4, LGR4/5/6 and ZNRF3/RNF43 [START_REF] Jiang | Novel regulation of Wnt signaling at the proximal membrane level[END_REF]Lehoczsky and Tabin, 2018). The canonical Wnt pathway (Wnt/ß-catenin) regulates target gene transcription through stabilization and nuclear accumulation of ß-catenin by inhibition of its destruction complex, consisting of Axin-GSK3ß-APC, whereas the non-canonical Wnt pathway (Wnt/planar cell polarity or Wnt/PCP) is implicated in polarized cellular orientation and asymmetric cell movements through activation of major regulators of the cytoskeleton. Aberrant signaling of both pathways leads to tumorigenesis and metastasis of multiple cancer types, as well as human birth defects [START_REF] Clevers | Wnt/β-catenin signaling and disease[END_REF][START_REF] Butler | Planar cell polarity in development and disease[END_REF][START_REF] Humphries | From instruction to output: Wnt/PCP signaling in development and cancer[END_REF], but how they are regulated in development and disease remains elusive.

Dishevelled (Dvl or Dsh in Drosophila) is a family of proteins that function as common intracellular conductors of both Wnt/ß-catenin and Wnt/PCP pathways [START_REF] Boutros | Dishevelled: at the crossroads of divergent intracellular signaling pathways[END_REF][START_REF] Wallingford | The developmental biology of Dishevelled: an enigmatic protein governing cell fate and cell polarity[END_REF][START_REF] Gao | Dishevelled: the hub of Wnt signaling[END_REF]. Drosophila dsh alleles were first identified in genetic mutants with disruptions of hair and bristle polarity [START_REF] Wallingford | The developmental biology of Dishevelled: an enigmatic protein governing cell fate and cell polarity[END_REF]. Vertebrates possess three highly conserved Dvl genes. Extensive studies in mouse, Xenopus and zebrafish have revealed their critical roles in germ layer specification and morphogenetic movements, which require Wnt/ß-catenin and Wnt/PCP signaling, respectively.

There is accumulating evidence that Dvl isoforms display both distinct and redundant functions [START_REF] Gentzel | Dishevelled paralogs in vertebrate development: redundant or distinct? Front[END_REF]. However, a number of important questions regarding Dvlmediated Wnt signaling during development remain enigmatic [START_REF] Mlodzik | The Dishevelled protein family: Still rather a mystery after over 20 years of molecular studies[END_REF], such as Dvlregulated switch of distinct Wnt pathways, the specific functions of Dvl isoforms in Wnt signaling and development, the post-translational modifications of Dvl functions, and the maternal contributions of Dvl to early developmental events. Fortunately, structure-function and mutational analyses have significantly advanced our understanding of Dvl-regulated Wnt signaling in development and disease. This review focuses on progresses made in this fascinating research field by using complementary vertebrate animal models.

Dvl functional domains in Wnt signaling

Dvl proteins contain several highly conserved domains required for activating different Wnt pathways, including in particular the N-terminal DIX (DIshevelled and aXin) domain, the central PDZ (Post-synaptic density protein-95, Disc large tumor suppressor, Zonula occludens-1) domain, and the C-terminal DEP (Dishevelled, Egl-10 and Pleckstrin) domain (Fig. 1). The DIX domain is involved in Wnt-induced dynamic Dvl homo-and hetero-oligomerization that is important for Wnt/ß-catenin signaling [START_REF] Kishida | DIX domains of Dvl and axin are necessary for protein interactions and their ability to regulate beta-catenin stability[END_REF][START_REF] Kan | Limited dishevelled/Axin oligomerization determines efficiency of Wnt/β-catenin signal transduction[END_REF][START_REF] Ma | Single-molecule dynamics of Dishevelled at the plasma membrane and Wnt pathway activation[END_REF]. The PDZ domain interacts with a conserved KTxxxW motif located immediately after the seventh transmembrane domain of Frizzled receptors [START_REF] Umbhauer | The Cterminal cytoplasmic Lys-Thr-X-X-X-Trp motif in frizzled receptors mediates Wnt/beta-catenin signalling[END_REF][START_REF] Wong | Direct binding of the PDZ domain of Dishevelled to a conserved internal sequence in the C-terminal region of Frizzled[END_REF], and with a wide variety of binding partners that function either as agonists or antagonists of Wnt signaling [START_REF] Wallingford | The developmental biology of Dishevelled: an enigmatic protein governing cell fate and cell polarity[END_REF][START_REF] Sharma | Dishevelled: A masterful conductor of complex Wnt signals[END_REF]. Biochemical and functional analyses suggest that it participates in both Wnt/ß-catenin and Wnt/PCP signaling [START_REF] Habas | Wnt/Frizzled activation of Rho regulates vertebrate gastrulation and requires a novel Formin homology protein Daam1[END_REF][START_REF] Lee | Conformational change of Dishevelled plays a key regulatory role in the Wnt signaling pathways[END_REF]. The DEP domain plays a major role in Dvl membrane recruitment by Frizzled receptors [START_REF] Axelrod | Differential recruitment of Dishevelled provides signaling specificity in the planar cell polarity and Wingless signaling pathways[END_REF][START_REF] Rothbächer | Dishevelled phosphorylation, subcellular localization and multimerization regulate its role in early embryogenesis[END_REF][START_REF] Wong | Structural basis of the recognition of the dishevelled DEP domain in the Wnt signaling pathway[END_REF][START_REF] Pan | Characterization of function of three domains in dishevelled-1: DEP domain is responsible for membrane translocation of dishevelled-1[END_REF][START_REF] Park | Subcellular localization and signaling properties of dishevelled in developing vertebrate embryos[END_REF]. It functions in the Wnt/PCP pathway either with the PDZ domain to activate Rho/ROCK or by direct interaction with Rac to trigger JNK activation. However, more recent evidence suggests that Dvl dimerization triggered by the N-terminal region of DEP domain is required for Wnt/ß-catenin signaling (Gammons et al., 2016a, b;Paclíková et al., 2017).

The C-terminal region beyond the DEP domain also shows evolutionary conservation but displays unique features among Dvl isofroms. Specifically, Dvl3 contains histidine-single amino acid repeats required for Wnt5a-stimulated activation of nuclear factor of activated T cells (NF-AT) and possesses proline-rich domains likely involved in the interaction with other Dvl isoforms [START_REF] Ma | Dishevelled-3 C-terminal His single amino acid repeats are obligate for Wnt5a activation of non-canonical signaling[END_REF][START_REF] Wang | Dishevelled C-terminus: prolyl and histidinyl motifs[END_REF]. The extreme 13 amino acids that are conserved in all Dvl isoforms binds to the third intracellular loop of Frizzled receptors and stabilizes Frizzled-Dvl interaction in Wnt/ß-catenin signaling [START_REF] Tauriello | Wnt/β-catenin signaling requires interaction of the Dishevelled DEP domain and C terminus with a discontinuous motif in Frizzled[END_REF]. Moreover, the last 3 residues represent a type II PDZ-binding (PDZ-B) motif that can occupy the peptide-binding pocket of the PDZ domain, inducing Dvl to adopt a closed conformation and an auto-inhibited state [START_REF] Lee | Conformational change of Dishevelled plays a key regulatory role in the Wnt signaling pathways[END_REF][START_REF] Qi | Autoinhibition of Dishevelled protein regulated by its extreme C terminus plays a distinct role in Wnt/β-catenin and Wnt/planar cell polarity (PCP) signaling pathways[END_REF]. Dvl variants with an opened conformation show efficient membrane recruitment and reduced activity in Wnt/ß-catenin signaling but display increased activity in Wnt/PCP signaling [START_REF] Qi | Autoinhibition of Dishevelled protein regulated by its extreme C terminus plays a distinct role in Wnt/β-catenin and Wnt/planar cell polarity (PCP) signaling pathways[END_REF][START_REF] Harnoš | Dishevelled-3 conformation dynamics analyzed by FRET-based biosensors reveals a key role of casein kinase 1[END_REF]. The function of Dvl C-terminal region in Wnt signaling is further demonstrated in autosomal-dominant Robinow syndrome caused by de novo frameshift mutations in human DVL1 and DVL3 genes, which delete and replace the Cterminal region after the DEP domain [START_REF] Bunn | Mutations in DVL1 cause an osteosclerotic form of Robinow syndrome[END_REF][START_REF] White | DVL1 frameshift mutations clustering in the penultimate exon cause autosomal-dominant Robinow syndrome[END_REF][START_REF] White | DVL3 alleles resulting in a -1 frameshift of the last exon mediate autosomal-dominant Robinow syndrome[END_REF][START_REF] Danyel | Autosomal dominant Robinow syndrome associated with a novel DVL3 splice mutation[END_REF]. In vitro analysis suggests that DVL1 lacking the C-terminal region displays reduced activity in Wnt/ß-catenin signaling [START_REF] Bunn | Mutations in DVL1 cause an osteosclerotic form of Robinow syndrome[END_REF]. These findings suggest an importance of the Cterminus in intra-or inter-molecular interaction, which may be subjected to regulation by other partners to switch pathway specificity. Indeed, recent studies show that casein kinase 1ɛ (CK1ɛ) and NIMA-related kinase 2 (NEK2) function as scaffold proteins and regulate the dynamics of Dvl conformational changes by phosphorylation of the PDZ domain and modulation of its interaction with the extreme C-terminal tail [START_REF] Harnoš | Dishevelled-3 conformation dynamics analyzed by FRET-based biosensors reveals a key role of casein kinase 1[END_REF][START_REF] Hanáková | Comparative phosphorylation map of Dishevelled 3 links phospho-signatures to biological outputs[END_REF].

Dvl dosage effect on morphogenetic movements and cell fate specification during development in mice

The three Dvl genes (Dvl1, Dvl2 and Dvl3) in mice are broadly expressed throughout early development. Extensive analyses of mutant phenotypes have uncovered both unique and redundant functions for these genes. Mice deficient in Dvl1 show reduced social interaction and abnormal sensorimotor gating [START_REF] Lijam | Social interaction and sensorimotor gating abnormalities in mice lacking Dvl1[END_REF]. This abnormal behavior is caused by defective Wnt/ß-catenin signaling that may impair central nervous system functions [START_REF] Belinson | Prenatal β-catenin/Brn2/Tbr2 transcriptional cascade regulates adult social and stereotypic behaviors[END_REF].

Mice deficient in Dvl2 exhibit more severe phenotypes, with defective cardiac morphogenesis, somite segmentation, and neural tube closure [START_REF] Hamblet | Dishevelled 2 is essential for cardiac outflow tract development, somite segmentation and neural tube closure[END_REF]. Dvl3 functions redundantly with Dvl1 and Dvl2 in several processes, including cardiac outflow tract, cochlea and neural tube development [START_REF] Etheridge | Murine dishevelled 3 functions in redundant pathways with dishevelled 1 and 2 in normal cardiac outflow tract, cochlea, and neural tube development[END_REF]. These works reveal a sensitivity of Wnt/PCP signaling to Dvl dosage because most defective phenotypes in Dvl mutants are related to impaired Wnt/PCP signaling, in particular the defective outflow tract morphogenesis [START_REF] Sinha | Disheveled mediated planar cell polarity signaling is required in the second heart field lineage for outflow tract morphogenesis[END_REF]. Thus, Dvl isoforms are critically required for morphogenetic movements. Particularly, they mediate Wnt/PCP in CE movements during neurulation. Dvl2 plays a predominant role in neural tube closure, but Dvl1 and Dvl3 are also involved in this process [START_REF] Wang | Dishevelled genes mediate a conserved mammalian PCP pathway to regulate convergent extension during neurulation[END_REF]. The functional importance of Dvl isoforms in neural tube formation has been confirmed by the identification of rare mutations in all three human DVL genes, which disturb normal functions of DVL isoforms in non-canonical Wnt signaling and cause neural tube defects [START_REF] De Marco | Genetic analysis of disheveled 2 and disheveled 3 in human neural tube defects[END_REF][START_REF] Liu | DVL mutations identified from human neural tube defects and Dandy-Walker malformation obstruct the Wnt signaling pathway[END_REF]. By comparison, low levels of Dvl expression from a single allele may be sufficient to normally support those developmental processes triggered by the Wnt/ß-catenin pathway [START_REF] Soares | Functional studies of signaling pathways in peri-implantation development of the mouse embryo by RNAi[END_REF][START_REF] Wynshaw-Boris | Dishevelled: in vivo roles of a multifunctional gene family during development[END_REF]. However, deletion of all six Dvl alleles causes absence of mesoderm gene expression and mesoderm formation that are dependent on Wnt/ßcatenin signaling [START_REF] Ngo | Deletion of the Dishevelled family of genes disrupts anterior-posterior axis specification and selectively prevents mesoderm differentiation[END_REF].

Distinct and redundant Dvl functions during Xenopus development

Dvl function in vertebrates was first studied during Xenopus development. Overexpression of Dvl2 (Xdsh) in the ventral region of early embryos induced the formation of a complete secondary axis reminiscent of activation of maternal Wnt/ß-catenin signaling [START_REF] Sokol | Dorsalizing and neuralizing properties of Xdsh, a maternally expressed Xenopus homolog of dishevelled[END_REF]. However, dorsal overexpression of Xdd1, a truncated form of Dvl2 that lacks the PDZ domain and interferes with Wnt/ß-catenin signaling triggered by Wnt ligands, did not affect dorsoventral axis formation [START_REF] Sokol | Analysis of Dishevelled signalling pathways during Xenopus development[END_REF]. Because maternal Wnt/ß-catenin signaling is required for dorsal fate specification by activating the transcription of target genes in the Spemann organizer [START_REF] Carron | Specification of anteroposterior axis by combinatorial signaling during Xenopus development[END_REF], the absence of an inhibitory effect by Xdd1 implies that Dvl function may be dispensable for the activation of maternal Wnt/ß-catenin signaling. Consistently, simultaneous depletion of maternally expressed Dvl2 and Dvl3 from oocytes did not affect the expression of maternal Wnt/ß-catenin target genes and the formation of dorsal axis [START_REF] Tadjuidje | The functions of maternal Dishevelled 2 and 3 in the early Xenopus embryo[END_REF]. However, it is possible that low levels of Dvl proteins are still present in the oocytes due to incomplete depletion of maternal Dvl mRNA. Thus, the requirement of maternal Dvl for dorsal axis formation in Xenopus requires a complete loss-of-function study. Nevertheless, a recent study suggests that activation of maternal Wnt/ß-catenin pathway and formation of dorsal axis may be achieved through a Dvl-independent mechanism [START_REF] Yan | Maternal Huluwa dictates the embryonic body axis through β-catenin in vertebrates[END_REF]. During organogenesis, it seems that Dvl isoforms display less functional redundancy in developmental processes that involve Wnt/ß-catenin signaling, which may be due to their differential expression patterns. For example, Dvl1 and Dvl2 are required for neural crest cell specification and somite segmentation, while Dvl3 maintains gene expression in the myotome [START_REF] Gray | Diversification of the expression patterns and developmental functions of the dishevelled gene family during chordate evolution[END_REF].

Maternal Dvl2 and Dvl3 likely display distinct and redundant functions during CE movements. In these coordinated processes, lateral cells converge toward the dorsal region to narrow the germ layers, while dorsal midline cells undergo mediolateral intercalation by polarized protrusive behaviors to lengthen the embryo along the anteroposterior axis [START_REF] Keller | Convergent extension in the amphibian, Xenopus laevis[END_REF]. Previous studies show that Dvl regulates CE movements through Wnt/PCP signaling [START_REF] Tada | Xwnt11 is a target of Xenopus Brachyury: regulation of gastrulation movements via Dishevelled, but not through the canonical Wnt pathway[END_REF][START_REF] Djiane | Role of frizzled 7 in the regulation of convergent extension movements during gastrulation in Xenopus laevis[END_REF][START_REF] Wallingford | Xenopus Dishevelled signaling regulates both neural and mesodermal convergent extension: parallel forces elongating the body axis[END_REF][START_REF] Habas | Wnt/Frizzled activation of Rho regulates vertebrate gastrulation and requires a novel Formin homology protein Daam1[END_REF][START_REF] Habas | Coactivation of Rac and Rho by Wnt/Frizzled signaling is required for vertebrate gastrulation[END_REF]. Moreover, inhibition of maternal Dvl2 or Dvl3 function suggests that they exhibit a non-redundant but an additive effect on CE movements [START_REF] Tadjuidje | The functions of maternal Dishevelled 2 and 3 in the early Xenopus embryo[END_REF]. More recent works reveal that Dvl1, but not Dvl2 or Dvl3, activates the Wnt/Ca 2+ pathway (another branch of the non-canonical pathway) during CE movements [START_REF] Gentzel | Distinct functionality of dishevelled isoforms on Ca2+/calmodulin-dependent protein kinase 2 (CamKII) in Xenopus gastrulation[END_REF], further supporting the distinct functions of Dvl isoforms in morphogenetic movements. As in mice, Dvl function is also required for neural tube closure by regulating CE movements of the midline and coordinating polarity among epithelial cells [START_REF] Wallingford | Neural tube closure requires Dishevelled-dependent convergent extension of the midline[END_REF][START_REF] Seo | Bimodal regulation of Dishevelled function by Vangl2 during morphogenesis[END_REF]. In addition, Dvl2mediated activation of Rac1 through the DEP domain also controls the protrusive activity of neural crest cells during migration [START_REF] Kratzer | The Rho guanine nucleotide exchange factor Trio is required for neural crest cell migration and interacts with Dishevelled[END_REF]. Because cellular polarization is tightly dependent on the asymmetric activation of Wnt/PCP signaling, both reduced and increased levels of Dvl disrupts asymmetric movements. However, the cellular behaviors are completely different. Increasing the activity of Dvl perturbs cell polarity by randomizing the formation of cellular protrusions, whereas reducing the activity of Dvl prevents cellular protrusions [START_REF] Wallingford | Dishevelled controls cell polarity during Xenopus gastrulation[END_REF][START_REF] Cheng | Leucine repeat adaptor protein 1 interacts with Dishevelled to regulate gastrulation cell movements in zebrafish[END_REF].

Maternal contributions of Dvl proteins to axis patterning and morphogenetic processes during zebrafish development

The zebrafish genome contains five Dvl paralogs: dvl1a, dvl1b, dvl2, dvl3a, and dvl3b. In the early embryos, dvl2 and dvl3a are maternally expressed and represent about 98% of the total pool, whereas the transcript levels of the other Dvl genes are negligible [START_REF] Harvey | Identification of the zebrafish maternal and paternal transcriptomes[END_REF]. Knockout of Dvl genes reveal both distinct and redundant functions in embryonic axis specification and morphogenetic movements [START_REF] Xing | Mutational analysis of dishevelled genes in zebrafish reveals distinct functions in embryonic patterning and gastrulation cell movements[END_REF]. Maternal-zygotic mutants for dvl1a, dvl1b, dvl3a, and dvl3b are phenotypically normal and fertile. In sharp contrast, maternalzygotic dvl2 mutants display strongly impaired CE movements during gastrulation and develop severe craniofacial defects with a "bulldog" facial phenotype, reminiscent of impaired Wnt/PCP signaling in midline structures [START_REF] Kimmel | Specification and morphogenesis of the zebrafish larval head skeleton[END_REF]. Zygotic dvl2 mutant embryos are essentially normal, but only about half of them can survive to adulthood, and all male individuals show absence of courtship behavior. This suggests that there may be defects in central nervous system functions as mice Dvl1 mutants [START_REF] Lijam | Social interaction and sensorimotor gating abnormalities in mice lacking Dvl1[END_REF]. The highest level of dvl2 expression during early development may explain at least partly the most severe phenotypes of dvl2 mutants.

The specification of dorsal axis in zebrafish also requires maternal ß-catenin signaling to trigger the expression of organizer genes [START_REF] Kelly | Maternally controlled (beta)-catenin-mediated signaling is required for organizer formation in the zebrafish[END_REF][START_REF] Bellipanni | Essential and opposing roles of zebrafish beta-catenins in the formation of dorsal axial structures and neurectoderm[END_REF][START_REF] Fuentes | The maternal coordinate system: Molecular-genetics of embryonic axis formation and patterning in the zebrafish[END_REF]. However, the involvement of upstream regulators has not been conclusively established. Importantly, maternal-zygotic dvl2 and dvl3a double mutants, which were unresponsive to stimulation by Wnt ligands, show normal specification of dorsal cell fate, suggesting that components of Wnt signaling upstream of ß-catenin may be dispensable for its stabilization. The activation of maternal Wnt/ß-catenin signaling independent of Dvl activity is further confirmed in zebrafish huluwa mutants. Maternal depletion of huluwa impairs Wnt/ßcatenin signaling and causes loss of dorso-anterior structures. Mechanistically, Huluwa protein accumulates in the cell membrane at the dorsal region and functions independently of Wnt ligands and Frizzled receptors to promote tankyrase-mediated degradation of Axin, thereby stabilizing ß-catenin [START_REF] Yan | Maternal Huluwa dictates the embryonic body axis through β-catenin in vertebrates[END_REF].

Mutational analyses of Dvl functions also reveal a major contribution of maternal Dvl to zygotic events and confirm the importance of Dvl dosage in Wnt/PCP signaling (Fig. 2). This dosage effect is particularly reflected by the requirement of Dvl2 and Dvl3a for CE movements.

Although Dvl2 plays a predominant role, Dvl3a exerts a strong synergistic effect on the loss of Dvl2 function, and progressive removal of Dvl2 and Dvl3a maternal or zygotic products increasingly aggravates CE defects and reduces the elongation of anteroposterior axis.

Furthermore, zygotic dvl2 and dvl3a double mutants only display a shortened body length.

However, maternal-zygotic dvl2 and dvl3a double mutants show most strongly impaired CE movements and completely lack axis extension. They also develop severe trunk and posterior deficiencies associated with down-regulation of zygotic Wnt/ß-catenin target genes [START_REF] Xing | Mutational analysis of dishevelled genes in zebrafish reveals distinct functions in embryonic patterning and gastrulation cell movements[END_REF]. Because zygotic Wnt/ß-catenin signaling has an opposite effect with respect to maternal Wnt/ß-catenin signaling and functions to specify the posterior region [START_REF] Carron | Specification of anteroposterior axis by combinatorial signaling during Xenopus development[END_REF], these findings highlight the importance of maternal Dvl in setting up zygotic morphogenetic and patterning processes. They support the view that maternally expressed gene products perform essential functions after zygotic genome activation [START_REF] Marlow | Setting up for gastrulation in zebrafish[END_REF][START_REF] Solnica-Krezel | Maternal contributions to gastrulation in zebrafish[END_REF]. Thus, both maternal and zygotic Dvl dosages are important for proper cell movements and embryonic axis patterning that occur during gastrulation.

Perspectives

Dvl conserved domains in Wnt pathways have attracted much attention. Although the function of DIX, PDZ and DEP domains is relatively understood, how they cooperate to switch pathway specificity remains elusive. Detailed analysis of other conserved domains or isoform-specific regions, such as the basic region preceding the PDZ domain, the proline-rich region and histidine-single amino acid repeats in the C-terminal region beyond the DEP domain, and the extreme C-terminus, may provide insights into Dvl-mediated signal transduction. Because Dvl post-translational modifications, in particular phosphorylation and ubiquitination, and Dvl interaction partners are important for subcellular localizations and specific functions of Dvl proteins [START_REF] Sharma | Dishevelled: A masterful conductor of complex Wnt signals[END_REF][START_REF] Harrison | Molecular mechanisms mediating asymmetric subcellular localisation of the core planar polarity pathway proteins[END_REF], it is of interest to understand how these modulate Dvl activity and dictate signaling outcomes in key developmental processes. Indeed, dysregulation of Dvl phosphorylation impairs both Wnt/ß-catenin and Wnt/PCP signaling during zebrafish and Xenopus embryogenesis [START_REF] Shimizu | Hipk2 and PP1c cooperate to maintain Dvl protein levels required for Wnt signal transduction[END_REF][START_REF] Rauschenberger | The phosphatase Pgam5 antagonizes Wnt/β-Catenin signaling in embryonic anterior-posterior axis patterning[END_REF]. The specificity of Dvl isoforms also merits investigations because tissue-specific expression patterns and differential biochemical properties may contribute to their particular functions. Another intriguing question is the sensibility of Wnt/PCP signaling, but not Wnt/ß-catenin signaling, to Dvl dosage, which is observed in all vertebrates. It suggests that Wnt/PCP-dependent developmental processes critically require Dvl function to activate downstream effectors. Consequently, moderate diminution of Dvl dosage could significantly affect polarized cellular behaviors and cell polarity. By comparison, ß-catenin may be stabilized independently of upstream Wnt signaling components, and target genes of Wnt/ß-catenin signaling may be regulated by other factors, at least during dorsal fate specification in zebrafish and Xenopus [START_REF] Li | Direct regulation of siamois by VegT is required for axis formation in Xenopus embryo[END_REF][START_REF] Yan | Maternal Huluwa dictates the embryonic body axis through β-catenin in vertebrates[END_REF].

Thus, tissue morphogenesis regulated by Wnt/PCP signaling is more sensitive to Dvl dysfunction, as a result, many human congenital disorders, such as neural tube defects and Robinow syndrome, are associated with mutations in DVL genes.

Fig. 1 .

 1 Fig. 1. Structure-function relationship of conserved Dvl domains in the activation of Wnt/ß-catenin and Wnt/PCP signaling (see text for details).

Fig. 2 .

 2 Fig. 2. Maternal and zygotic Dvl dosage effect on embryonic axis elongation in zebrafish embryos at 5 days post-fertilization (Xing et al., 2018). (A) Wild-type sibling embryo. (B) Zygotic dvl2 mutant. (C) The "bulldog" facial phenotype of maternal-zygotic dvl2 mutant with cyclopia. (D) Maternal-zygotic dvl3a mutant. (E) Triallelic dvl2 and dvl3a mutant. (F) Zygotic dvl2 and dvl3a double mutant. Scale bar (A-F): 400 µm.

Acknowledgments

This work was supported by the Centre National de la Recherche Scientifique (CNRS) and the Sorbonne University.

Author contributions

DLS performed the literature analysis, prepared the figures, and wrote the manuscript.

Conflicts of interest

The author has no conflict of interest to declare.