Z Z E L L I

Sediment dynamics. Part 2. Dune formation in pipe flow

 accounts reasonably well for the experimental observations.

Introduction

A very common feature that arises when bed constituted of particles are submitted to shearing flows is the formation of ripples, i.e. small waves on the bed surface having wavelengths scaling with the particle size, or of dunes, i.e. larger mounds or ridges having wavelengths comparable to or larger than the fluid height.

Since the seminal work of Kennedy (1963), the prevalent mechanism for dune or ripple formation is the fluid inertia or more precisely the phase-lag between the bottom shear stress and the bed waviness generated by the fluid inertia. In that case, the shear stress, the maxima of which are slightly shifted upstream of the crests, drags the particles from the troughs up to the crests. This destabilising mechanism seems to be robust enough to apply to any steady flow, either turbulent, see e.g. Engelund (1970), [START_REF] Richards | The formation of ripples and dunes on an erodible bed[END_REF], [START_REF] Sumer | On the formation of ripples on an erodible bed[END_REF], [START_REF] Colombini | Revisiting the linear theory of sand dune formation[END_REF] and [START_REF] Claudin | A scaling law for aeolian dunes on Mars, Venus, Earth, and for subaqueous ripples[END_REF], or viscous, see e.g. [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF], and [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF].

A few stabilising mechanisms have been proposed, among which a well identified stabilising effect due to the inclination of the bed, see e.g. [START_REF] Fredsøe | On the development of dunes in erodible channels[END_REF], [START_REF] Richards | The formation of ripples and dunes on an erodible bed[END_REF], [START_REF] Sumer | On the formation of ripples on an erodible bed[END_REF], [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF], [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF] and [START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF]. For nonzero slope, the gravity force parallel to the bed surface favours the downhill motion of the particles and conversely impedes the uphill motion. Another stabilising mechanism related to particle inertia was also suggested to arise from the phase-lag between the bottom shear stress and the particles flow rate. This effect can be expressed in term of an inertial length, see e.g. [START_REF] Andreotti | Selection of dune shapes and velocities. Part 2: a two-dimensional modelling[END_REF], or a deposition length coming from a stabilising crest-erosion process, see e.g. [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF] and [START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF].

However, a complete description of the bed instability is still lacking as the coupling between the granular media and the fluid is poorly understood. Usually, one first calculates the fluid flow as if the wavy bottom were fixed by considering the superposition of a base flow on the flat bed and a perturbation induced by the wavy bottom. Then, one needs to relate the particle flow rate to the calculated shear stress at the top of the bed. Table 1. Stability analysis prediction for the instability threshold. The particle Reynolds number is Re p (with d the particle size, u * the shear velocity, and ν the kinematic viscosity of the fluid). The Shields number, θ, is constructed as the ratio of the shear stress at the top of the bed to the apparent weight of a single particle. The critical Shields number for incipient particle motion is θ c 0 , the friction coefficient µ, and the fluid height is h f . The Galileo number Ga = d 3 (ρ pρ f )g/ν (with ρ p and ρ f the particle and fluid density, respectively) is the Reynolds number based on the Stokes settling velocity of the particles.

Particle mass conservation equation is finally solved to provide the linear growth rate of the instability. Several algebraic law relating the particle flow rate to the bottom shear stress have been proposed in the literature, see tables 1 and 2 of [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF], leading to different control parameters for the instability threshold, such as the Reynolds, Shields, or Galileo numbers, see table 1 where some predictions of linear stability analyses, which are of particular interest for the present study, are presented.

There is not yet a complete experimental proof that this type of modelling captures the essence of the instability. The stability analyses of [START_REF] Sumer | On the formation of ripples on an erodible bed[END_REF] and [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] having the particle transport described by an algebraic law, a Bagnold-type law and a cubic law respectively, have found preferred initial wavelengths to be an order of magnitude smaller than the observed dune lengths. A recent stability analysis by [START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF] seems to improve the predictions by including a phase-lag which erodes the peaks, i.e. the additional stabilisation mentioned above.

The present contribution aims at investigating dune formation in a confined, wellcontrolled flow, namely a flow in a closed pipe. In § 2, we present an experimental study of the evolution of a particle bed in a pipe flow. Different dune patterns are observed as the flow rate is increased from the laminar to the turbulent regimes. We then focus on the threshold for destabilisation of the flat bed leading to dune formation. In § 3, we adapt the stability analysis of [START_REF] Charru | Phase diagram' of interfacial instabilities in a two-layer Couette flow and mechanism of the long-wave instability[END_REF] to a Poiseuille flow and choose to use simply the algebraic law relating the particle flux to the cube of the Shields number found by [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF]. The control parameter of the instability is the Reynolds number and the predicted wavelength at onset scales with the fluid thickness. In § 4, a phase diagram for the different dune patterns observed is proposed and comparison with theoretical predictions provided.

Experimental observations

Experimental set-up

The experimental test section was a horizontal glass tube of inner diameter D = 3 cm and length L = 1.8 m. The measurement zone was located at ≈ 50 cm from the entrance. This length corresponded to the entry length necessary for the laminar flow to develop fully inside the tube at Re pipe ≈ 250, where Re pipe = 4Q pipe ρ f /πDη is the tube Reynolds number with Q pipe the flow rate, ρ f the fluid density, and η the viscosity. Experiments were performed in the range 10 -1 Re pipe 10 4 . For Re pipe 250, the laminar flow was then fully developed in all the measurement zone while it was not for 250 Re pipe 2500. The transition toward turbulence occurred for Re pipe ≈ 2500.

Three different batches of spheres and five different mixtures of distilled water and UCON oil 75H-90000 were used in the experiments as indicated in tables 2 and 3. First, the tube was filled with fluid and the particles were carefully introduced to form an uniform flat bed. Secondly, in order to avoid flow perturbations from a pump, the flow was driven by gravity. The fluid was delivered to the tube by continuous overflow from an overhead tank of variable height. At the outlet of the tube, the particles were captured by a mesh while the fluid was run into a reservoir. A pump carried the fluid back to the overflowing tank. The temperature of the fluid was maintained at 35 ± 1 • C by using a thermostated bath as a fluid reservoir in the fluid circulating loop. To measure Q pipe , we collected a given volume of fluid at the outlet of the tube in a given time. It is worth mentioning that, even though the flow is controlled by gravity, the pressure losses in the flow loop result in an imposed flow-rate in the test section despite the changes in the tube section as particles are lost.

Classification of dune patterns

Five different regimes can be observed as summarised in figure 1 for two different batches of particles (batches A and C) and using fluids of varying viscosity to explore the full Re pipe range. For convenience, we have chosen to represent their domain of validity as a function of Re pipe . Below the critical Shields number for incipient particle motion, θ c 0 ≈ 0.12 (the 0 subscript indicates that this threshold corresponds to a flat bed), the bed shape does not evolve, see [START_REF] Charru | Erosion and deposition of particles on a bed sheared by a viscous flow[END_REF], Loiseleux, Gondret, Rabaut & Doppler (2005), and [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. Note that this first threshold is determined by the Shields number and not the Reynolds number. Above this first threshold, we observe a second regime where the bed shape evolves but without dune formation. In this regime, a few layer of particles are set in motion at the surface of the bed by the fluid flow. The bed shape evolves, becomes slightly tilted and the bed height decreases with increasing time as the test section is not fed in with particles. When the experiment is run for a long enough time, the bed shape eventually freezes, exhibiting a flat surface.

Increasing again Re pipe , three regimes exhibiting different dune patterns can be found. For the first dune regime, which only exists in laminar flow, the dunes denoted 'small dunes' present small amplitudes. The second dune regime is observed either in laminar or turbulent flow. In this regime, the dunes are characterised by the existence of vortices located at their front as can be seen in online movie 1. The vortices erode the particle bed (see figure 1 top view of the vortex dunes). This thus leads to dunes of very large amplitude denoted 'vortex dunes'. Finally, when the flow is turbulent (and only in that case), we observe a new dune pattern that we called 'sinuous dune'. The bottom photograph of figure 1 shows a top view of a 'sinuous dune'. The particle bed is eroded asymmetrically leading to the formation of a pattern having a double periodicity. This may be explained by a destabilisation of the initial vortices observed in the 'vortex dunes'.

Evolution of dune amplitude, wavelength, and phase velocity

To quantify more precisely the time-evolution of the dunes, we recorded the evolution of the bed height as a function of time, for a combination of fluid and particles, a given flow rate Q pipe , and an initial height of the bed, h start p . Using the same experimental technique as in [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF] where further details can be found, the bed height was measured by imaging the upper layer of the bed illuminated by a laser sheet aligned with the tube length in its middle. The measurement zone spanned over 45 or 75 cm and started at ≈ 50 cm from the entrance of the tube. The use of such a large measurement zone provided a precise study of the wavelength evolution over a long time but leaded to a poor resolution in dune amplitude as the accuracy in height was ≈ 2 -3 mm.

Here, we focus on the two first dune regimes, i.e. the regimes of small and vortex dunes. The spatio-temporal evolution of the bed height for three experiments are plotted in figure 2, in the small dune regime (a) and in the vortex dune regime (b) and (c), with bed heights coded in grey scale. In both regimes, we observe the formation of initial dunes in the measurement zone. These dunes have a relatively uniform wavelength and move in the flow direction inside the tube. However, the initial dunes preceding the vortex dunes differ from the small dunes as will be evidenced later.

In the case of vortex dunes, the flow perturbation at the entrance of the tube generates a first vortex dune whose vortices create a second dune downstream and so on. The vortex perturbation by propagating downstream creates dunes downstream over the entire length of the tube which progressively take over the initial dunes, see the top plots of figure 2 (b) and (c). At the same time, the dune pattern moves downstream as the dunes themselves are moving in the flow direction. With increasing time, the vortices between two adjacent dunes erode the particle bed and the dunes are separated by regions emptied of particles. Even though no dune coalescence is observed, the pattern wavelength increases with time whereas the dune velocity decreases. The dune motion eventually stops as the tube is not fed in with particles and the dunes are completely separated, see the bottom plots of figure 2 (b) and (c). We can note that the length-scale depends on the particles. The time-scale also differs. For batch C in fluid 1, the dunes stop after less than one week, while, for batch A in fluid 1 with a flow rate of the same order of magnitude, they move during approximately two weeks before stopping.

In the case of the small dunes, there is also a propagation of dunes created by the entrance perturbation. But, as these dunes have no vortex, the velocities of dunes due to the entrance perturbation and of the initial dunes are similar and thus these dunes do not overtake the initial dunes created in the measurement zone, see the top plot of figure 2 (a). As time is increased, the wavelength increases without showing any saturation, the dune velocity decreases, and the dune amplitude increases, see the bottom plot of figure 2 (a). The small dunes never stop moving downstream and eventually leave the tube when they reach the tube outlet. When all the dunes have left the tube, the bed becomes flat again. A strong change in the pattern can be observed at the bottom of figure 2 (a) as a ripple instability is observed.

Small dunes and vortex dunes exhibit different behaviours. While the particle bed eventually freezes in both cases, their final states also differs. For vortex dunes, due to presence of vortices, the final state corresponds to immobile dunes separated by regions emptied of particles. For small dunes, as there are no vortices to dig out the particle bed, the small dunes move downstream until they leave the pipe and the remaining state is a flat bed. These two behaviours are shown in online movies 2 and 3. Beside, their length-scales and time-scales are dissimilar even at the initial stage.

To obtain a quantitative comparison, we have measured the amplitude, wavelength, and phase velocity of the dunes in both regimes. The bed height measurements described above are further analysed to determine the local minima and maxima of the height profile. For each acquisition time, the wavelength is given by the average of the distance between two maxima over the measurement zone. In a similar way, the amplitude is given by the averaged difference in height between the minima and maxima. The velocity is given by averaging the displacements of maxima over typically 5 acquisition times.

Figure 3 shows the evolution of amplitude, wavelength, and phase velocity with time of dunes appearing in the measurement zone in the small ( ) and vortex (•) dune regimes. The quantities have been made dimensionless by using D as a length scale and the mean velocity of the flow U = 4Q pipe /πD 2 as a velocity scale. The time scale of dune formation differs by a factor ten between the two dune regimes. The amplitude and wavelength are similar in the first instant. However, the phase velocity of the initial dunes is about ten times larger in the vortex dune regime than in the small dune regime as shown in figure 3 (b). The behaviour also differs as it increases initially in the vortex dune regime while it presents a slow decrease in the small dune regime. This clearly shows that the dynamics of the initial dunes which precede the fully developed vortex dunes differ from those of the small dunes. At later instant, for a dimensionless time ≈ 900 in figure 3, there is a dramatic increase of both amplitude and wavelength as well as a maximum in the phase velocity in the vortex dune regime. This corresponds to the take-over of the fully developed vortex dune having vortices which intensify the dynamics. For further times, the amplitude and wavelength of the vortex dunes saturate as the dunes progressively become detached entities. Their speed decreases and the dunes eventually stop moving. Small dunes have an amplitude about five times smaller than vortex dunes at their later stage. Their wavelength shows a slow and continuous increase.

A simple linear stability analysis

Poiseuille flow on a wavy bottom

To determine the fluid flow over a wavy bottom in a two-dimensional channel, we follow the approach of [START_REF] Charru | Phase diagram' of interfacial instabilities in a two-layer Couette flow and mechanism of the long-wave instability[END_REF] initially undertaken for a Couette flow that we adapt for a Poiseuille flow. We consider a fluid layer lying between a flat upper wall and a wavy bottom which is assumed to be perturbed sinusoidally as ξ = ξ 1 cos(kx) with wave-number k and amplitude ξ 1 as sketched in figure 4. Following the previous studies of [START_REF] Charru | Phase diagram' of interfacial instabilities in a two-layer Couette flow and mechanism of the long-wave instability[END_REF] and [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF], we assume that the time scale of the fluid flow is much shorter than the time scale of the bed evolution. The fluid flow can then be calculated as if the wavy bottom were fixed, by considering the superposition of a base flow, u 0 along the x direction, over a flat bed and a disturbance, u 1 and v 1 along the x and y directions respectively, induced by the wavy bottom.

We decide to make all the values dimensionless by scaling the length by the channel thickness h f , the velocity by q f /h f where q f is the fluid flow rate, and the pressure by a viscous pressure η f q f /h 2 f . Therefore the time is scaled by h 2 f /q f . We note Re 2D = q f /ν the Reynolds number of the channel. The dimensionless velocity and pressure are decomposed into a base Poiseuille flow, ū0 = 6ȳ(1ȳ) and ∂ p0 /∂ ȳ = -ρ f gh 3 f /η f q f , and a disturbance [ū 1 (ȳ), v1 (ȳ), p1 (ȳ)]e i k x. Here, the upper bar indicates dimensionless values, the 0 subscript the base state, and the 1 subscript the perturbation.

Substituting this flow into the dimensionless linearized mass and momentum conservation equations, we obtain a set of linear ordinary differential equations

i kū 1 + ∂v 1 ∂ ȳ = 0, Re 2D i6 kū 1 ȳ(1 -ȳ) + 6v 1 (1 -2ȳ) = -i k p1 -k2 ū1 + ∂ 2 ū1 ∂ ȳ2 , i6Re 2D kv 1 ȳ(1 -ȳ) = - ∂ p1 ∂ ȳ -k2 v1 -i k ∂ ū1 ∂ ȳ2 . (3.1)
The kinematic boundary conditions become ū1 (0) = -6 ξ1 and ū1 (1) = v1 (1) = v1 (0) = 0.

(3.2) Equations (3.1) with the boundary conditions (3.2) can be solved numerically using a Chebychev spectral collocation method, see [START_REF] Gottlieb | The mechanics of dunes and antidunes in erodible-bed channels[END_REF]. It is also interesting to find an analytical solution in the shallow viscous regime, k ≪ 1, by performing an asymptotic expansion for the small dimensionless wave-number k. Assuming Re 2D = O(1), the fluid velocity components can be expended in powers of k which gives the shear-rate

dū 1 dȳ = -6 ξ1 (6ȳ -4) + i ξ1 Re 2D 70 504ȳ 5 -1260ȳ 4 + 840ȳ 3 -108ȳ + 12 k + O( k2 ). (3.3)
This exhibits the phase-lag of the bottom shear rate disturbance kRe 2D /140 due to inertia similar to that found by [START_REF] Charru | Phase diagram' of interfacial instabilities in a two-layer Couette flow and mechanism of the long-wave instability[END_REF] for Couette flow.

Dune formation

The dune growth is determined by the conservation equation for the particle flux, q p , which is obtained by integrating the particle mass conservation equation over the bed height

∂ qp ∂ x + φ 0 ∂ ξ ∂ t = 0, (3.4)
where φ 0 is the particle volume fraction inside the bed. We assume that the dynamics of the particle is well accounted by the algebraic law relating the dimensionless particle flux to the Shields number found by [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF]

qp = φ 0 η η e Ga 24Re 2D θ 3 θ c 2 , (3.5)
where Ga = ρ f ∆ρgd 3 /η 2 is the Galileo number, θ = (ηdu/dy)/∆ρgd the Shields number, and η e an effective viscosity of the mixture of the particles and fluid that was found to equate well to the Einstein viscosity η(1 + 5/2φ 0 ). This algebraic law has been shown to be valid for moderate Shields numbers 0.5 θ 1.5, i.e above the threshold for incipient motion (the moving thickness is larger than one particle diameter) but for shearing flow not substantially perturbed by the motion of the granular media, see [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF].

We now introduce the time evolution of the bed surface, ξ = ξ1 e i( kx-ω t) . We decompose the Shields number into a base Shields number θ 0 corresponding to a flat bed and a perturbation θ 1 e i( kx-ω t) . In the same way, we write the critical Shields number as θ c = θ c 0 + θ c 1 e i( kx-ω t) . Linearising equation (3.5), we obtain

qp = η η e φ 0 Ga 24Re 2D θ 3 0 θ c 0 2 (1 + 3 θ 1 e i( kx-ω t) θ 0 -2 θ c 1 e i( kx-ω t) θ c 0 ). (3.6)
The local inclination of the bed surface modifies the critical Shields number, see e.g [START_REF] Fredsøe | On the development of dunes in erodible channels[END_REF], [START_REF] Richards | The formation of ripples and dunes on an erodible bed[END_REF], [START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF][START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF][START_REF] Charru | Selection of the ripple length on a granular bed sheared by a liquid flow[END_REF], which becomes

θ c = θ c 0 (1 + ∂ ξ/∂x µ ) = θ c 0 + θ c 0 µ i k ξ, (3.7)
giving by identification

θ c 1 = θ c 0 µ i k ξ1 , (3.8)
where µ is a friction coefficient which is the tangent of the angle of repose of the grains and mainly depends on the grain geometry.

We suppose that the time evolution of the bed surface is slow enough to relate the Shields number to the shear rate calculated at the top of the fixed wavy bottom found in the preceding section

θ = Re 2D Ga ( d h f ) 2 dū 0 dȳ ( ξ) + dū 1 dȳ ( ξ)e i( kx-ω t) .
(3.9)

The linearized shear rate calculated at the top of the fixed wavy bottom is given by the two equations

dū 0 dȳ ( ξ) = dū 0 dȳ (0) + ξ d 2 ū0 d 2 ȳ (0) + O( ξ2 ), dū 1 dȳ ( ξ) = dū 1 dȳ (0) + O( ξ2 ), (3.10) exhibiting θ 0 = Re 2D Ga ( d h f ) 2 dū 0 dȳ (0), θ 1 = Re 2D Ga ( d h f ) 2 ξ1 d 2 ū0 d 2 ȳ (0) + dū 1 dȳ (0) + O( ξ2 1 ) . (3.11)
Considering the time evolution of the bed surface and introducing the linearized flowrate given by equation (3.6), the conservation equation of the particle flux at leading order in ξ1 becomes

ω ξ1 = k η η e Ga 24Re 2D θ 3 0 θ c 0 2 3 θ 1 θ 0 -2 θ c 1 θ c 0 .
(3.12)

Clearly, the perturbation in critical Shields number due to the local inclination of the bed surface (which is purely imaginary) has always a stabilising effect while the imaginary part of the perturbation in Shields number contains the destabilising effect due to inertia. The frequency ω is directly given by equation (3.12) where θ 1 is obtained using the numerical computation of equations (3.1) with the boundary conditions (3.2) mentioned in the preceding section. The frequency can be separated into an imaginary part and a real part, giving the growth rate ωi and the phase velocity ωr / k. In the shallow viscous regime, k ≪ 1, equation (3.12) can be solved analytically by using the shear rate given by equation (3.3). At leading order in k, this gives

ω = k η η e Ga 24Re 2D θ 3 0 θ c 0 2 (6 + i 3Re 2D 35 k -2i k 1 µ ), (3.13)
which leads to

ωr k = η η e Ga 4Re 2D θ 3 0 θ c 0 2 , (3.14) ωi = k2 η η e Ga 24Re 2D θ 3 0 θ c 0 2 ( 3Re 2D 35 -2 1 µ ).
(3.15)

The instability threshold, corresponding to ω i = 0, occurs at

Re c 2D = 70 3µ . (3.16)
The important finding is that the control parameter for the dune instability is the Reynolds number, Re 2D , and not the Shields number, θ. It is worth mentioning that the threshold for having a monolayer of particles in motion on a flat bed has been found to be θ c 0 = µφ 0 /2 which is proportional to the friction coefficient µ and to the volume fraction of the particle inside the bed φ 0 , see [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF]. The threshold for dune instability involves the Reynolds number and is also related to the friction coefficient. Using the above expression for θ c 0 , it can be expressed as

Re c 2D = 35φ 0 3θ c 0 .
(3.17)

These thresholds, θ c 0 and Re c 2D , differ but are related. Consequently, the onset for the instability may not coincide with the onset for particle motion. In this simple modelling, these thresholds only depend on two physical parameters, the particle volume fraction inside the bed and the friction coefficient.

Figure 5 shows the dimensionless growth-rate, ωi , versus dimensionless wave-number, . The different regimes represented are: 'no motion' (×), 'flat bed in motion' ( ) and ( ) when outside the domain of validity of the model, 'small dunes' (•), 'vortex dunes' ( ), and 'sinuous dunes' (△). The dashed line is the predicted threshold for particle motion, the horizontal solid line is the predicted instability threshold, and the dotted line indicates the domain of validity of the particle flux law. The dashed-dotted line is the prediction for instability threshold of [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] and the vertical solid line is that of [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF].

k, for φ 0 = 0.55 and θ c 0 = 0.12 as found experimentally by [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF] (or equivalently µ = 0.43). Equation (3.15) valid in the shallow viscous regime, k ≪ 1, shows a long wavelength instability with a threshold at Re c 2D = 54.3. The numerical predictions present a good agreement with this asymptotic case for k 0.1. However, the numerical solution indicates that the instability is not a long-wave instability at threshold but presents a finite value k ≈ 1.7 for a slightly different threshold Re c 2D = 49.3, see blow-up in figure 5. This is an interesting finding as it indicates that the wavelength at onset is of the order of the fluid thickness.

Discussion and conclusion

Figure 6 presents the phase diagram of the dune patterns in the plane Re pipe , Ga(h f /d) 2 . We choose this plane to exhibit both the threshold for incipient particle motion controlled by the Shields number and that for dune instability predicted to be controlled by the Reynolds number in the linear stability analysis of § 3. In this plane, the threshold for particle motion is given by the dashed line Re pipe = (2θ c 0 /3βπ)Ga(h f /d) 2 ≈ 0.014 Ga(h f /d) 2 with θ c 0 = 0.12 as found experimentally by [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF] with the same experimental apparatus. The predicted instability threshold is the horizontal solid line Re pipe = 140φ 0 /3βπθ c 0 ≈ 37. Indeed, in the case of a pipe flow, the Reynolds number is Re pipe = 4ρ f Q pipe /πηD and the threshold equation (3.17) obtained for a two dimensional channel needs to be modified by replacing Re 2D by βπRe pipe /4 (or equivalently q f by βQ pipe /D) where β = 1.85 has been found from numerical analysis in the limit 0.2 h f /D 0.8, see [START_REF] Ouriemi | Determination of the critical Shields number for particle erosion in laminar flow[END_REF]. Note that this predicted threshold only depends on a single parameter (µ or θ c 0 ). The dotted line indicates the domain of validity of the algebraic law relating the dimensionless particle flux to the Shields number found by [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF] and thus indicates the domain of validity of the instability threshold prediction of § 3. The three regimes of 'no motion' (×), 'flat bed in motion' ( and when outside the domain of validity of the model), and 'small dunes' (•) are well delineated by these boundaries in the given limit of validity. Clearly, the threshold for incipient particle motion and that for small dune instability are observed to differ as there is a large region of 'flat bed in motion' without any dune formation. The instability threshold is well described by Re pipe as a control parameter and not by θ which would be a line parallel to the dashed line nor by Ga which would be a vertical line in figure 6. The threshold prediction of the simple linear stability analysis of § 3 (horizontal solid line Re pipe ≈ 37) provides a correct boundary for the 'small dune' instability. Furthermore, the regimes of 'vortex dunes' ( ) and 'sinuous dunes' (△) seem separated and their thresholds also well described by Re pipe as a control parameter.

We have also tested in figure 6 the prediction for instability threshold of Charru & Mouilleron-Arnould (2002) (dashed-dotted line) and [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF] (vertical solid line) listed in table 1. Below both of these predicted thresholds, no dune is observed as expected but these boundaries do not properly delimitate the regime of small dune formation. Note that we did not calibrate the adjustable coefficients of these two models but used the original calibration of [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] and [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF] respectively. Different coefficients would produce boundaries having the same slope as the dashed-dotted and vertical solid lines. These would still be unable to describe the small dune threshold in figure 6. The difference between the stability analysis of § 3 and that of [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] lies into the algebraic law relating q p to θ. The analysis of § 3 uses q p ∝ θ 3 /θ c2 while [START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] takes q p ∝ (θθ c ) 3 . Using a power law involving the excess Shields number yields a threshold depending both on the Reynolds and Galileo numbers (see table 1) and having the slope of the dashed-dotted line in figure 6. The analysis of [START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF] differs from the present analysis by the use of a different model for particle transport which is meant to be applied only to a mobile monolayer and which introduces the new stabilising effect of crest erosion. The competition between this new stabilising effect and the destabilising fluid inertia depends on the Galileo number. If the Galileo number is smaller than a critical number (given in table 1 and vertical solid line plotted in figure 6), crest erosion overcomes and the bed is stable. In the opposite case, the bed is unstable above a critical Shields number. This last threshold would correspond to a line with a slope similar to that of the dashed line in figure 6 which clearly does not delimitate the observed instability.

Further quantitative comparisons with the predictions of the stability analysis of § 3 can be attempted for the onset of small dunes. Figure 7 (a) compares the initial dune velocity with equation (3.14). Despite large error bars and some dispersion of the data, the agreement is good. Figure 7 (b) presents the initially observed wave-numbers as well as the most amplified numerical wavenumber as a function of Re pipe . Clearly, a long-wave instability is not observed in the experiments and a finite value ≈ 1.2 h -1

f is obtained at threshold. The experimental wave-numbers seem rather independent of Re pipe with values ≈ 1 -5 h -1 f while the numerical wavenumber is ≈ 1.7 h -1 f at threshold and presents an increase with increasing Re pipe . Nonetheless, same order of magnitude is recovered close to threshold. It should be, however, mentioned that experimental measurements at instability onset are very delicate and thus do not permit further detailed comparisons. Note that, as expected, the stability analysis developed here is unable to account for the onset of the vortex dunes. Equation (3.14) underestimates by a factor 10 4 the experimental velocity of the vortex dunes. Note also that, while in figure 6, the experimental data for small dunes are obtained with particles of batches A to C with fluids 1 to 5, in figure 7 only the data for particles of batch A and fluids 2 to 5 are used.

In conclusion, we have given the experimental phase diagram of the different dune patterns observed when a bed composed of spherical particles is submitted to a shearing flow in a pipe. 'Small dunes' present small amplitudes and only exist in laminar flow. 'Vortex dunes' are characterised by the existence of vortices at their front and are found either in laminar or turbulent flow. 'Sinuous dunes', showing a double periodicity, appear in turbulent flow. While the threshold for incipient motion is determined by the Shields number, that for dune formation seems to be described by the Reynolds number and not by the Shields or Galileo numbers. Moreover, the dune instability is not a longwave instability at threshold but do present a finite vawelength of the order of the fluid thickness.

To predict the small dune formation, we have performed a simple linear stability analysis where inertia in the fluid produces a phase-lag in the shear stress which is destabilising, while the component of gravity down an incline stabilises the perturbations. We first calculated the perturbed fluid flow over a wavy bottom considered as if fixed. Then, we used the particle flux found by [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF] to relate the bed height evolution to the shear stress at the top of the bed through the particle mass conservation. The threshold for dune formation is found to be controlled by the Reynolds number. This threshold prediction with a single adjustable parameter that we have taken to be realistic provides a correct boundary for the onset of 'small dunes'. The predicted wavelength at instability threshold is of the order of the fluid thickness in agreement with the experiments.

This simple stability analysis containing the basic ingredient of the destabilising fluid inertia and stabilising gravity is found sufficient to provide realistic predictions. It is worth pointing out that the two-phase nature of the problem has been only accounted for in the particle conservation equation in the present study and in most of the studied found in the literature. The interesting questions are whether there is any delay in the flux adapting to a change in the shear stress and whether this additional stabilising mechanism as well as particle inertia and feed back production are significant. Undertaking a stability analysis using the complete two-phase modelling developed by [START_REF] Ouriemi | Sediment dynamics. Part 1. Bed-load transport by shearing flows[END_REF] should enable to answer these issues in the future.

FlatFigure 1 .

 1 Figure 1. The different dune patterns.

Figure 2 .Figure 3 .

 23 Figure 2. Spatio-temporal plots for (a) batch A in fluid 2 with Q pipe = 3.40 10 -5 m 3 s -1 and h start p = 22 mm (small dune regime), (b) batch A in fluid 1 with Q pipe = 3.79 10 -5 s -1 and h start p = 21 mm (vortex dune regime), and (c) batch C in fluid 1 with Q pipe = 3.96 10 -5 m 3 s -1 and h start p = 15 mm (vortex dune regime). The grey scale represents the height h p of the bed in mm.

Figure 4 .

 4 Figure 4. Poiseuille flow on a wavy bottom.

Figure 5 .

 5 Figure 5. (a) Numerical ωi as a function of k for Re 2D = 35 (♦), Re 2D = 49.3 (•), Re 2D = 54.3 (×), and Re 2D = 70 ( ). The lines represent equation (3.15) for Re 2D = 35 (small-dashed line), Re 2D = 49.3 (dotted line), Re 2D = 54.3 (solid line), and Re 2D = 70 (long-dashed line) and (b) blow-up.

Figure 6 .

 6 Figure6. Phase-Diagram of the dune patterns in the plane Re pipe , Ga(h f /d) 2 . The different regimes represented are: 'no motion' (×), 'flat bed in motion' ( ) and ( ) when outside the domain of validity of the model, 'small dunes' (•), 'vortex dunes' ( ), and 'sinuous dunes' (△). The dashed line is the predicted threshold for particle motion, the horizontal solid line is the predicted instability threshold, and the dotted line indicates the domain of validity of the particle flux law. The dashed-dotted line is the prediction for instability threshold of[START_REF] Charru | Instability of a bed of particles sheared by a viscous flow[END_REF] and the vertical solid line is that of[START_REF] Charru | Ripple formation on a particle bed sheared by a viscous liquid[END_REF].

Figure 7 .

 7 Figure 7. (a) Dimensionless initial dune velocity against equation (3.14) (the solid line corresponds to the slope one) and (b) dimensionless initial wavenumber as a function of Re pipe (the solid line represents the numerical prediction of the most unstable wave-number and the dotted line the instability threshold) for batch A in fluids 2 (•), 3 (△), 4 ( ), and 5 (×).

Table 2 .

 2 Particle characteristics. The particle density, ρ p , was measured with a pycnometer and a fluid of measured density. The particle size distributions were determined with a digital imaging system. The mean diameter is noted d and the error corresponds to one standard deviation.

	Fluid η (cP) ρ f (g/cm 3 )
	1	0.7 ± 1 0.999 ± 0.001
	2	6.7 ± 1 1.016 ± 0.001
	3	8 ± 1 1.018 ± 0.001
	4	10 ± 1 1.022 ± 0.001
	5	12 ± 1 1.023 ± 0.001

Table 3 .

 3 Fluid characteristics at T = 35 • C. The viscosity, η, was measured with a falling ball viscometer and the fluid density, ρ f , with a pycnometer.

Acknowledgement

We would like to thank Y. Forterre and O. Pouliquen for discussions regarding the dune stability analysis, F. Charru and E. J. Hinch for comments and critical reading, and P. Cervetti, S. Martinez, and F. Ratouchniak for technical assistance. Funding from the Institut Français du Pétrole and Agence Nationale de la Recherche (Project Dunes ANR-07-3 18-3892) are gratefully acknowledged. This work is part of the thesis of M. Ouriemi on ' Érosion, transport et instabilités d'un lit de particules dans un tube' (Université de Provence, 2007).


Leïla Perié
email: leila.perie@curie.fr

Ken R Duffy

Lianne Kok

Rob J De Boer

Ton N Schumacher
email: t.schumacher@nki.nl.

The branching point in erythro-myeloid differentiation

Development of mature blood cell progenies from hematopoietic stem cells involves the transition through lineage-restricted progenitors. The first branching point along this developmental process is thought to separate the erythro-myeloid and lymphoid lineage fate by yielding two intermediate progenitors, the common myeloid and the common lymphoid progenitors (CMPs and CLPs). Here, we use single cell lineage tracing to demonstrate that socalled CMPs are highly heterogeneous with respect to cellular output, with most individual CMP yielding either erythrocytes or myeloid cells after transplantation. Furthermore, based on the labeling of earlier progenitors, we show that the divergence between the myeloid and erythroid lineage develops within multipotent progenitors (MPP). These data provide evidence for a new model of hematopoietic branching in which multiple distinct lineage commitments occur in parallel within the MPP pool.

Introduction

During hematopoiesis, hematopoietic stem cells (HSCs) self-renew or differentiate into all blood cell types through successive stages of lineage commitment, a process that has become a prototype of multi-lineage diversification from a stem cell pool. The prevailing model of hematopoiesis predicts a step-by-step process of lineage commitment in which HSCs give rise to multipotent progenitors (MPPs) that subsequently generate intermediate lineage restricted progenitors. Strongest support for this model has come from the identification of the common lymphoid progenitor (CLP) [START_REF] Kondo | Identification of clonogenic common lymphoid progenitors in mouse bone marrow[END_REF]) that produces lymphoid cells (i.e. T and B lymphocytes, NK cells) and the common myeloid progenitor (CMP) [START_REF] Akashi | A clonogenic common myeloid progenitor that gives rise to all myeloid lineages[END_REF] that gives rise to granulocyte-macrophage (GM) progenitors and the megakaryocyte-erythroid (MkE) progenitors. In this model, the lymphoid-myeloid split forms the first step in lineage commitment downstream of MPPs [START_REF] Reya | Stem cells, cancer, and cancer stem cells[END_REF]. However, a number of studies have questioned the CMP-CLP split as the first step of commitment. Reports of unequal output of GM or MkE progenitors by subsets of CMPs [START_REF] Iwasaki | Distinctive and indispensable roles of PU.1 in maintenance of hematopoietic stem cells and their differentiation[END_REF][START_REF] Nutt | Dynamic regulation of PU.1 expression in multipotent hematopoietic progenitors[END_REF][START_REF] Terszowski | Prospective isolation and global gene expression analysis of the erythrocyte colony-forming unit (CFU-E)[END_REF] suggest a separation between the erythroid and myeloid lineages earlier than the CMP stage. Additionally, earlier bypasses towards lymphoid, myeloid (GM) or MkE commitment have been reported in HSCs [START_REF] Benz | Hematopoietic stem cell subtypes expand differentially during development and display distinct lymphopoietic programs[END_REF][START_REF] Dykstra | Long-term propagation of distinct hematopoietic differentiation programs in vivo[END_REF][START_REF] Sanjuan-Pla | Platelet-biased stem cells reside at the apex of the haematopoietic stem-cell hierarchy[END_REF]Yamamoto et al., 2013) and MPPs [START_REF] Adolfsson | Identification of Flt3+ lympho-myeloid stem cells lacking erythro-megakaryocytic potential a revised road map for adult blood lineage commitment[END_REF][START_REF] Cabezas-Wallscheid | Identification of regulatory networks in HSCs and their immediate progeny via integrated proteome, transcriptome, and DNA methylome analysis[END_REF][START_REF] Igarashi | Transcription from the RAG1 locus marks the earliest lymphocyte progenitors in bone marrow[END_REF][START_REF] Lai | Asymmetrical lymphoid and myeloid lineage commitment in multipotent hematopoietic progenitors[END_REF][START_REF] Mansson | Molecular evidence for hierarchical transcriptional lineage priming in fetal and adult stem cells and multipotent progenitors[END_REF][START_REF] Miyawaki | CD41 marks the initial myelo-erythroid lineage specification in adult mouse hematopoiesis: redefinition of murine common myeloid progenitor[END_REF][START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF][START_REF] Pietras | Functionally Distinct Subsets of Lineage-Biased Multipotent Progenitors Control Blood Production in Normal and Regenerative Conditions[END_REF][START_REF] Takano | Asymmetric division and lineage commitment at the level of hematopoietic stem cells: inference from differentiation in daughter cell and granddaughter cell pairs[END_REF]Yoshida et al., 2006). Based on some of these data, a myelo-lymphoid versus erythro-megakaryocytic separation has been proposed as the first commitment step [START_REF] Adolfsson | Identification of Flt3+ lympho-myeloid stem cells lacking erythro-megakaryocytic potential a revised road map for adult blood lineage commitment[END_REF], but other results have also contradicted this model [START_REF] Boyer | All hematopoietic cells develop from hematopoietic stem cells through Flk2/Flt3-positive progenitor cells[END_REF][START_REF] Forsberg | New evidence supporting megakaryocyte-erythrocyte potential of flk2/flt3+ multipotent hematopoietic progenitors[END_REF]. Overall, which lineage commitment occurs first between the lymphoid, myeloid or erythroid lineages remains to be determined. Importantly, most of the results that bear on this question have been obtained using in vitro clonal assays or population-based tracking approaches that can miss in vivo cellular heterogeneity and thereby influence our interpretation of lineage commitment. To address this issue, we here have utilized cellular barcoding technology that allows the in vivo tracking of single cell fate [START_REF] Gerlach | Heterogeneous Differentiation Patterns of Individual CD8+ T Cells[END_REF][START_REF] Gerrits | Cellular barcoding tool for clonal analysis in the hematopoietic system[END_REF][START_REF] Lu | Tracking single hematopoietic stem cells in vivo using high-throughput sequencing in conjunction with viral genetic barcoding[END_REF][START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF][START_REF] Naik | Cellular barcoding: a technical appraisal[END_REF][START_REF] Schepers | Dissecting T cell lineage relationships by cellular barcoding[END_REF] to describe the steps of myeloid and erythroid commitment during successive stages of hematopoietic development.

Results

To understand the clonal output of individual CMPs, we obtained mouse CMPs based on their original description (ckit + Sca1 -CD16/32 low CD34 + ) [START_REF] Akashi | A clonogenic common myeloid progenitor that gives rise to all myeloid lineages[END_REF], labeled these with unique genetic barcodes in a 6 hr in vitro process, and assessed their ability for myelopoiesis (for the GM lineage) and erythropoiesis (for the MkE lineage). At different time points after transplantation of 2⋅10 3 barcode-labeled CMPs into irradiated recipients, CMP-derived myeloid cells (GFP + CD11b + ) and erythroblasts (GFP + Ter119 + ) were sorted and assessed for their barcode identity using PCR and deep-sequencing (Fig S .1.A, C and D, and Supplemental Table 1). Erythroblasts were measured as a proxy for erythropoiesis, as they retain the nucleus required for barcode analysis. Subsequently, barcode-labeled CMPs were assigned to a class of (lineage biased) progenitor, using a hand-tailored classifier that is based on cellular output of each individual progenitor towards the examined lineages [START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF]. Six days after transplantation, we found that 95% of the engrafted CMPs showed highly biased (>99% of cellular output) production of myeloid cells (69±17% of presumed CMPs) or highly biased (>99% of cellular output) production of erythroblasts (26±13% of presumed CMPs) (Fig 1 .A, quality controls Fig S2 .A and B, effect of classification threshold in Fig S .2.D and F). Furthermore, this large pool of myeloid biased or erythroid biased CMPs accounted for 82±13% and 53±27% of total myeloid cell and erythroblast production respectively (Fig 1.B and Fig S2.D). Only 5±4% of CMPs produced appreciable numbers of both myeloid cells and erythroblasts, overall accounting for 11±12% of total myeloid cell reads and 40±27% of total erythroblast reads (Fig 1.B and Fig S2.D). The observation of these high frequency of unioutcome CMPs could not be explained by limits of detection or sampling issues, as determined by the following detection controls: I) the per progenitor output (as measured in read counts) of uni-outcome CMPs towards one specific lineage (i.e. myeloid cells or erythroblasts) was as high as the per progenitor output of bi-outcome CMPs toward the same lineage (Fig S2 .H) and II) the absence of output of uni-outcome CMPs towards the other lineage was reproducible in technical replicates before applying the replicate filter (not shown). In other words, lack of detectable output of a given CMP towards either the myeloid cell or erythroblast compartment was a reproducible feature.

Similarly, when assessed 14 days after transplantation, 74% of the engrafted CMP produced essentially only myeloid cells (35±3% of presumed CMPs) or only erythroblasts (39±12% of presumed CMPs) (Fig 1 .C and Fig S2 .E), accounting for 74±16% and 53±35% of total myeloid cell and erythroblast production, respectively (Fig 1.D and Fig S2.G). Also at this point in time, only a minority of CMPs showed substantial production of both myeloid cells and erythroblasts. Finally, heterogeneity in output toward either the myeloid and erythrocyte lineage is not specific to the spleen, but was also observed in the bone marrow, where 90% of the engrafted CMP also produced essentially only myeloid cells (35±9% of presumed CMPs) or essentially only erythroblasts (55±12% of presumed CMPs) (Fig 1.G,and Fig S.4.D and E). The observation of the somewhat higher fraction of CMPs with detectable output in both the myeloid and erythrocyte lineage over time (Day 6: 5±4%; Day 14: 26±15% of presumed CMPs) may either be interpreted as a more prolonged engraftment potential of the small subset of CMPs producing both erythroblasts and myeloid cells or contamination of CMPs by Lin -Sca1 + ckit + (LSK) cells (Fig S2 .C), whose contribution may also be expected to increase over time after transplantation. A preponderance of uni-output was also observed at two CMP doses tested (Fig S.4.F and G).

In addition to heterogeneity in output towards the myeloid and erythrocyte lineage, further disparity in cellular output between individual CMPs was observed at day 14 after transplantation (Fig 1.E and Fig S.3), with 16±4% of CMPs producing detectable output only in the dendritic cell (DC) lineage (Fig S.4.A, note that DC production could not be assessed at day 6). These DC uni-outcome CMPs accounted for 24±25% of total DC production (Fig S .4.B) and their presence could not be explained by limits of detection (Fig S.4.B). While these data do not exclude the possibility that some DC uni-outcome CMPs might have produced myeloid cells or erythroblasts at an earlier point in time, this would then imply that individual CMPs differ in the kinetics with which they contribute to DC versus myeloid/erythroid lineage. Collectively, these data show a striking heterogeneity within the CMP pool, greater than previously suggested by in vitro or in vivo experiments that examined output at the population level [START_REF] Iwasaki | Distinctive and indispensable roles of PU.1 in maintenance of hematopoietic stem cells and their differentiation[END_REF][START_REF] Nakorn | Characterization of mouse clonogenic megakaryocyte progenitors[END_REF][START_REF] Nutt | Dynamic regulation of PU.1 expression in multipotent hematopoietic progenitors[END_REF][START_REF] Terszowski | Prospective isolation and global gene expression analysis of the erythrocyte colony-forming unit (CFU-E)[END_REF]. Importantly, this heterogeneity could not be explained by detection issues and it was relatively insensitive to the threshold used for output classification .

To determine at which developmental point the observed divergence between the myeloid and erythroid lineage occurs, we subsequently generated barcode-labeled HSCs (LSK + CD150 + ) and MPPs (LSK + CD150 -Flt3 + , excluding LMPPs) (Fig S1 .B) and, after transplantation, analyzed their production of myeloid cells and erythroblasts, together with their lymphoid cell production (measured by analysis of CD19 + B cells), and CD11c + dendritic cell production (Fig S1 .D). As previously described using cellular barcoding [START_REF] Gerrits | Cellular barcoding tool for clonal analysis in the hematopoietic system[END_REF][START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF], a small number of HSC were the major contributors to total cellular output at day 27 after transplantation, and these HSC were multi-outcome (Fig 2.A and B,and Fig S6.A). With respect to contribution towards the erythroid and myeloid lineage, >99% of the total cell output of myeloid cells and erythroblasts was derived from HSCs (45±7%) that produced both cell types at this time point ( .C). Furthermore, the occurrence of these biased HSCs may well be explained by insufficient sensitivity of detection, as these HSCs have a per progenitor output (as measured in read counts) towards the erythroid or the myeloid lineage that is several fold lower than that of multi-outcome HSCs for the same lineage (Fig S .5.D). Also, when assessed at day 42 after transplantation (Fig S .6.B), the vast majority of myeloid cells and erythroblasts (≥90%) were derived from HSCs that contributed to both lineages. In conclusion, the large majority of erythroid and myeloid cells are produced by multi-outcome HSCs.

In contrast to this, substantial variability in the ability to contribute to erythroid and myeloid cell production was apparent amongst MPPs (Fig 3 .A andB). At Day 14 after transplantation, only 20±9% of the MPPs produced both erythroblasts and myeloid cells, irrespective of their capacity to produce dendritic cells or B cells (Fig 3.C and Fig S.8, quality controls and effect of classification threshold in Fig S .7.A and B). While these cells still contributed to a large fraction of erythroblast output (84±23%), approximately 50% of total myeloid output was derived from the large pool of MPPs (51±9%) that produced myeloid cells without producing erythroblasts (Fig 3.D and Fig S7.C). Furthermore, contrary to the uni-outcome biased HSCs, the observation of these myeloid-not-erythroid biased MPPs could not be explained by limits of detection, as I) these cells consistently had a per progenitor output (as measured in read counts in myeloid cells) as high as the per progenitor output (read counts in myeloid cells) of MPP producing both myeloid cells and erythroblasts (Fig S.7.D); II) their lack of output in erythroblasts was reproducible in technical replicates before applying the replicate filter (not shown). Note that we also detected a significant fraction of MPPs that only show dendritic cell output (Fig 3 .B), as has also been observed for LMPPs [START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF]. These data show the existence of a subpopulation of true common progenitors to myeloid and erythroid lineages within the MPP pool, potentially within the recently described CD41 + cell population [START_REF] Miyawaki | CD41 marks the initial myelo-erythroid lineage specification in adult mouse hematopoiesis: redefinition of murine common myeloid progenitor[END_REF]. However, within the MPP population, divergence of a myeloid biased subpopulation is already apparent.

Discussion

The main finding of this work is that the population classically defined as common myeloid progenitors is highly heterogeneous, with most individual CMPs only yielding erythrocytes or only yielding myeloid cells after transplantation, and with some CMPs contributing towards the DC lineage only. A minor subpopulation is observed that could represent a rare but true subset of bipotent CMP, but may also reflect a contamination of LSK+ cells.

As the existence of a "common" erythro-myeloid progenitor has been the basis for the model in which the first step of commitment is formed by the separation between CMP and CLP, our findings suggest a revised model of hematopoiesis. Miyawaki et al. have recently identified a more immature common myeloid progenitor within the MPP population [START_REF] Miyawaki | CD41 marks the initial myelo-erythroid lineage specification in adult mouse hematopoiesis: redefinition of murine common myeloid progenitor[END_REF]. Existence of such a common erythro-myeloid progenitor is still consistent with a model in which the first step of commitment is formed by the myeloid-lymphoid split. Likewise, our data also show that a true common erythro-myeloid progenitor exists within the MPP population. Importantly though, 40% of the total output of myeloid cells was derived from progenitors that neither produced erythrocytes or lymphoid cells (Fig S.7.C). These data indicate that, next to the classical erythro-myeloid versus lymphoid split, an additional myeloid-only commitment already manifests itself at the MPP stage.

Alternative models for hematopoiesis have been proposed in which the megakaryocyte-erythroid lineage branches first, followed by a myeloid-lymphoid split [START_REF] Adolfsson | Identification of Flt3+ lympho-myeloid stem cells lacking erythro-megakaryocytic potential a revised road map for adult blood lineage commitment[END_REF][START_REF] Arinobu | Reciprocal activation of GATA-1 and PU.1 marks initial specification of hematopoietic stem cells into myeloerythroid and myelolymphoid lineages[END_REF]. In this model, direct erythropoiesis from HSC or MPP is predicted to occur independently of myelopoiesis and lymphopoiesis. We did not observe evidence for a significant role of such an erythroid bypass from either HSC or MPP. Specifically, of total erythroblast production from HSC or MPP respectively, only 0.1±0.1% and 9±0.01% was derived from cells that did not yield myeloid or B cells. Elegant recent in vivo lineage tracing analysis suggests that transplantation can skew output of HSCs towards a few dominant clones [START_REF] Sun | Clonal dynamics of native haematopoiesis[END_REF]. However, at least under transplantation conditions, an early erythroid commitment (Yamamoto et al., 2013) is at best a minor pathway.

The experimental measurements that we make in our work are the in vivo output of single cells. As such, the strong bias in single cell output that we observe could result from the commitment pre-injection of the cells that are used for reconstitution, or a dominant effect of the niche in which an individual progenitor cell finds itself after transplantation. We strongly believe that the biased output that we observe for individual progenitor cells originates from pre-commitment for two reasons. First, a pre-commitment of for instance individual CMP to produce either erythrocytes or myeloid cells implies an 'encoding' of this commitment at the genetic level, and the observation of substantial transcriptional heterogeneity between individual cells within the murine CMP pool by Amit and coworkers (I. Amit, pers. comm.) provides strong support for such a model. Second, if the commitment towards the production of a single type of output is pre-existing rather than induced by the niche that is encountered, one would expect to also observe heterogeneous behavior under in vitro conditions in which individual cells all experience identical conditions. In support of this, Notta and colleagues have provided strong evidence that human bone marrow CMPs and MPPs are highly heterogeneous and composed of subpopulations with uni-outcome myeloid, erythroid, and megakaryocyte potential in in vitro assay systems [START_REF] Notta | Distinct routes of lineage development reshape the human blood hierarchy across ontogeny[END_REF].

Collectively with these data, our single cell tracing data provide evidence for a model in which multiple combinations of lineage commitments (e.g. erythroid-myeloid, myelo-lymphoid) are already manifest during the transition of the MPP stage of hematopoietic development. The combination of single cell transcriptomics and single cell lineage tracing, preferably under conditions of unperturbed hematopoiesis [START_REF] Sun | Clonal dynamics of native haematopoiesis[END_REF], Busch 2015) should help reveal how sibling cells commit to different types of output at this stage.

Experimental procedures:

Mice, progenitor isolation and barcode labeling C57BL/6J (CD45.2) donor mice, and C57BL/6 Pep3b (CD45.1) recipient mice irradiated with 5 gy were used in all experiments. Labeling of cells with the lentiviral barcode library was performed as described in [START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF]. In brief, BM was harvested from femurs, tibias and ilia and enriched using anti-CD117 magnetic beads (Miltenyi). The c-kit + fraction was stained with antibodies against CD16/32 (FITC, clone 24G2, BD Pharmingen or PerCPCy5.5, clone 24G2, BD Biosciences), CD117 (c-kit APC, clone 2B8, Biolegend), Sca-1 (Pacific Blue, clone D7, Biolegend), CD135 (Flt3 PE, clone A2F10, ebiosciences), CD34 (alexa700, clone RAM34, ebioscience) and CD150 (Slam PEcy7, clone TC15-12F12.2, Biolegend). MPPs, CMPs and HSCs were sorted using a strategy similar to that described previously [START_REF] Akashi | A clonogenic common myeloid progenitor that gives rise to all myeloid lineages[END_REF]

(Fig S.1.B & C).
Per experiment, up to 10 5 progenitors were transduced for 6 h in stem cell medium (stempan SFEM, Stem Cell Technologies) with 50ng/ml stem cell factor (SCF). For the first 90 min of transduction, progenitors were centrifuged at 200rpm (low brake) and cells were subsequently incubated at 37°C for 4.5 hours. Following transduction, cells were washed and injected intravenously into at least two recipient mice in saline solution. Efficiency of transduction ranged between 10-40%.

Cell isolation

Cell suspensions were derived from spleen or bone marrow of killed mice and positively selected for Ter119 expression using biotinylated anti-Ter119 antibody (BD Pharmingen) and anti-biotin beads (Miltenyi). The Ter119 -fraction was stained at 4°C using a combination of antibodies 1.

Polymerase chain reaction and sequencing

Sorted cell samples were lysed in Viagen buffer (direct PCR, Viagen) and split into technical replicates [START_REF] Schumacher | Mapping the life histories of T cells[END_REF]. Barcodes were amplified by nested polymerase chain reaction (PCR) as previously described in [START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF]. In brief, barcode sequences were first amplified using top-LIB (5'TGCTGCCGTCAACTAGAACA-3') and bot-LIB (5'GATCTCGAATCAGGCGCTTA-3') primers (30 cycles: 15'' 95°C; 15'' 57.2°C; 15'' 72°C). PCR products were then subjected to a 2 nd PCR (30 cycles: 5'' 94°C; 5'' 57.2°C; 5'' 72°C) in which the P5 and P7 (5'-CAAGCAGAAGACGGCATACGAGAT-3') adaptors required for sequencing on Illumina platforms, as well as a unique 8 bp index and a sequencing primer annealing site were attached. Primers used were as follows: forward (5'-P5-seq_prim-index-CAGGCGCTTAGGATCC-3'), reverse (5'-P7-TGCTGCCGTCAACTAGAACA-3'). Up to 192 tagged PCR products were pooled and sequenced by next-generation sequencing using a HiSeq2000 or HiSeq2500 platform (Illumina). Total read lengths were 50-60 bp, of which the first 22 bp consisted of the index and constant regions. The following 15 bp were used to distinguish between different barcodes.

Filtering procedure

Sequence reads were first screened for quality by the software packages that are part of the Illumina pipeline and only high-quality reads were selected (~10-70x10 6 per run). Subsequently, only reads showing a 100% match to the expected scheme: index-CAGGCGCTTAGGATCC-random_sequence were selected (typically >90% of reads), whereby a 100% match to one of the 192 index sequences was required. Of the random sequence (the barcode), only the first 15 bp were considered. To filter out reads that reflect PCR or sequencing errors within the 15 bp barcode region, all filtered reads were compared to a barcode reference list [START_REF] Gerlach | Heterogeneous Differentiation Patterns of Individual CD8+ T Cells[END_REF], and only sequences that showed a 100% match to one of the sequences within the barcode reference file were maintained (typically 80-93%). Using a customized script in R [START_REF] Team | R: A language and environment for statistical computing[END_REF], we subsequently applied a filtering procedure composed of 3 steps, similar to the procedure used in [START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF]. Specifically, samples for which insufficient read counts were obtained during deep sequencing (average of the two technical replicates less than 10 3 , compared to the expected 10 5 ) were first excluded. Samples having passed this step were then normalized to 10 5 reads for each sample. To control the quality of the analysis of the barcode content of a given sample, we additionally excluded samples where the two technical replicates displayed a Pearson correlation coefficient of <0.8. These technical replicate controls are essential to give confidence whether barcode expansion and cell type recovery were sufficient for informative lineage relationship analysis. Thirdly, after passing these technical replicate quality controls, all reads of barcodes present in only one of the two technical replicates of a given sample -an indication that there was a low confidence for inclusion of that barcode -were changed to zero (0) reads for that sample. Resulting data were collated in a table of normalized barcode read counts in each cell type.

Barcode analysis and barcode classification

For heatmap analysis, read counts were transformed using the hyperbolic arcsine function that is similar to a logarithmic function, but can accommodate barcode reads with a value of 0. After confirming that technical replicates yielded similar data, as described above and by visual inspection (Fig S .2.A,S.5.A and S.7.A), their average was taken for further analysis. To classify individual progenitors by their lineage bias, we used a previously published hand tailored classifier [START_REF] Naik | Diverse and heritable lineage imprinting of early haematopoietic progenitors[END_REF]. In summary, an additional normalization per progenitor was applied in each mouse, thereby enabling categorization into classes of biased output towards the analyzed lineages (Myeloid cells, Erythroblasts, Dendritic cells or B cells depending on the experiment). Potency to contribute to a given lineage was assigned if more than 1% of the output of a given CMP, MPP, or HSC was towards that lineage (Fig 123). Other thresholds were also analyzed (Fig S .2, S.4, S.5,S.7).

Quality controls

To illustrate the quality of the samples, the barcodes find in both duplicates of the same sample are always compared and representative plots of each experiment are provided in the supplemental figures (Fig S .2, S.4, S.5,S.7). Furthermore, by cross-comparing two mice that received barcoded progenitors from the same transduction batch, the extent of 'repeat use'instances where two progenitors are transduced with virus particles that harbor the same barcode -was always compare and was found to be minimal (Fig S .2, S.4, S.5,S.7).

Uni-outcome progenitors were further analyzed for potential detection issues by: I). Determining whether the distribution of clonal outputs of uni-outcome progenitors towards one specific lineage was similar to the clonal output towards the same lineage of multi-outcome progenitors. II). Determining whether the lack of output of uni-outcome progenitors towards the other lineages was reproducible in technical replicates before applying the replicate filter.

Yamamoto, R., Morita, Y., Ooehara, J., Hamanaka, S., Onodera, M., Rudolph, K.L., Ema, H., and Nakauchi, H. (2013). Clonal analysis unveils self-renewing lineage-restricted progenitors generated directly from hematopoietic stem cells. Cell 154, 1112-1126. Yoshida, T., Ng, S.Y., Zuniga-Pflucker, J.C., and Georgopoulos, K. (2006). Early hematopoietic lineage restrictions directed by Ikaros. Nature immunology 7, 382-391. Based on their relative output to the different cell types analyzed, barcodelabeled CMPs were assigned to a class of (lineage biased) progenitor using a hand-tailored classifier (see methods). The effect of different thresholds for lineage bias classification is depicted for the proportion of CMPs assigned to each class, for output observed at day 6. (E). Same as D for CMPs output observed at day 14. (F). Quantitative contribution of the classes to the myeloid lineage (left panel) and erythroid lineage (right panel) for different thresholds (x-axis) for CMP output observed at day 6 (mean)/ (G). Same as (F) for CMP output observed at day 14. '0', '1', and '2.5' refer to the minimal percentage of total reads of a progenitor that is observed within a given cell type in order to classify this progenitor as contributing to this lineage. As an example, a CMP of 98% of the output was in E and 2% was in M would be classified as E biased using the 2.5% classification threshold, but not the 0% or 1% classification thresholds. Note that a large fraction of CMPs produced >99% of its output in either the myeloid or the erythroid lineage. (H). Cumulative read distribution of barcoded CMPs pooled from 4 mice classified as uni-outcome (only producing myeloid cells (M, left panel) or erythroblasts (E, right panel)) versus barcoded CMPs classified as bi-outcome (producing both E and M, M+E, with the E and M read distributions respectively) at day 6. Classification was done with a threshold of 1%. The graph is representative of one experiment and shows that the read distribution of uni-outcome CMPs was similar to that of the bi-outcome CMPs. (I) Same as (H) for day 14. Showing that also at this time point the cumulative read distribution of uni-outcome CMPs was similar to that of bi-outcome CMPs. Together with sampling controls as shown in panel A, these data indicate that the frequent observation of uni-outcome CMPs is not explained by the low cellular output of these cells.

!"# $!"# %!"# &!"# '!"# (!!"# ($!"# (%!"# )# *# *+)# )# *# B !"# $!"# %!"# &!"# '!"# (!!"# ($!"# (%!"# )# *# *+)# )# *# E 39% % output

B cells Erythroblasts Myeloid cells Dendritic cells

!"# $!"# %!"# &!"# '!"# (!!"# ($!"#
# !"# $!"# %!"# &!"# '!"# (!"# )!"# *!"# +!"# ,!"# $!!"# !# $# %-(# !"#$%&$%"'#(%)*+$,#("%#"-" ./0# 0# .# !"# $!"# %!"# &!"# '!"# (!"# )!"# *!"# +!"# ,!"# $!!"# !# $# %-(#
!"#$%&$%"'#(%)*+$,#("%#"-"

  Fig 2.C and D, quality controls and effect of classification threshold Fig S.5.A-C). Some HSCs produced only myeloid cells (10±4%) or erythroblasts (12±7%), without detectable production of any of the other cell types (Fig 2.B). However, these HSCs collectively contributed to less than 0.1% of the total cell output of erythroblasts and myeloid cells (Fig 2.D and Fig S.

  5

  against CD45.2 (Pacific blue, clone 104, Biolegend), CD11C (APC, clone HC3, BD biosciences), CD11b (PercPCy5.5, clone M1/70, ebioscience), CD19 (APC-Cy7, clone 1D3, BD Pharmingen), Ter119 (PEcy7, BD Pharmingen) to sort Myeloid cells, B cells and Dendritic Cells (Fig S1.D). The Ter119 + fraction was stained at 4°C with antibody against Ter119 (Pecy7, BD Pharmingen) to sort the erythroblast fraction (Fig S1.D). Output cell types were sorted following the gating strategy in Fig S1.D and cell numbers are depicted in Supplemental Table

Fig. 1 .

 1 Figure legends:

Fig. 2 .

 2 Fig. 2. Common origin of Erythrocytes and myeloid cells within HSC Barcode-labeled CD150 + HSCs were transplanted into 5gy-irradiated mice. At day 27 posttransplantation, Myeloid cells (CD11b + ), Erythroblasts (Ter119 + ), B cells (CD19 + ) and dendritic cells (CD11c + ) were recovered from spleen and barcode content was analyzed. (A) Heatmap representation of the output of individual CD150 + HSC (rows) to different cell types (columns), arcsin transformed data. (B) Proportion of HSCs classified using a 1% classification threshold in the different possible categories (averaged over mice). (C) Same analysis as in (B) but examining erythroblast and myeloid cell output only (D) Quantitative contribution of the different classes to Myeloid cell and Erythroblast production (mean+sd). All the results comes from 4 mice/experiment. Quality controls and other thresholds in Fig S5.

Fig. 3 .

 3 Fig. 3. Divergence between erythrocyte and myeloid output within the MPP pool Barcode-labeled MPPs were transplanted into 5gy-irradiated mice. At day 14 posttransplantation, Myeloid cells (CD11b + ), Erythroblasts (Ter119 + ), B cells (CD19 + ) and dendritic cells (CD11c + ) were recovered from spleen and barcode content was analyzed. (B) Proportion of MPPs classified using a 1% classification threshold in the different possible categories (average over mice). (C) Same analysis as (B) but examining erythroblasts and myeloid cell output only (D) Quantitative contribution of the classes to Myeloid cell and Erythroblast production from

Figure

  Figure 1

  Figure 2 A

  Figures S1-S8Legends of Supplemental Figures Fig S.1. Summary of experimental setup and sort profiles of progenitors and output cell types (A). Progenitors (HSC, MPP or CMP) were transduced with low efficiency (5-15%) to minimize the chance of 2 different barcode virus particles integrating into the same progenitor. Progenitors from a single transduction batch were then transplanted into 2-8 5gy irradiated mice (to allow later comparison of barcodes between mice as a control for random barcode sharing). At the indicated time points after transplantation, Myeloid cells (CD11b + ), Erythroblasts (Ter119 + ), Dendritic Cells (CD11c + ) and B cells (CD19 + ) were isolated from the spleen of recipient mice. Each sample was split into two replicates; barcodes were amplified by PCR and analyzed by deep-sequencing. (B). Bone marrow progenitors were enriched using anti-CD117 (c-kit) magnetic beads and subsequently stained with antibodies against CD16/32, CD117, Sca-1, CD135 (Flt3), CD34 and CD150. MPP and HSC were sorted using the gates shown. (C). Gating strategy used to isolate CMP. (D). Splenic cells were separated into a Ter119 + and Ter119 -fraction by magnetic bead-based selection. The Ter119 -fraction was used to sort B cells, Myeloid cells and Dendritic Cells as shown. The Ter119 + fraction was used to sort Erythroblasts as shown.

  Fig S.1. Summary of experimental setup and sort profiles of progenitors and output cell types(A). Progenitors (HSC, MPP or CMP) were transduced with low efficiency (5-15%) to minimize the chance of 2 different barcode virus particles integrating into the same progenitor. Progenitors from a single transduction batch were then transplanted into 2-8 5gy irradiated mice (to allow later comparison of barcodes between mice as a control for random barcode sharing). At the indicated time points after transplantation, Myeloid cells (CD11b + ), Erythroblasts (Ter119 + ), Dendritic Cells (CD11c + ) and B cells (CD19 + ) were isolated from the spleen of recipient mice. Each sample was split into two replicates; barcodes were amplified by PCR and analyzed by deep-sequencing. (B). Bone marrow progenitors were enriched using anti-CD117 (c-kit) magnetic beads and subsequently stained with antibodies against CD16/32, CD117, Sca-1, CD135 (Flt3), CD34 and CD150. MPP and HSC were sorted using the gates shown. (C). Gating strategy used to isolate CMP. (D). Splenic cells were separated into a Ter119 + and Ter119 -fraction by magnetic bead-based selection. The Ter119 -fraction was used to sort B cells, Myeloid cells and Dendritic Cells as shown. The Ter119 + fraction was used to sort Erythroblasts as shown. Fig S.2. Quality controls of CMP barcode data (related to Fig 1) (A). CMPs from a single transduction batch were transferred into 2-4 mice per experiment. The barcode repertoire observed in technical replicates of CMP-derived myeloid cells is depicted, both within and between mice. Plots representative of the quality controls for accuracy of barcode sampling (comparison of barcodes seen in technical replicates of a mouse, left panel) and of the inferred chance of random sharing of barcodes between different progenitors in one mouse (comparison of barcodes observed in two mice that received CMPs from the same transduction batch, right panel). (B). Quantification of average repeat usage of barcodes. The proportion of barcodes shared between mice (m1 & m2) or present in only one of the mice (m1 or m2) from the same transduction batch is averaged over 3 different experiments at day 6 (left panel) or 2 experiments at day 14 (right panel). (C). Possible contamination of CMP by LSK cells. The first and the second panel show the definition of the CMP gate in order to exclude LSK, and later MEPs and GMPs. The third panel shows the position of LSK cells when not excluded by the gate in the first panel, revealing a potential source of contamination at the border between the LS-K and LSK gates as indicated by the red arrow. One representative experiment shown. (D).Based on their relative output to the different cell types analyzed, barcodelabeled CMPs were assigned to a class of (lineage biased) progenitor using a hand-tailored classifier (see methods). The effect of different thresholds for lineage bias classification is depicted for the proportion of CMPs assigned to each class, for output observed at day 6. (E). Same as D for CMPs output observed at day 14. (F). Quantitative contribution of the classes to the myeloid lineage (left panel) and erythroid lineage (right panel) for different thresholds (x-axis) for CMP output observed at day 6 (mean)/ (G). Same as (F) for CMP output observed at day 14. '0', '1', and '2.5' refer to the minimal percentage of total reads of a progenitor that is observed within a given cell type in order to classify this progenitor as contributing to this lineage. As an example, a CMP of 98% of the output was in E and 2% was in M would be classified as E biased using the 2.5% classification threshold, but not the 0% or 1% classification thresholds. Note that a large fraction of CMPs produced >99% of its output in either the myeloid or the erythroid lineage. (H). Cumulative read distribution of barcoded CMPs pooled from 4 mice classified as uni-outcome (only producing myeloid cells (M, left panel) or erythroblasts (E, right panel)) versus barcoded CMPs classified as bi-outcome (producing both E and M, M+E, with the E and M read distributions respectively) at day 6. Classification was done with a threshold of 1%. The graph is representative of one experiment and shows that the read distribution of uni-outcome CMPs was similar to that of the bi-outcome CMPs. (I) Same as (H) for day 14. Showing that also at this time point the cumulative read distribution of uni-outcome CMPs was similar to that of bi-outcome CMPs. Together with sampling controls as shown in panel A, these data indicate that the frequent observation of uni-outcome CMPs is not explained by the low cellular output of these cells.

Fig S. 3 .

 3 Fig S.3. Heatmaps of CMP contribution to different cell types for individual mice (related to Fig 1). Each heatmap represents the contribution of individual CMPs (rows) to the indicated cell types

Fig S. 7 .

 7 Fig S.7. Quality controls of MPP data (related to Fig.3) (A). As in Fig.S.2, plots are representative of the quality controls for accuracy of barcode sampling (left panel) and of the inferred chance of having different progenitors with the same barcode in one mouse (repeat use), by comparing barcode content in two mice that received MPPs from the same transduction batch (right panel). (B). The effect of different value of the threshold (t) used for lineage bias classification is shown for the proportion of MPPs assigned to each class, for output recovered at day 14. (C). Quantitative contribution of the classes to the myeloid lineage (left panel) and erythroid lineage (right panel) for different thresholds (mean) for MPPs at day 14 Here the classification of MPPs was performed using either only erythroblasts and myeloid cells output or output for all cell types analyzed (erythroblast, myeloid, dendritic and B cells). The results shown in (C) are from the latter classification. Color-coding indicates the correspondence between the two classifications: shade

Fig S. 8 .

 8 Fig S.8. Heatmaps of MPP contribution to different cell types for individual mice (related to Fig 3). Each heatmap represents the contribution of individual MPPs (rows) to the indicated cell types (columns) of each mouse shown in Fig.3, day 14 after transplantation. Arcsin transformed data clustered by complete linkage using Euclidean distances.

  Figure S5

  Figure S8

  293 62 60 167 59 308 122 307 173 309 61 170 362 361 234 161 85 73 254 252 248 253 264 267 11 12 271 249 251 247 266 273 246 152 242 192 355 193 204 183 278 206 210 97 53 94 92 182 187 26 219 341 180 337 356 336 203 319 321 330 208 347 217 310 323 191 334 335 45 303 313 77 314 47 95 349 188 324 83 348 338 325 225 72 89 227 68 320 100 207 213 353 87 270 101 222 216 189 19 238 233 38 50 119 118 27 29 286 236 256 239 241 277 257 258 255 288 292 282 281 280 263 272 245 274 260 240 275 276 358 112 108 230 231 346 312 315 265 200 333 185 190 139 212 71 107 117 106 259 294 34 284 235 237 42 114 110 37 105 8 8699 298 289 296 88 283 345 46 184 148 156 290 360 104 109 116 57 243 302 115 141 129 123 133 128 124 120 137 154 166 144 145 136 5 41 52 155 30 21 25 23 58 4 13 153 279 174 285 16 3 17 15 43 9 28 54 22 20 44 318 18 159 125 81 291 24 138 351 96 316 51 198 158 48 195 326 90 179 209 211 202 297 39 262 221 223 177 205 226 304 93 201 215 33 82 299 64 78 75 224 350 80 331 56 186 301 300 218 32 343 352 329 79 103 194 342 344 214 151 98 340 332 74 91 66 84 69 70 295 339 196 76 67 322

63 1 36 146 287 268 163 261 162 140 35 178 130 40 199 220 49 160 311 102 14 2 244 181 250 157 126 6 135 143 10 55 142 134 147 121 131 132 149 127 229 111 228 164 113 232 359 171 357 31 269 354 327 328 197 317 305 7 150 165 65 306 176 172 169 175 168

  9 3 20 46 27 77 85 75 50 110 70 113 28 116 76 64 86 87 1 63 56 21 58 109 60 90 91 99 81 54 98 101 103 59 102 111 88 79 112 29 83 108 107 97 84 95 96 115 106 78 80 89 62 73 72 69 82 93 94 38 44 65 52 67 68 66 24 92 61 2 43 42 104 57 48 55 18 53 49 119 100 114 37 117 118 105 71 74 120 121 45 36 16 17 14 47 51 41 40 39 22 31 15 4 33 19 34 35 13 11 7 12 26 23 30 8 10 6 25 32 5 56 49 48 54 51 3 52 78 59 23 15 28 2 12 13 111 112 30 8 110 109 45 113 119 114 97 70 117 121 120 122 17 36 34 55 65 29 37 116 118 20 93 22 10 414 9 5 1 18 35 47 25 26 81 66 71 87 74 85 64 53 50 105 67 77 32 99 75 95 33 73 84 89 98 91 100 61 63 68 72 108 106 69 103 88 94 115 42 76 19 40 6 39 79 60 90 92 83107 82 96 102 104 58 86 27 44 43 21 11 101 41 80 62 38 7 16 24 33 32 28 39 36 29 34 35 27 14 55 3 50 112 107 92 81 70 56 47 52 22 31 25 10 20 2 124 19 123 125 8 12 13 11 15 5 18 16 41 26 4 1 84 115 76 46 121 65 105 117 108 77 74 61 71 118 111 82 114 120 42 80 90 66 110 100 88 91 85 73 119 79 44 68 60 93 94 72 95 98 109 89 122 63 97 67 86 102 116 101 45 87 75 106 9 103 104 7 126 21 54 6 23 17 40 30 64 78 99 57 62 83 49 96 48 51 43 59 53 113 69 58

	B cells 57 31 46 B cells erythroblasts myeloid cells dendritic cells dendritic cells erythroblasts myeloid cells 38 37 24
				65 59 61 62 60 113 111 55 31 39 33 82 78 75 64 79 103 56 80 32 98 74 91 76 67 66 84 69 70 112 108 71 107 117 106 34 46 88 8 86 99 42 114 110 37 105 115 104 109 116 57 19 38 50 119 118 27 29 26 97 53 94 92 83 95 101 100 45 77 47 68 89 72 87 7 11 12 23 58 4 13 5 41 52 30 21 25 16 3 17 15 43 9 28 54 22 20 44 36 63 1 14 2 6 10 49 96 102 40 90 85 73 35 81 24 18 93 51 48
	dendritic cells	B cells	erythroblasts	myeloid cells

Acknowledgments:

The authors thank A. Pfauth and F. van Diepen for cell sorting, J. Urbanus for lentivirus production, A. Velds and R. Kluin for computational assistance and S. Naik for discussion. The data and the code presented in this manuscript are available upon request. The authors declare no conflict of interest. This work was supported by Human Frontier Science Program grant RGP0060/2012 (to T.S. and K.D.). Netherlands Organization for Scientific Research grant TOPGO.L.10.042 (to T.S. and R.d.B.), ERC grant Life-His-T (to T.S.), and Science Foundation Ireland grant 12 IP 1263 (to K.D.), an ATIP-Avenir grant from CNRS and Bettencourt-Schueller Foundation (to L.P.) and two grants from the Labex CelTisPhyBio (No. ANR-10-LBX-0038) and Idex Paris-Science-Lettres Program (ANR-10-IDEX-0001-02 PSL) (to L.P.).