Symmetry breaking in the embryonic skin triggers directional and sequential plumage patterning

The development of an organism involves the formation of patterns from initially homogeneous surfaces in a reproducible manner. Simulations of various theoretical models recapitulate final states of natural patterns, yet drawing testable hypotheses from those often remains difficult. Consequently, little is known about pattern-forming events. Here, we surveyed plumage patterns and their emergence in Galliformes, ratites, passerines, and penguins, together representing the three major taxa of the avian phylogeny, and built a unified model that not only reproduces final patterns but also intrinsically generates shared and varying directionality, sequence, and duration of patterning. We used in vivo and ex vivo experiments to test its parameter-based predictions. We showed that directional and sequential pattern progression depends on a species-specific prepattern: an initial break in surface symmetry launches a travelling front of sharply defined, oriented domains with self-organising capacity. This front propagates through the timely transfer of increased cell density mediated by cell proliferation, which controls overall patterning duration. These results show that universal mechanisms combining prepatterning and self-organisation govern the timely emergence of the plumage pattern in birds.


Homogeneous steady states
Spatially homogeneous equilibria are solutions of the equation (1) (n(t, r), u(t, r), v(t, r)) that are independent of time (steady state, or equilibrium) and of the position r on the tissue (spatially homogeneous). These spatially homogeneous solutions are thus solutions of the system of algebraic equations: that has two relevant solutions: • the trivial state (0, 0, 0), corresponding to an unrealistic situation of complete absence of cells in the tissue (n = 0), in turn leading to an absence of activator and inhibitor u = v = 0. Because of the proliferation term, we expect this solution to be always unstable: a solution of equation (1) starting with a small number of cells will be able to proliferate at a rate equal to α n (while n is small compared to the tissue capacity) and thus the system will move away from that steady state.
• A non-trivial equilibrium (n s , u s , v s ) where the number of cells reaches the carrying capacity of the tissue n s = β n , and v s = αvβn δv u s 2 and u s is a solution to the degree five polynomial: This degree-five polynomial can have up to 5 solutions. However, we found numerically that only one solution is real, and the other four are formed of two pairs of complex conjugated numbers. Only the real solution is relevant for our study.

Linear Stability -General Analysis
To assess stability of the equilibria found, we consider whether small perturbations about one of solution are initially damped or amplified. When small perturbations vanish, the equilibrium is said to be stable, and it is otherwise unstable.
Consider a small initial perturbation of the equilibrium of size ε, and let us denote the solution n = n 0 + εφ n , u = u 0 + εφ u and v = v 0 + εφ v for some functions (φ n , φ s , φ v ) to be determined and (n 0 , u 0 , v 0 ) one of the spatially homogeneous solutions computed above. Using the fact that n is solution to the equation, we obtain: and using the fact that u 0 is spatially homogeneous (so the gradient is 0), that n 0 is spatially homogeneous (thus treated as a constant for the divergence in the chemotaxis term) and neglecting terms of order ε 2 (and possibly higher), one obtains the linearized equation on φ n : Proceeding in the same fashion for the other two variables, we obtain the linearized system: of the linear perturbation equation decays to zero in time. To this purpose, we will decompose the initial shape of the perturbation on the Fourier basis, which is given by the functions W k 1 ,k 2 = cos k 1 πx L 1 cos k 2 πy L 2 to ensure that our boundary conditions are satisfied. Mathematically, this set of functions is said to form a Hilbert basis (in the L 2 sense) of the set of functions Ω such that ∂ and any initial perturbation can be decomposed univocally in this Fourier basis: Stability analysis thus amounts to investigating whether Fourier modes are amplified and damped. If at least one Fourier mode is amplified, the solution is unstable. The choice of decomposing the initial condition on the Fourier basis is relevant to study equation (5) because Fourier modes are eigenfunctions of the diffusion operator: ∆W k 1 ,k 2 = µ k 1 ,k 2 W k 1 ,k 2 associated with the eigenvalues: Moreover, because the equation (5) is linear, from the principle of superposition of solutions for linear equations, it is sufficient to consider the evolution of a single mode. To assess whether solutions are damped or amplified, we thus look for solutions φ z (for z = n, u, v) of the form 1 for C z some constant and λ a complex number whose real part is the rate at which the perturbation is damped (real part of λ < 0) or amplified (real part of λ > 0). We note Φ = (φ n , φ u , φ v ) and C = (C n , C u , C v ). In this equation, we observed that if a mode is amplified on one coordinate, it is also on the other, and we thus make the simplifying assumption k n 1 = k u 1 = k v 1 and k n 2 = k u 2 = k v 2 . The linearized equation for this perturbation reads: implying, when inserting the specific form of Φ in the equation, that: In other words, λ is an eigenvalue of the matrix on the righthand side.

Stability of the trivial steady state
For the trivial spatially homogeneous state (0, 0, 0), the matrix associated with stability is given by: so the expressions of eigenvalues are straightforward: Because µ k 1 ,k 2 is non-positive, λ 2 , λ 3 are strictly negative, but because α n > 0, λ 1 is positive for k 1 = k 2 = 0. This implies that the trivial state is always unstable: a small cell density will progressively grow because of the logistic growth term.

Stability of the non-trivial equilibrium
To determine the values of λ in the case of the non-trivial equilibrium, we use the classical fact that eigenvalues are the roots of the polynomial: = λ 3 + a 2 (k 1 , k 2 )λ 2 + a 1 (k 1 , k 2 )λ + a 0 (k 1 , k 2 ) where The Routh-Hurwitz criterion provides three conditions for ensuring that eigenvalues have strictly negative real part (i.e. that the equilibrium is stable). In detail, all eigenvalues have strictly negative real part and only if we have simultaneously a 0 , a 2 ≥ 0 and a 2 a 1 > a 0 . An instability thus arises when one of these conditions is not satisfied, i.e.: a 0 (k 1 , k 2 ) < 0 or While these expressions may look complex, they are relatively simple functions of µ k 1 ,k 2 . We found that for a variety of parameters in the vicinity of the central parameter chosen, only the term a 0 can break the condition of stability of Routh-Hurwitz.
We plotted in S4 and S6-9 Figs the corresponding dispersion relations showing, for the reference (red), highest (orange) and lowest (blue) displayed parameter sets, the evolution of the real part of largest eigenvalues of the stability matrix (10) as a function of the wavenumbers k (in these plots, we took k = k 1 = k 2 , allowing a simple visualisation). These curves highlight the existence of Turing instabilities in cases where a finite range of positive λ exist: the associated modes become unstable and are amplified, yielding the emergence of a pattern.
We note that the non-trivial steady state is independent of the diffusion and chemotaxis coefficients. We thus considered its stability as D u and D v (respectively, D u and κ) are varied in combination, as is shown in S6 Fig (respectively, S7 Fig). We found that instabilities require to have D u relatively small and D v relatively large (respectively, D u relatively small and κ relatively large). This is perfectly consistent with the loss of patterning when varying these parameters, as observed in S6-7 Fig.