Origin of the bidirectionality of cerebrospinal fluid flow and impact on long-range transport between brain and spinal cord

The circulation of cerebrospinal fluid (CSF) plays pivotal roles for body axis formation and brain development. During embryogenesis, CSF is rich in particles and proteins and flows bidirectionally in the central canal. The origins of bidirectional flow and its impact on development are unknown. Experiments combined with modeling and simulations demonstrate that the bidirectionality of CSF flow is generated locally by caudally-polarized motile cilia confined to the ventral wall of the central canal. Such active bidirectional flow of the CSF accelerates the long-range transport of particles propagating rostrally and caudally. In addition, spontaneous muscle contractions increase local CSF flow and consequently enhance long-range transport of extracellular lipidic particles. Focal ablation of the channel connecting brain ventricles to the central canal reduces embryo length, indicating that long-range transport contributes to embryonic growth. Our study also demonstrates that at this early stage, motile cilia ensure the proper formation of the central canal.


Introduction
Precise control of flow in biological systems is essential to transport critical molecules to the different cells in the organism. In order to respond to the cells demands, passive diffusion in tissues imposes that every cell in the organism is close enough to a supply source. For organisms larger than a few hundreds of microns, diffusion alone is too slow to achieve supply of nutrients to all cells in the body. Therefore, organisms have developed a large variety of flows to accelerate and orient the transport of key molecules. Blood flow is the main carrier of glucose and oxygen in most organisms and is mostly driven by the pressure gradient caused by the heart pump 1 . Another large category of biological flows, including intestine, esophageal or lymph flows, relies on peristalsis, where muscles contractions drive the fluid in the direction of propagation of the contraction wave [2][3][4][5][6] . Another major class of fluid flows relies on the asymmetric beating of motile cilia or flagella, which are efficient motors to generate directional flows at small scales and to mix fluids in biological systems [7][8][9][10] . For instance, ciliadriven flows were demonstrated to be critical for left/right asymmetry in the developing embryos of many species including humans 11 , in the Kupffer's vesicle in zebrafish 10,12 , and in the node in mammals [13][14][15] . Ciliadriven flows can also be observed in the pronephros 16,17 , in the respiratory tract 18 , and in the nasal cavity 19 . The motility of cilia in brain ventricles is also correlated with the cerebrospinal fluid (CSF) flow in zebrafish 17,20,21 , African clawed frog 22,23 and rodents 24 .
Understanding CSF circulation throughout the body is crucial as this biological fluid plays numerous roles in the development of brain and spinal cord. The CSF provides the hydromechanical protection of the brain 25 , enables the transport of nutrients and waste 25,26 , exosomes [27][28][29][30][31] , and signaling molecules impacting neurogenesis and neuronal migration in the brain [32][33][34][35] . The CSF circulation has recently been linked to body axis formation in the embryo and spine morphogenesis in juveniles/adults 21,36,37 . During embryogenesis, mutations affecting ciliogenesis, cilia polarity or motility lead to defects of transport along the rostrocaudal axis and body axis formation 36,38 .
Historically, it has been generally assumed that CSF flow along the walls of all cavities must be driven by motile cilia. In the brain ventricles, the direction of beating cilia correlates indeed with the direction of CSF flow 20,24 . However, in the long cylinder of the central canal, the dynamics of CSF flow is strikingly different than in the brain ventricles. We recently showed in 24 hours post fertilization (hpf) zebrafish embryos that the CSF flows bidirectionally 31 , 39 . To our knowledge, these observations report for the first time a constant bidirectional flow in a single channel in biological systems. Interestingly, bodyaxis defects observed in mutants with defective cilia appear by 30 hpf, i.e. few hours before cilia paving the walls of the brain ventricles become motile and a directed CSF flow emerges in the brain ventricles 20,40 . It suggests that proper dynamics of CSF flow in the central canal, not in the brain ventricles, is critical for embryonic body axis formation.
Since the mechanisms leading to the bidirectionality of CSF flow in the central canal were unknown, we developed an automated method for quantifying CSF flow in the central canal of zebrafish embryos. We show that motile cilia are critical to form the lumen of the central canal, as the canal collapses in mutants with defective cilia. We combined experimental, numerical and theoretical approaches to test the role of the spatiallyasymmetric distribution of motile cilia in the generation of a bidirectional flow. We developed a general and simple hydrodynamic model to account for the contribution of cilia, and show that this model can be applied as a tool in many ciliadriven flows in confined environments. Additionally, we show experimental and theoretical evidence for an intricate relation between local flow and longrange transport of particles in the CSF. We demonstrate that this bidirectional CSF flow in the central canal is enhanced by muscle contractions, which accelerate bidirectional transport of particles both rostrally and caudally. We show that the CSF flows unidirectionally from the diencephalic/mesencephalic ventricle to the entrance of the central canal via a long and thin funnel. In vivo twophoton ablation of this funnel leads to a reduction in growth, suggesting that the exchange of critical particles between brain ventricles and central canal is necessary for optimal embryogenesis.

Quantification of CSF flow in the central canal
Understanding the mechanisms determining CSF flow in the central canal requires an in vivo quantification and a theoretical framework model with minimal parameters to explain it. To establish quantitative and reproducible flow profiles along the dorsoventral (DV) axis of the central canal (CC), we develop an improved automated analysis inspired by our previous manual kymograph measurements 31,39 . By coinjecting 3,000 molecular weight (MW) red dextran with 20 nmdiameter green fluorescent beads in the diencephalic/mesencephalic ventricle (DV), we measure the volume of cavities and particle trajectories. The dextran rapidly diffuses over all fluidfilled cavities, including the CC ( Figure 1a ), enabling its visualization and the assessment of the injection quality. Fluorescent nanobeads show random Brownian motion superimposed to a smooth displacement along the streamlines ( Supplementary Movie S1 ). Bead trajectories reveal the bidirectionality of CSF flow along the rostrocaudal axis and the presence of a few regions of recirculation, commonly referred to as vortices ( Figure 1b ). The flow analysis relies on the permutation of the 3D (2D in space and time) stack depicting trajectories of fluorescent beads ( Figure  1c1 ) to form kymographs at various DV positions in the CC ( Figure 1 c2, c3, c4, and Supplementary Movie S2 ). The typical velocity profile is bimodal with rostrocaudal flow in the ventral side and caudorostral flow in the dorsal side, and reverses close to the center of the CC (normalized DV position = 0.53, Figure 1d from 110 wild type (WT) embryos each imaged in 2 rostrocaudal positions). Automated kymographs show better reproducibility compared to particle tracking velocimetry (PTV) or particle imaging velocimetry (PIV), as discussed in Methods section. Velocity profiles show ventrally a maximal speed of 4.78 ± 0.79 μm.s 1 (DV = 0.29 on average) and dorsally a minimal speed of 4.80 ± 0.82 μm.s 1 (DV = 0.82 on average) ( Figure 1d1 and 1d2 ) and is nearly symmetrical (DV =0.53 on average). Similarly, the average speed throughout the CC is zero (0.02 ± 0.15 μm.s 1 ), refuting the hypothesis of a global pressure gradient as the unique driver of the flow in the CC.

Central canal geometry and properties of the motile cilia
To build a theoretical model of CSF flow, we quantify essential parameters such as the canal geometry and the local cilia dynamics. The central canal geometry, as estimated by Dextran injections, corresponds to a cylinder ( Figure 2a1, 2a2 ) of diameter 8.9 ± 0.9 μm (129 WT embryos, Figure 2a3 ) with a very small dorsal extension showing low fluorescence level in the uppermost position, as previously observed with antibody staining for tightjunctions 41 . The dynamics of motile cilia in the CC are investigated in 30 hpf Tg(βactin:Arl13bGFP) embryos in which GFP labels all cilia 31,39,42 . Density and motility of cilia are higher in the ventral CC as previously reported 31,39,42 . Movies recording the cilia beating illustrate the high density of cilia and diversity in ciliary length and beating frequency ( Figure 2b , Supplementary Movie S3 and Supplementary Figure S1 ). While some cilia beat at high frequency (at least 45 Hz, see Methods ) with a single frequency pattern, others beat at lower frequency (~13 Hz) with more complex and irregular patterns. An automated quantification of cilia beating frequency using temporal Fourier transform on each pixel ( Figure 2b1 ) reveals distinct spots of constant frequencies ( Figure 2b2 ), likely corresponding to the displacement of a single cilium. For each spot, we extract the main frequency, orientation, length and height ( Figure 2c14 ), indicating that on average cilia beat at 37.3 Hz with a caudal tilt of 65.0°towards the tail (median of the absolute value, see Methods ), and are 6.3 μm long, thereby occupying a 2.9 μm height in the ventral canal, representing a bit less than half of the central canal. The caudal tilt is consistent with previous reports 39,42 . The measured frequencies are higher than the 1215 Hz previously reported in zebrafish central canal 17 (see Methods explaining the discrepancies), but are consistent with beating frequencies of ependymal cilia located in the brain ventricles of rats 43 and zebrafish larvae 20,43 .

The asymmetric distribution of beating cilia generates a bidirectional flow
The CC can be considered as a cylinder ( Figure 2a ) containing numerous motile cilia along the ventral wall ( Figure 2b ). Modeling the effect of each individual cilium beating would be extremely challenging. Since the CSF bidirectional flow does not vary over time and over rostrocaudal axis ( Figure 1c4 ), we build a twodimensional model that homogenizes cilia contribution as a constant force. The model divides the CC in two regions of equal thickness ( Figure 3a1, To solve this system of equations, one needs a third relation given by the constraint of zero flux across the diameter of the central canal as observed in vivo : With a unique fitting parameter chosen as 0.5, (leading to N/m 3 ) we obtain a ⍺ 000 f v = 4 symmetric bimodal flow with extreme velocities equal to ±5 μm/s ( Figure 3a3 ), and an axial pressure gradient N/m 3 . We additionally solved the complete Stokes equation dP dz 000 / = 2 numerically, using FiniteElements Methods simulations (COMSOL software) to confirm the hypotheses of the theoretical model. Theoretical and numerical predictions are in good agreement with 57 velocity profiles measured in vivo ( Figure 3b ). The model is further refined to account for the spatial inhomogeneity in the cilia distribution in the ventral region, that we design as spots of constant volume force f v regularly distributed along the ventral wall of CC with a period ( Figure 3c1 ). With this w "scarce cilia" model, numerical simulations show that recirculation regions appear between these spots ( Figure 3c2 ). The amplitude of the recirculation regions is more pronounced as increases. For large passive regions of width larger than the diameter of the CC ( ), w w > d all the streamlines are located within the vortices. Counterintuitively, this model predicts that vortices may originate not from cilia dynamics, but rather from the local absence of motile cilia in the ventral side on a distance larger than d. Note that the presence of passive regions decreases the pressure in the CC proportionally to the fraction of passive regions ( Supplementary Figure S2a ). We can further generalize these results by showing that differences in beating frequency between active neighboring cilia lead to the emergence of recirculation regions ( Supplementary Figure S2d ).

Cilia activity controls CSF longrange transport, CC structure, and local flow
To verify that asymmetric cilia beating is solely responsible for the bidirectional flow i n vivo , we injected beads in the diencephalic ventricle of several mutants for which either ciliogenesis or ciliary motility is affected. Consistently with previous reports 36,39 , longrange transport of beads down the central canal is altered in mutants with defective cilia ( Figure  4a1, 4a2 ). In paralyzed mutant embryos several hours post injection, no particles entered the CC, as shown here in traf3ip tp49d/tp49 , commonly referred to as elipsa 44 . Besides, we notice that the CC geometry is dramatically altered in elipsa as well as in dnaaf1 tm317b/tm317b , (referred to as lrrc50 45 ), in the newlygenerated mutant allele foxj1a nw2 (see Methods and characterization in Supplementary Figure S3 ) , and in cfap298 tm304/304 , referred to as kurly 46 ( Figure 4b,c ). While lrrc50 embryos show complete loss of cilia motility 45 , foxj1a nw2 mutants retain normal primary cilia but show a decreased density of glutamylated cilia in the central canal ( Supplementary Figure S3 ). Glutamylated cilia in the CC of WT embryos are caudallypolarized and ventrallylocated and most likely correspond to motile cilia ( Supplementary Figure S3 e1,e2 ). The severe reduction of glutamylated cilia foxj1a nw2 mutants is consistent with the previous report showing the necessity of foxj1 genes for formation of motile cilia 55 . In the mutant elipsa 39,44 , cilia partially degenerate. The mutant allele cfap tm304 is characterized by a thermosensitive loss of cilia motility 46 . In all four mutants with defective cilia that we tested here, the CC forms but is flattened along the DV axis and tortuous compared to WT siblings ( Figure 4b1, b2, 4c) . In contrast, we measure no difference in the CC diameter in the mutant scospondin icm15/icm15 , which also exhibits a curled down phenotype but with normal cilia properties and activity 39 , suggesting that cilia integrity and motility is a critical parameter to open the CC. It also suggests that previous observations describing defects of transport in mutants with defective cilia 36 could be simply explained by the collapse of the CC, and not only by CSF flow defects. Because exogenous beads are not transported down the CC in paralysed mutants, one cannot directly test the role of motile cilia on generation of the local CSF flow. Instead, we investigate whether letting mutant embryos with ciliary defects spontaneously twitch would facilitate the transport of beads. Indeed, we find that for mutants with a moderate reduction of the CC diameter, such as elipsa and kurly , some beads can enter the CC. These beads exhibit Brownian motion and some exhibit a residual flow with numerous recirculation zones ( Supplementary Movie S4 ), in good agreement with our "scarce cilia" model.
Our model also predicts that the bidirectional CSF flow in the CC is generated locally and does not rely on a global pressure gradient generated in ventricles. This is confirmed by multiple observations: i) even after large regions without cilia and local flow, the CSF flow in more caudal regions can show a bidirectional flow, as long as motile cilia beat locally in the ventral region (see Figure 3c ). ii) When some portion of around 200 μm of the CC is isolated by using photoablation to block flow on both sides ( Figure 4d1 ), the canal remains wideopen and the CSF still flows bidirectionally ( Supplementary Movie S5 ) with a similar velocity profile ( Figure 4d2 ).

Muscle contraction transiently modifies CSF flow
So far, we mainly restricted our study to paralyzed zebrafish embryos. However, muscle contractions are also at the origin of several biological flows 47 , and were recently shown to create a massive transient flow in the brain ventricles of zebrafish embryos 20 . Consequently, we investigate in details the CSF flow and transport on spontaneously twitching embryos restrained in agar where contraction can be identified from the notochord displacement (white arrow) on the transmitted image ( Figure 5a ). Following each contraction, rapid displacements of CSF (too fast to be properly captured here) are observed in either caudal or rostral direction, followed by a slower counterflow in the opposite direction, whose velocity decays until the original bidirectional flow is recovered ( Figure 5b . However, we do not find any obvious correlation between the flow direction and our estimates of the contraction direction or strength. Body contractions are also shown to enhance transport of beads in mutants where cilia are defective ( Supplementary Movie S4 ). As no significant difference in the CC diameter is detected during the muscle contractions, we propose that the muscledependent enhancement of CSF flow is not due to a local pinching of the canal, but rather to a change in brain ventricles volume recently observed 20 .

Bidirectional longrange transport in the central canal is accelerated by the flow
Now that the CSF flow has been fully characterized, its ability to transport particles on long distances needs to be demonstrated. The longrange transport follows an apparent diffusive law, as the propagation front of nanoparticles progresses with the square root of time 36 and depends on particle size ( Figure 6a1, 6a2 ). The propagation front is also accelerated by muscle contractions ( Figure 6a3 ). In practice, 20nm particles injected in the DV travel through the entire CC in around two hours ( Figure 6a2 ), far from the timescales expected for both passive diffusion (more than one day) and sole active transport of 5 μm/s (less than 10 minutes). This suggests that the transport of particles should be modelled by including an advection term in the diffusion equation. In a bidimensional canal ( Figure 3a2 ), the transport equation can be written as: where is the local concentration of particles, and , and are respectively the time c ∂ t ∂ z ∂ y and spatial derivatives in the rostrocaudal ( ) and dorsoventral ( ) directions. is the z y (y) v velocity profile ( eqs. 1, 2 ), and is the diffusion coefficient of particles in CSF. For a D spherical particle of radius , the StokesEinstein relation gives: r where is the temperature and is Boltzmann's constant.
T k B As seminally described by Taylor 48 and Aris 49 , the diffusivity of solutes is amplified in presence of velocity gradients, due to a coupling between pure diffusion and flowinduced transport. They demonstrated that on long distances, the transport still follows a diffusive process, with an effective diffusivity enhanced by the shear flow. Following analogous developments of Bruus 50 , we perform the theoretical derivation of the effective diffusivity of particles submitted to a bidirectional flow ( Supplementary Information ), which leads D eff to: where is the dimensionless Péclet number, comparing the characteristic times d D  concentration profile is recovered, with a wider extent in presence of a flow ( Figure 6b ). This diffusive transport is numerically obtained for three different maximal speeds of the flow: 0, 6 μm/s, and 20 μm/s ( Figure 6c1 ), and shows that transport is accelerated by the flow speed, as predicted by our theoretical derivation. Another interesting feature of this diffusivity enhancement is its weak dependence on the size of the transported particles, as opposed to classical diffusion for which the diffusivity is inversely proportional to the particle size ( Figure  6c2, eq. 7 ).
In vivo front tracking experiments illustrate this flowinduced accelerated transport of 20 nm nanoparticles in paralyzed and active zebrafish embryos ( Figure 6d ). While both active and paralyzed embryos exhibit a diffusivelike transport, the effective diffusivity is higher in the case of contracting embryos. It confirms that muscle contractions can increase the longrange transport efficiency, which act similarly to a constant increase of the local flow velocity ( Figure 6c1 ). Interestingly, our model also predicts that transport is enhanced in both rostrocaudal and caudorostral directions. We indeed demonstrate that beads injected locally in the caudal part of the central canal are transported rostrally towards the brain ventricles ( Figure 6e1 ). Similarly to rostraltocaudal transport, the reverse transport is sped up for contracting embryos ( Figure 6e2 ).

The connection from brain ventricles to central canal in the spinal cord is gated by a funnelshaped channel and is critical for embryonic growth
Now that the ability of the CSF flow to transport particles across the whole CSF zebrafish body has been demonstrated, its consequences for embryonic development have to be discussed. We first summarize the entire CSF circulatory system of zebrafish embryos ( Figure 7a ). The DV is connected to the CC via a small channel making a funnel connection of small constant diameter (0.55 ± 0.12 times the CC diameter and 233 ± 24 μm long, Figure  7b1, b2, b3 , 11 embryos). At the caudal end of the CC, a secondary channel thinner than the CC closes the loop by connecting the CC to the rhombencephalic ventricle (RV) ( Figure 7a ). In 24 hpf to 30 hpf embryos, primary cilia are found in the brain ventricles but they are mostly nonmotile. Therefore brain CSF flow is mostly brownian with a pulsatile component at the heart beat frequency, which is hindered if the heart beat is suppressed ( Supplementary Movie S7 ), indicating that brain ventricles marginally participate to the global circulation of CSF between 24 and 30 hpf. On the contrary, a strong unidirectional flow towards the caudal region is driven through the funnel ( Supplementary Movie S8 ), with a few smaller channels connecting the funnel to the bottom of the RV. Through the funnel, this fast unidirectional flow is characterized by parabolic Poiseuillelike velocity profile of maximal speed 15 μm/s ( Figure 7c1, 7c2 ). We observe motile cilia inside the funnel (data not shown), but the constrained geometry prevents us to analyze cilia distribution and properties in this channel.
The funnel allows a connection between the DV and the CC during development, that we investigate by performing 2photon mediated photoablations of the connection between the funnel and the CC ( Figure 7d1 ), and control ablations in the yolk extension ( Figure  7d2 ). 30 hours after ablation, half of funnelablated larvae were shorter than siblings undergoing control ablation or without ablation ( Figure 7e1, 7e2, 7f1 ) and they often exhibited an increase in bodyaxis curvature ( Figure 7f2 ). We thus investigate whether the CSF contained extracellular vesicles, using labelfree imaging technique, called full field optical coherence tomography ( Figure 7g and Supplementary Movie S9 ). We identify many vesicles with high lipids or proteins content floating in central canal, in agreement with previous observations 31 . The presence of dense exosomes is confirmed using fluorescence labeling ( Supplementary Movie S10, Methods ), and electron microscopy with immunogold labeling of exosomes ( Figure 7h, Methods ) 30 . Note that both optical techniques show similar bidirectional flow of endogenous extracellular vesicles in the CC ( Supplementary Movies S10, S11 ). Altogether, these experiments indicate that CSF circulation contributes to embryonic growth, and suggest that the transport of extracellular vesicles between brain and spinal cord is a critical process.

Discussion
In this study, we have performed a novel quantification of cerebrospinal fluid (CSF) circulation in zebrafish embryos between 24 and 30 hpf, demonstrating for the first time that a constant bidirectional CSF flow is driven by ventral motile cilia and that body movements create additional transient flows. These two driving forces accelerate the long range transport of all particles larger than 10 nm, including macromolecules and extracellular vesicles present in the CC. This permits fast transport both from brain to spinal cord, and from spinal cord to the brain, apparently critical for embryonic development. In addition, we reported two intriguing phenomena : 1) in mutants with ciliary defects, the CC collapses a phenomenon that had been missed in previous characterizations; 2) the density of extracellular lipidic vesicles in CC, possibly containing exosomes, is particularly high at the embryonic stage, as already reported in mammals using fixed tissue or CSF punctures 27,29 . We showed using 2 live imaging techniques that extracellular lipidic vesicles are transported in both directions along the rostrocaudal axis.
Little is known about CSF embryonic flow as its investigation represents a daunting imaging challenge. Indeed, CSF is a clear fluid with little endogenous contrast, and the CC is a small channel difficult to access inside the spinal cord. In healthy young humans and mammals where the CC is open, the canal is often too small to be imaged with phasecontrast MRI 51 , in contrast to brain ventricles or subarachnoid spaces. Here, we take advantage of the transparency of zebrafish embryos to characterize CSF bidirectional flow as well as cilia dynamics in the CC confined geometry. A simple hydrodynamic model accounts for our observations by modeling cilia contribution as a homogeneous volume force . The latter fits experimental data for an amplitude of 4000 N/m 3 , corresponding to a force f v of 10 pN generated by a single cilium occupying a volume of characteristic size 10 μm, in good agreement with the values reported in the literature 52 . This modeling of the force is f v reminiscent to the study of Siyahhan et al. 53 in ventricles of the human brain, where a linearly varying force of comparable magnitude is used instead . We report that including a linear force profile does not significantly alter the predicted velocity profile ( Supplementary Figure  S2b ), which shows the robustness of these averaged models. We acknowledge here that we cannot model the complexity of many interacting cilia beating at many different frequencies, and that only the average flow can be inferred by our model. Interestingly, our model can efficiently predict flow profiles for many other aspect ratios between the cilia extension and the channel diameter and heterogenous ciliary beating frequency ( Supplementary Figure  S2c ). Notably, it convincingly predicts the flow previously reported by Shields et al. in the context of artificial magnetic cilia in microchannels 8 , where this aspect ratio is low. Our model should also be relevant to account for ciliainduced flows in brain ventricles or in their connecting channels where the confinement is moderate.
Our model also predicts the establishment of an overpressure in the CC, increasing linearly towards the caudal direction ( Supplementary Figure S2a ). This inner pressure is dramatically decreased as the density of motile cilia is lower. These findings may be directly related to the collapse of the CC reported for mutants with defective cilia ( Figure 4c ). This suggests that, among other roles related to transport of chemical cues or mechanosensing, the ciliainduced flow ensures the integrity of the CC, necessary for a proper body axis formation. Interestingly, similar transitions from cylindrical to collapsed geometries due to inner pressure decrease were already reported in small veins 1 . It may also account for the apparent fragility of the CC that easily collapses during fixation (Wyart lab, unpublished observations). We report that mutants with defective cilia were unable to maintain an open CC. Elipsa and kurly embryos with partial loss of cilia motility show intermediate phenotypes, with some regions open and other collapsed, and with open regions showing bidirectional flow with a higher density of recirculation zones. We believe this intermediate phenotype is predicted by our "scarce cilia" model ( Figure 3c ), since a few cilia remain motile in kurly embryos 54 and in elipsa between 24 hpf and 30 hpf 39 . The model also accounts for the reduced longrange transport previously observed in kif3b mutants, where cilia density is reduced 36 . In contrast, foxj1a mutant embryos showing ependymal cells with defective cilia 55 ( Supplementary Figure S3 ) have a completely collapsed CC. This observation strongly suggests that the motile cilia generating the flow and maintaining the CC geometry originate from ependymal cells only, while remaining cilia probably belonging to CSF contacting neurons 56 are too short and not dense enough to generate a stable CC ( Supplementary Figure S3 ).
We provided the first description of the narrow channel linking the CC with brain ventricles with a funnelshaped entrance from the diencephalic ventricle. We found another channel connecting the caudal part of the CC to the rhombencephalic ventricle, forming a closed loop ( Figure 7a ). The pressure distribution in the CC induced by the motile cilia could thus also be responsible for a flow in these narrow channels even if they would not contain any motile cilia. The presence of motile cilia, observed in the funnel, may accelerate the flow even further. The cilia in the CC have consequently the potential to drive a large scale pressuredriven flow as well, playing the role of a pump for the whole CSF circulatory system. Unfortunately, flow measurements in these narrow channels are extremely challenging and will be the scope of a future study. Based on our ablation experiments we propose here a critical role for the passage of CSF from the brain to the spinal cord through the funnel, during embryonic development. This correlates well with previous reports of body axis straightening and spine morphogenesis defects arising when the CSF flow in ventricles is altered at later stages 21,36,37 . Another interesting followup, today limited by our imaging technology, is to investigate the role of body movements in freely behaving embryos, knowing that they normally twitch every few seconds.
In this study, we neglected the presence of the Reissner's fiber (RF) in the central canal, although it was shown to play a critical role in body axis straightening 39 . As it does not compartmentalize the fluid in 3D, the presence of the RF does not qualitatively modify the nature of the flow and its bidirectionality, as confirmed by velocity profiles measured in scospondin mutants deprived of RF 39 . However, we expect an increase in the viscous stress at the vicinity of the fiber thereby likely changing the mechanical force applied by CSF flow on cells contacting the CC. As CSFcontacting neurons are known to be mechanosensitive 31,57 , a change in activity of these neurons in the absence of the RF may also occur.
Finally, our study focused on embryogenesis as this is a critical period when CSF shows its highest concentration in secreted extracellular lipidic vesicles 27,29 , and when blood flow along the body and ciliary beating in the ventricles are not yet implemented. As the structure of the central canal is well conserved among vertebrates 58 but very challenging to access with optical imaging in most species, we expect that our findings in zebrafish will guide future investigations in mammals. The thorough study of CSF flow and role in larval, juvenile and adult stages will be the object of future developments, as we expect the CSF flow to be different at later stages since both structure and cilia dynamics in the CC change with time. Altogether, we believe our work to be a major first step into the understanding of the role of CSF transport in central canal to ensure a longrange communication between cells throughout the body that may control development, inflammatory response, and modulation of locomotor and postural functions carried by spinal circuits.  , we filtered any region whose orientation is close to 0 or 90°± 180°(excluding regions whose cosine or sine of orientation are below 0.1). Note that the intensity threshold is not so strong, but that most of the particle filtering is made through these 4 last conditions.

Materials and Methods
A few regions corresponding to individual particle rostrocaudal trajectory are thus obtained.
The tangent of these lines orientation is calculated, and once corrected by pixel size and frame time, the individual particle velocity is obtained. All particles from one kymograph were aggregated to the average velocity and standard error to infer the local flow (which we are interested in rather than individual particles speeds). Note that we refer to the number of events, or number of lines tracked, rather than the number of particles, as this algorithm cannot guarantee that one particle is not tracked several times if its dorsoventral position changes over time.  Comparison between methods to measure local flow (PTV/PIV/and kymograph approach).
Besides the PTV and kymograph approaches, we also tried to use a particle imaging velocimetry (PIV) algorithm, which is generally used to measure local flow 61,62 . This section aims at comparing all these techniques qualitatively with respect to our imaging conditions.
Before describing why we mainly used the automatic kymograph approach, we have to recall that we cannot inject the particles directly in the central canal, and realized that exogenous particles above 100 nm in diameter stayed confined inside the brain ventricles. In turn, we could only image particles of diameter below 100 nm, and mainly used 20 nm particles with fairly low signal, especially at high imaging frequencies. Besides, the Brownian motion of 20 nm beads is important, and can generate instantaneous displacements of the same order of magnitude than the CSF flow itself. With that in mind, PIV should be efficient to measure the local flow even in low SNR conditions, as it uses a cross correlation between successive frames. Nevertheless, because of the importance of Brownian motion, the crosscorrelation fails at producing a peak of high SNR. Also, the geometry of the canal itself, as well as the bidirectional nature of CSF flow, makes that the use of standard large squared regionofinterest on which the cross correlation is typically calculated in PIV is not well adapted here. As particles can move in opposite directions in a 6x6 μm² (32x32 pixels), the cross correlation results in several small competing peaks. Asymmetric region of interests should be more adapted to take a full advantage of the constant flow along the rostrocaudal direction, but performing a PIV with asymmetric regions is not standard and not trivial. On the other hand, PTV tracks every individual particle and should produce the trajectory of each particle versus time. One can expect that a given trajectory would then have a Brownian component, plus a directed component, each of them being separable in principle.
Unfortunately here, PTV is not very robust to low SNR, and not very efficient when particle density is quite high. The most critical part is to efficiently detect the particles without error. A first solution is to have a high threshold that selects only the brightest particles, as what was done in Figure 1. Nevertheless, only a few particles could be tracked for a given dataset, and the profile could not be inferred from only a few particle positions. If the threshold is lowered, it increases the calculation time, and the probability of false detection. Since the PTV algorithm tend to minimize the distance traveled by particles, as the number of false detection increases, the false associations increases even more leading to an important underestimation of particle speed. In contrast, the automatic kymograph seemed more robust to us in this context because the flow is mostly constant in the rostrocaudal direction, which allows us to perform efficient averaging in this direction. Brownian motion is average out by averaging on 35 consecutive dorsoventral positions, so that the particles produce traces long enough to be identified. The main strength of the kymograph approach is that many particles can be detected at every position, so that even if several mistakes are made, the standard error giving the error to find the correct average speed, becomes small. The longer the total acquisition time, the more particles can be averaged, and the more precise the technique is, which is not true for both PTV and PIV in which all images are independent.
Nonetheless, we must admit that most of our efforts have been focused on the automatic kymograph approach, and that probably PTV and PIV algorithms may be adapted to our configuration. We could not find any available PTV and PIV algorithm that could work out directly with our data, and theories behind both techniques become more complicated in non optimal configurations. In contrast, the kymograph just requires to fit lines in an image, for which many solutions are available and just required small adjustments.  3i and Leica, and one 2photon scanning microscope using the scanning artifact, as previously described by 7 , in order to validate our results) . However, all the results displayed in this report were obtained with the same Leica commercial microscope to ensure data homogeneity. We can only postulate that the difference found in beating frequency reflects a bias for long cilia with possible slower kinematics when using differential interference contrast (DIC) to identify cilia.
Quantification of parameters related to cilia structure and dynamics . Previously described First, the central frequency within each region is calculated. Then, each region is associated to an ellipse, whose long axis is assumed to be cilia length. After geometrical calculation, the ellipse orientation gives the angle between the cilia and the dorsoventral axis. The values corresponding to the negative orientations may principally originate from the few beating dorsal cilia, also oriented towards the caudal end, but from dorsal to ventral (leading to a negative angle). To account for this, we calculated the median of the absolute value of the cilia orientation in the results section. Note that it gives the same value than the median of only the positive tilted cilia (65.1°versus 65.0°for the median of absolute values). Finally, we calculate a parameter called beating height, as the cilia length times the cosine of orientation (which therefore does not depend on the sign of the orientation), which gives an estimate of which portion of central canal is occupied by beating cilia. Although some reports 43 used the power spectral density to quantify the main beating frequency, no difference was found here between the main frequency found with the Fourier transform versus the power spectral density calculation.

Declaration of Interest
The authors declare no competing interests.  Figure  2B1. Note the diversity of beating patterns of neighboring cilia. Not only cilia beat at many frequencies, but some of them also show unrepeatable patterns. Horizontal scale bar is 20 µm, and vertical scale bar is 0.5 second. B)Corresponding power spectral density map showing the position of the different maximal frequencies for each cilium. Note that some cilia show a power spectrum with several local maxima, accounting for the asynchronicity of the beating patterns observed. In order to obtain the frequency map as in Figure 2B2, we extract only the global maximum. The overpressure represents the excess of pressure with respect to the inlet of the CC (at the junction with the brain ventricle). In the case of a homogeneous cilia distribution (red curve), modelled as a constant force f v = µf /h, the pressure increases linearly in the caudal direction. In the case of a loose distribution of cilia, represented by an alternation of active regions with the aforementioned force f v and passive regions without any bulk force, the pressure follows a piecewise linear increase, whose frequency is controlled by the distance between active and passive regions : w = d for the blue curve and w = d/2 for the green curve. The cilia activity is thus directly correlated to the excess of pressure in the central canal, as dramatically put in evidence by the magenta flat curve in the absence of any cilia activity. This excess of pressure in the central canal induced by the cilia activity could help maintaining the lumen of the channel open, explaining why the central canal is found to be collapsed in cilia mutants ( Figure  4C).B) In the main manuscript, we modelled the contribution of the cilia as a constant bulk force f v = µf /h in the ventral region occupied by the cilia. In a prior publication, Siyahhan et al. (2014) modelled it as a force linearly increasing with the distance to the wall (also located in the ciliary region). We show in this figure that this choice (blue triangles) does not substantially alter the velocity profile predicted by our model (red squares), which simply turns into a 3 rd degree polynomial in the ventral half of the channel. The y-position for vanishing local velocity is slightly shifted towards the dorsal wall. C) Our model still holds for variations of the relative thickness of the ciliary region and the dorsal region. Following the derivation detailed in the manuscript, the velocity profiles can be obtained for arbitrary values of h, at fixed t = 10 µm. The green curve corresponds to the zebrafish CC features, and the red curves corresponds to very thin cilia with respect to the diameter of the tube. This situation is representative to that encountered by Shields et al. (2010) with artificial magnetic cilia. They decided to model the effect of cilia as a moving wall entraining fluid. While this moving wall modelling correctly accounted for their observations in these high length ratio conditions, we point out that it would not successfully predict the flow in more common length ratios as it fully neglects the flow in the ciliary region. Conversely, our model accounts for the flow in the whole channel, which renders it very versatile to different geometries. D) -In Figure 3C, we have shown by numerical simulations that the alternation of active and passive regions can give birth to vorticity in the flow profile. Here, we show that recirculation regions can also be obtained from the succession of neighboring cilia with different beating frequencies. And the larger the difference in frequency, the larger the vortices. The black lines correspond to the trajectory of the flow, and the colormap in the channel corresponds to the rostro-caudal velocity. The colored rectangles under the ventral wall indicate the frequency for each spot in the simulation for the bulk force f v = µf /h. This choice is consistent with the in vivo measurements detailed in the Figure  2B of the main manuscript.   The gRNA sequence is highlighted in green. B) The Foxj1a protein consists of 458 amino acids and includes a DNA-binding forkhead domain (from amino acid 140 to 230, indicated in red). The Foxj1a nw2 mutant allele results in a truncated protein due to a frame shift and early stop codon and lacks the forkhead domain. C) -Foxj1a nw2/nw2 homozygous mutant larvae display a curved body axis as seen on images of a control sibling (heterozygous or WT) versus a homozygous mutant larva taken with a stereomicroscope at 4 dpf. D) Z-Projection stack of a lateral optical sections (Z-step = 1 µ over 3 sections) of the spinal cord immunostained with DAPI and against acetylated tubulin in a 30 hpf control sibling embryo (D1) and a foxj1a nw2/nw2 curled-down embryo (D2). E) Z-Projection stack of lateral optical sections (Z-step = 1 µm over 3 sections) of the spinal cord stained with DAPI and immunostaining against glutamylated tubulin in a 30 hpf control sibling (E1) and a foxj1a nw2/nw2 curled-down embryo (E2). Scale bars are 15 µm.

Supplementary Information
Both immunostaining for acetylated tubulin and glutamylated tubulin show a decrease in density, and in length of cilia in central canal in foxj1nw2/nw2 homozygous mutant embryos. Glutamylation of cilia is usually associated with active motility, which is coherent here with the position (ventral) and polarity (caudally-oriented) of glutamylated cilia in CC. The phenotype observed in the immunostaining against glutamylated tubulin usually labeling motile cilia is more drastic and suggests that most motile cilia from ependymal cells do not form in the foxj1a nw2/nw2 homogygous mutant embryos.

IV. 2D MODEL FOR THE DIFFUSIVE-CONVECTIVE TRANSPORT OF SOLUTES IN THE CENTRAL CANAL WITH A BIDIRECTIONAL FLOW-ASSOCIATED WITH FIGURE 6
We here provide further details concerning the derivation of the effective diffusivity D ef f of particles in the central canal in presence of a bidirectional flow.
In the CC , the motile cilia occupy the ventral half (Fig. 2 of the main manuscript). From the flow model derived earlier (eqs. (3) and (4) of the main manuscript), one can rewrite the velocity profile changing the y-coordinate origin (y = 0 at the center of the channel) and using the maximal velocity V max = 6 µm/s: in the dorsal region (positive values of y), and : v(y) = 16V max y 2d + y 2 d 2 (2) in the ventral region (negative values of y). The full transport equation of particles of local concentration c writes: where v(y) corresponds to the two velocity profiles written previously (eqs. (1) and (2)), as the local velocity field is uniquely oriented in the rostro-caudal z-direction (see also Fig. 3 of the main manuscript). D is the molecular diffusion coefficient, here taken as D = 10 −11 m 2 /s, corresponding to particles of diameter 40 nm using the Stokes-Einstein law (eq. (7) in the main manuscript). Several terms of equation (3) can be neglected considering their order of magnitude. With L the length of the channel, and d its diameter, the equation (3) may be written in dimensional analysis for times of order T L 2 /D = τ dif f (characteristic spanwise diffusion time) : i.e.
, or : Both V max L/D and L 2 /d 2 terms are much larger than unity, due to the high aspect ratio of the channel L/d. This allows to simplify equation (3) in the two dorsal (y > 0) and ventral (y < 0) regions, respectively, as follows: and 16V max y 2d + y 2 d 2 ∂ z c = D∂ 2 y c.
Adapting the derivation ofBruus (2008), that was done for a Poiseuille flow in a cylinder, we first assume that the axial derivative of the concentration ∂ z c is independent on y. This assumption will be verified a posteriori and confronted to numerical simulations. Consequently, eqs. (7) and (8) are ordinary differential equations for c(y), that can be solved. A first integration between y and y = d/2 for equation (7) for y > 0 leads to 16 V max d 2 Similarly, integrating equation (8) between y = −d/2 and y for y < 0 leads to By antisymmetry, ∂ y c | y=d/2 = −∂ y c | y=−d/2 . The equality of the concentration y-derivative in y = 0 then leads to ∂ y c | y=d/2 = −∂ y c | y=−d/2 = 0, which simplifies equations (9) and (10). A second integration between y = 0 and y of eq. (9) for y > 0 leads to Similarly, a second integration between y = 0 and y of the equation (9) for y < 0 leads to 4V max d 2 3D We now aim at expressing the concentration profile with respect to the average concentrationc over a cross sectionc which, using eqs. (11) and (12) dy + c(0) (14) and finally This theoretical prediction is confirmed by numerical simulations solving the complete transport equation without the assumptions made for the theoretical model (Fig. S3A), for different times, positions in the canal and flow velocities. Going back to our theoretical model, one can now derive the condition for having a y-independent axial gradient of concentration (condition required for recovering a diffusivelike concentration profile) by differentiating the equations (11) and (12) with respect to z. The latter condition is consequently fulfilled only if ∂ z c = ∂ z c(0), which requires the third term to be negligible. In terms of order of magnitude, this implies that where L is the length of the central canal. This condition can be rewritten as a condition on the Péclet number Pé: This condition on the Péclet number is always fulfilled in the CC of zebrafishes for nanobeads or lipidic vesicles (D ∼ 10 −11 m 2 /s), where Pé does not exceed 10.
We now aim at deriving the effective diffusivity of particles in the CC. We calculate the average current densityJ (z) through the cross-section, using the velocity profiles (eqs. (1) and (2)) and the concentration profiles (eqs. (11) and (12)): ρc(z, y)v(z, y)dy.
The zero average flux over the cross section d/2 −d/2 v(y)dy leads to the vanishing of the constant componentc of eqs. (11) and (12). The previous expression may then be simplified :J i.
which may be approximated toJ The effective diffusivity D ef f may then be directly from the effective Fick's lawJ (z) = D ef f ρ∂ zc : A more detailed calculation in a cylindrical channel was done by . It generalizes this result at intermediary Péclet numbers Pé 1. Adapted to our geometry and flow, it yields .
For particles of radius 20 nm (corresponding to the size of injected nanobeads), we have D 5 · 10 −12 m 2 /s, which leads for V max = 6 µm/s to an effective diffusivity D ef f 2.5D. Similarly, the effective diffusivity of the same particles for a faster flow of amplitude V max = 20 µm/s is D ef f 18D. These two theoretical predictions are in good agreement with the numerical simulations (Fig. 6C of the main manuscript).
Relevance of the initial hypothesis. One strong assumption of the model that enabled to readily integrate eqs. (7) and (8) is that the axial derivative of the concentration ∂ z c is only weakly varying over a cross section, such that we consider it constant. The Fig. S4B plots the relative variation of the normalized concentration gradient ∂ z c/∂ z c(0) as a function of y. It shows that this quantity exhibit minor variations of about 10% over a cross section, demonstrating the relevance of the simplifying assumption.  varies of about 10% around its mean value over a cross section, which justifies the assumptions done in the model (from eqs. (7) and (8) to eqs. (9) and (10)) to neglect its variations in comparison with the strong variations of the other integrand y 3 /d 3 + y 4 /d 4 − y/4d .