Breaking constraint of mammalian axial formulae

The vertebral column of individual mammalian species often exhibits remarkable robustness in the number and identity of vertebral elements that form (known as axial formulae). The genetic mechanism(s) underlying this constraint however remain ill-defined. Here, we reveal the interplay of three regulatory pathways (Gdf11, miR-196 and Retinoic acid) is essential in constraining total vertebral number and regional axial identity in the mouse, from cervical through to tail vertebrae. All three pathways have differing control over Hox cluster expression, with heterochronic and quantitative changes found to parallel changes in axial identity. However, our work reveals an additional role for Hox genes in supporting axial elongation within the tail region, providing important support for an emerging view that mammalian Hox function is not limited to imparting positional identity as the mammalian body plan is laid down. More broadly, this work provides a molecular framework to interrogate mechanisms of evolutionary change and congenital anomalies of the vertebral column. https://doi.org/10.1038/s41467-021-27335-z OPEN

T he mammalian vertebral column is comprised of seriallyrepeating vertebrae that, while individually-identifiable, are grouped based on morphological and functional similarity into cervical (C), thoracic (T), lumbar (L), sacral (S) and caudal (Ca) elements. Vertebrae arise from developmentally-transient structures called somites, that are in turn generated from a progenitor pool located in the posterior embryo through the coordinated processes of axial elongation and segmentation 1 . The genetic mechanisms dictating vertebral identity are well understood, and centre on the temporally-controlled activation of Hox cluster genes 2 within vertebral progenitors 3 . In contrast, uncovering the genetic mechanisms that control vertebral number of each axial region in a meristic rather than homeotic manner (as defined in Bateson, 1894 4 ), and consequently total vertebral number (TVN) of a mammalian species, has proved incredibly challenging.
The changes that constitute diversity of mammalian axial formulae are not uniform along the anterior-posterior (A-P; head-to-tail) axis 5,6 . Cervical number is almost exclusively fixed at 7, a trait that has rarely changed over 200 million years of mammalian evolution, spanning back to the mammalian ancestor of the late Triassic 7 . Two rare exceptions to this rule, manatee and sloth, have sparked intense interest in the advantage and developmental basis of these break in constraint [8][9][10] , with their genetic basis unresolved to date. Thoraco-lumbar (T-L) count, while less rigid, is still highly conserved at 19 or 20 across most mammals, though again with examples of relaxation in both absolute T-L number and intraspecies robustness seen in clades such as Xenarthra and Afrotheria 11,12 . Sacral count varies considerably across mammals and along with the extreme diversity seen in caudal/tail morphology reinforces a graded decrease in constraint along the A-P axis.
As mouse displays the likely ancestral mammalian presacral formulae 6,13 , and Rodentia exhibits the lowest deviation from median vertebral counts of all mammals 12 , it should serve as an excellent model with which to dissect genetic mechanisms constraining axial formulae. However to date, although many truncating mouse mutants exist, often being highly dysmorphic, the ability to increase TVN in a meristic not simply homeotic manner has been limited to perturbation of only two genetic pathways [14][15][16] (see Supp Information for note on Hox13 mouse mutant phenotypes). The largest of these two meristic increases was observed in, perhaps intuitively, the evolutionary more plastic tail region; ectopic Lin-28 was found to increase the number of tail elements by 5 while let-7 opposed this function 14,15 . Within the less variable presacral region, microRNA(miR)-196 activity has been shown to constrain thoracic number by two, constituting one homeotic and one meristic change 16 . The signalling molecule Gdf11 has a far greater role than miR-196 in constraining T-L count as seen in conditional or complete knockout mice 17,18 , and partial loss of Cdx transcription factors can increase presacral number by 1 19 , but each does so at the expense of the caudal skeleton. Whether more subtle manipulations such as timing of signal onset 20 or exact expression levels provide a genetic basis for divergent mammalian morphologies, and more broadly, the degree to which genetic compensation or pleiotropic effects 21,22 underly evolutionary constraint of axial formulae, remains unresolved.
Here, we show that in vitro modelling of axial progenitors identifies unexpected cooperation between transcriptional and post-transcriptional regulatory mechanisms that act to terminate presacral expression signatures. Applying this knowledge in vivo, we are able to experimentally prolong the natural process of axial elongation in the mouse and considerably increase total vertebral number within both presacral and caudal regions. We reveal a high level of redundancy and synergism between Gdf11, miR-196 and retinoic acid (RA) in the constraint of axial formulae, with Gdf11 and RA cooperativity impacting a surprisingly large extent of the vertebral column including the constraint of cervical number to 7. All regulatory mechanisms converge in their ability to control spatio-temporal Hox expression, a parameter that our data supports having a substantive role in axial elongation, at least in the tail region.

Results
Synergistic constraint of trunk Hox expression signatures. To dissect mechanisms constraining regionally-restricted Hox expression, as a proxy for regional axial morphology, we initially focused on modelling the trunk-to-tail (T-to-T) transition in vitro since progression through this critical moment in body axis formation correlates with the switch from an expanding to a depleting and eventually exhausted progenitor pool in species as diverse as mouse and snake 23 . miR-196 has a non-redundant role in timing the T-to-T transition, and thus we generated induced pluripotent stem cells (iPSCs) from wild-type (WT) and miR-196a1 gfp/gfp ;miR-196a2 gfp/gfp ;miR-196b −/− (miR-196-TKO) isogenic mice using standard protocols 24 (Supp. Figure 1a). Both iPSC lines were differentiated to model the developmental kinetics of a posterior growth zone by sequential addition of FGF2 on Day (D) 0 and the Wnt pathway agonist CHIR99021 on D2 (Fig. 1a, based previous work [25][26][27][28]. miR-196-TKO iPSCs were able to generate axial progenitors normally, however, exhibited marked differences in the kinetics of Hox cluster progression when compared to WT (Fig. 1b). Trunk Hox genes Hoxb8 and Hoxc8 displayed an identical timing of activation for both genotypes, but as differentiation proceeded, their quantitative level and temporal persistence were significantly increased in miR-196-TKO cells (Fig. 1b), consistent with the presence of functional miR-196 binding site(s) within their 3′UTRs 16,29 . Addition of Gdf11 to culture conditions was able to suppress trunk Hox gene expression and induce posterior Hox activation in both genotypes as anticipated. However, in the absence of Gdf11, both Hoxb8 and Hoxc8 displayed an unrestrained upwards trajectory in miR-196-TKO cells but not in WT (Fig. 1b). This key result indicated that miR-196 and Gdf11 are not simply redundant in this context, but act in concert in shutting down a trunk Hox code, a parameter suggested to support axial elongation under certain circumstances 30 .
Gdf11 and miR-196 synergistically constrain TVN. To determine whether the combined function of Gdf11 and miR-196 identified in vitro influences axial elongation in vivo, we interbred miR-196-TKO 16 and Gdf11 −/− 17 mice and characterised TVN and vertebral identity across the allelic deletion series (Fig. 2a- Table 1). Consistent with previous reports, complete loss of miR-196 alone resulted in an increase of 2 thoracic elements (C7, T15, L5, median TVN 62), while complete loss of Gdf11 alone resulted in a dramatically expanded presacral region (C7, T18, L8/9) that truncated soon after sacral elements formed. In contrast to the Gdf11 −/− truncation phenotype, we found that Gdf11 +/− embryos presented with an additional 1 thoracic and 2 caudal vertebrae (C7, T14, L6, median TVN 64), a surprising dose-dependent effect of Gdf11 on TVN that has not previously been appreciated. Crossing miR-196-TKO onto this Gdf11 +/− background increased the number of thoracic vertebrae further, producing viable offspring with a total of 5 additional elements (C7, T16, L6, median TVN 66). Gdf11 −/− ;miR-196-TKO compound mutant embryos demonstrated a remarkable increase in the number of presacral elements with an additional 9 thoracic and 4 lumbar elements compared to WT (C7, T22, L10; Fig. 2b), the greatest expansion of presacral number reported in mice to date. Together, these data provide 3 major advances: (i) Gdf11 activity is required under wild-type conditions to constrain total vertebral number; (ii) The vertebral column is exquisitely sensitive to dosage of Gdf11, with opposing outcomes revealed dependent on allele number; (iii) miR-196 and Gdf11 do not constitute entirely parallel mechanisms in the control of presacral vertebrae number but rather, act synergistically, since the effect of their combined loss is greater than the sum of individual mutant phenotypes (Fig. 2d).
Retinoic acid has pleiotropic effects on axial formulae. While TVN cannot be accurately quantified in Gdf11 −/− embryos due to the highly dysmorphic nature of post-sacral elements, the loss of miR-196 activity did not qualitatively appear to rescue Gdf11 −/− truncation (Fig. 2a). In this context, the partial rescue of Gdf11 −/− truncation by reducing endogenous levels of Retinoic acid (RA) signalling 31 led us to examine the potential for more elaborate redundancy or synergy between the three signalling/regulatory pathways. Indeed, oral gavage of pregnant dams across the genetic deletion series with pan-RA receptor inhibitor AGN193109 (AGN) revealed further combinatorial changes in vertebral number and identity, spanning each major subdivision of the vertebral column  Table 2). Notably, we showed that RA constrains TVN in mouse, with an additional 1-2 elements observed following RA depletion for all genotypes where TVN can be quantified, including WT, Gdf11 +/− , miR-196-TKO and Gdf11 +/− ;miR-196-TKO embryos. This increase manifested as an expansion in the number of presacral elements, the identity of which (see 'Methods') was dependent on the exact allelic combination (Compare Supp. Figure 2a with Fig. 2a). AGN-exposed WT embryos yielded anteriorising transformations of the cervical region (C2 → C1, C7 → C6) at low penetrance, no change in positioning of cervico-thoracic transition (Fig. 3c) and one additional T element (Supp. Figure 2a). All phenotypes were consistent with observed phenotypes in Retinoic acid receptor γ (RARγ) knockout mice 32 . The deletion of even a single Gdf11 allele under these RA-depletion conditions lead to a higher penetrance of C7 → C6 transformation, and in 40% of embryos, serial transformation resulted in a shift in cervico-thoracic positioning such that 8 cervical elements formed (Fig. 3c). This striking shift in cervico-thoracic positioning, widely considered as one of the strongest evolutionarily constrained traits of mammals 6,33 , became almost fully penetrant in Gdf11 −/− embryos treated with AGN, and where present, usually formed at the expense of a thoracic element and thus was homoetic in nature. Vertebral identity at the cervico-thoracic junction is known to involve the action of Hox5/6 paralogous groups 34 early in development. Consistent with this, and with observed phenotypic changes, we found the anterior boundary of Hoxc6 was shifted caudally by 1 somite in E10.5 Gdf11 −/− mutants treated with AGN (Supp. Figure 2b), demonstrating redundancy between RA and Gdf11 in the timely activation/spatial regulation of Hoxc6, and likely other trunk Hox genes. Countering this displacement of the cervico-thoracic junction, deletion of miR-196 paralogs returned this major morphological transition back to normal (Fig. 3d, Supp. Table 2 Table 2). Collectively, this extensive genetic/chemical deletion series has revealed that RA, Gdf11 and miR-196 all act to constrain TVN individually, and synergistically. In addition, each factor differentially shapes positional identity along the A-P axis, with individual, synergistic and, in this case, antagonistic interactions revealed. It is important to note that in these embryos, and all embryos assessed in this study, progression through the T-to-T transition always occurred, albeit late. This indicates that the regulatory synergism promoting this key transition was yet to be fully depleted, with prime candidates supporting the eventual T-to-T in compound mutant embryos being Gdf8 36 and potentially FGF signalling 37 .
Factors constraining TVN differentially shape Hox codes. One shared feature of the extrinsic and intrinsic mechanisms shown here to constrain TVN and shape regional identity is their capacity to influence spatio-temporal Hox expression 16,17,38 . Thus, to reveal the full extent to which miR-196 and Gdf11 individually and collectively regulate Hox codes during the T-to-T transition (somites 21 +/− 1 [~E9.5]) or after this transition is complete (somite 32 +/− 1 [~E10.5]), we quantified expression of all 39 Hox genes within tailbud tissue across the allelic deletion series (Fig. 4). Removal of miR-196 repressive activity (miR-196a2 −/− ;b −/− , phenotypically equivalent to miR-196-TKO) led to a robust upregulation of trunk Hox genes and concomitant downregulation, or failure to timely activate, posterior Hox genes  of all four clusters at E9.5 ( Fig. 4) 16 , a signature that began to resolve by E10.5 (Fig. 4). From a patterning perspective, either an increased trunk Hox code 39 or a delay in activation of Hox10 ribsuppression activity 3,40 have the potential to drive the expanded thoracic identity seen in miR-196-TKO embryos. In contrast, loss of one Gdf11 allele had minimal impact on trunk Hox expression, indicating that the mildly expanded thoracic phenotype of Gdf11 +/− embryos is more likely to result from the modest heterochronic delay in posterior Hox code activation seen at E9.5, which largely resolved to WT levels by E10.5 (Fig. 4). This is further supported in Gdf11 −/− embryos, where the level of trunk Hox upregulation was not consistent with the dramatically expanded thoracic phenotype, however, a striking reduction in expression of all posterior Hox genes was seen at E9.5 (Fig. 4),

+5
- Gain of vertebra within axial region Loss of vertebra within axial region  Table 1). Collectively, the altered Hox codes observed across the allelic deletion series allow us to rationalise the molecular basis for observed patterning changes ( Fig. 2d; Supp. Table 1). Moreover, the ability for ectopically-expressed trunk Hox genes Hoxa5 or Hoxb8 to drive axial elongation under certain circumstances 30 and the suggested importance of timely activation of Hox13 paralogs in terminating mouse axial elongation 15,30 , raised the question as to whether the heterochronic and quantitative changes in global Hox signatures observed here may also underlie changes in TVN.
Posterior Hox expression supports tail formation. The understanding of Hox12 and Hox13 paralog function during vertebral column formation is incomplete at best due to the current lack of paralog group mouse knockouts, and recently, the view of Hox13   Fig. 3 Manipulation of multiple regulatory mechanisms is required to further break constraint of mouse axial formulae. a-d In vivo characterisation of miR-196 and Gdf11 individual and compound mouse mutant skeletal phenotypes following treatment with RA-receptor inhibitor AGN193109. a Representative embryonic day (E)18.5 skeletal preparations across genotypes. C = cervical; T = thoracic; L = lumbar; total vertebral number indicated in brackets. b Quantification of total vertebral number across genotype. Raw data is presented in the upper plot (vertical error bar = mean and standard deviation). Mean differences against the respective untreated control are presented in the lower plot as bootstrap sampling distributions. Mean differences are depicted as dots and 95% confidence intervals are indicated by the ends of the vertical error bars. n refers to the number of individual animals used for this analysis. Source data are provided as a Source data file. c Altered skeletal identity surrounding the cervico-thoracic boundary in AGN193109-treated WT and Gdf11 mutant embryos. *=anterior arch of the atlas; TA = tuberculum anterior; yellow arrow indicates a rib anlage on the 8th vertebral element. d Schematic summary of vertebral alterations observed across the AGN193109-treated Gdf11;TKO allelic deletion series, relative to WT. Numbers represent the maximum phenotypic difference observed from the untreated wild-type control. Diff = difference, S-Ca = sacral-caudal vertebrae, question marks indicate dysmorphic and non-quantifiable elements. White drawn line = partially penetrant, often incomplete or unilateral, alterations observed at the first thoracic element (Supp. genes as endogenous axis terminators has been challenged in fish 41 . To interrogate posterior Hox function, and specifically, to dissect whether the delay in/loss of posterior Hox activation seen in Gdf11 +/− and Gdf11 −/− mutant embryos is of any phenotypic consequence, we sought to restore timely posterior Hox expression through the creation of transgenic mice expressing either Hoxd11 (Hoxd11 OE ) or Hoxd12 (Hoxd12 OE ) under the control of Cdx2 regulatory elements 42 (Supp. Figure 4a- Initial characterisation of these lines on an otherwise WT background revealed a reduction in lumbar count by 1 element and concomitant reduction in TVN, with both qualitative (Hoxd12 > Hoxd11) and quantitative (copy number) enhancement of phenotype observed (Fig. 5a- Table 3). This very mild reduction phenotype of the presacral region was surprising when compared with full tail truncation observed following ectopic expression of Hox13 paralogs 15,30 , which may stem from quantitative differences between transgenics or likely functional differences between these Hox proteins 43 . A near-identical reduction in lumbar count/TVN was observed when either Hoxd11 OE or Hoxd12 OE was bred onto a Gdf11 +/− background (  Table 3), supporting the view that the temporary delay in posterior Hox activation in these heterozygous embryos (Fig. 4) may indeed contribute to a shifted T-to-T transition and elongation phenotype, or at a minimum, that ectopic posterior Hox expression is dominant over elongation mechanism(s).
Each of the above genetic crosses represents a cumulative Hox scenario that is ectopic over WT levels, and thus we next cross-bred either Hoxd11 OE or Hoxd12 OE onto the Gdf11 −/− background, which is greatly depleted for posterior Hox expression (Fig. 4). Relative to truncated Gdf11 −/− embryos, restoration of a single posterior Hox gene yielded striking restoration of a segmented Uncx4.1 + tail-like structure at E12.5 and E13.5, albeit ventrally displaced ( Fig. 6a; Supp. Figure 5a-b). Detailed characterisation of Gdf11 −/− embryos via tissue sections (Supp. Figure 6a) revealed in many cases that the notochord underlying the primary neural tube had turned abnormally relative to WT, projecting ventrally at the level of the cloaca. There, notochordal cells were seen interspersed with Sox2 + -neural and enveloped by Foxa2 + -endodermal cells, with no observable patterning or structure to this previously noted ventral mass 44,45 (Supp. Figure 6a [v-vii]). The restoration of either Hoxd11 or Hoxd12 expression in Gdf11 −/− embryos did not prevent ventral projection of the notochord, however, these cells no longer became trapped and extended to the tip of the ventral tail in an organised manner (Fig. 6b, Supp. Figure 6a [ix-xii], [Supplementary movies [1][2][3]). Adjacent to the notochord in these Gdf11 −/− ;Hoxd11/12 OE embryos, an organised Foxa2 + -tailgut formed, extending from the cloacal orifice to the tip of the ventral tail, with distal cells coexpressing Foxa2 and Sox2 (Fig. 6b, Supp. Figure 6b, [Supplementary movies 1-3]) 46 . One signature component of a wild-type tail that could not be clearly delineated in this ventral tail was a Sox2 + secondary neural tube, despite the tailbud housing all major progenitors required for posterior body formation based on marker analysis (Fig. 6c, Supp. Figure 6b). Collectively, the presence of a patterned tail structure in Gdf11 −/− ;Hoxd11/12 OE embryos, and the maintenance of axial progenitors in this ventrally-displaced tailbud (Fig. 6c) at stages when these progenitors should normally be close to exhaustion 47 , strongly supports a minimum requirement for posterior Hox expression as an integral component of the gene regulatory network necessary to construct the tail. Despite this striking result at mid-gestation, skeletal analysis of E18.5 Gdf11 −/− ;Hoxd11 OE or Gdf11 −/− ;Hoxd12 OE embryos revealed no vertebrae arising from ventrally-displaced Uncx4.1 + -somites (Supp. Figure 5b). In fact in many Gdf11 −/− ;Hoxd12 OE embryos, the last-formed vertebra was found more rostrally than in Gdf11 −/− embryos alone (Supp. Figure 5c), likely due to the prolonged expression of posterior Hox genes resulting in increased apoptosis, decreased proliferation and degeneration of more mature tail structures 48 and/or a lack of secondary neurulation. In this light, these mutant embryos were more akin to human development, where a patterned tail is initially generated but fully regresses over time 49 .

Discussion
Historical hereditary studies have long supported a multigenic contribution to even mild examples of intraspecies vertebral variation 50 . Here, we reveal the highly integrated manner by which multiple regulatory layers constrain the modular logic of the vertebral column. Importantly, we were able to produce viable offspring with 5 additional vertebral elements spanning both trunk and tail regions, and with a maximal expansion of 7 additional elements. Regarding the latter, we expect that spatial and/or temporal conditional deletion of these regulatory mechanisms would circumvent embryonic lethality, enabling a vastly altered mammalian body plan, though the potential secondary consequences of this manipulation on locomotion or internal organ formation cannot be predicted. Combined with the recent identification of the Lin28/let-7 axis in constraining vertebral number in the tail 14,15 , these results demonstrate the depth of genetic redundancy that normally acts to constrain each vertebral region in the mouse and provide robustness to the system as a whole.
The comprehensive allelic deletion series performed allowed us to investigate how the breaking of regional constraint within this experimental model aligns with current hypotheses of natural variation. For example, the break in C7 constraint observed in sloths has been suggested to result from a mis-alignment of somite-derived (primaxial) vertebral elements with that of lateral plate mesoderm-derived (abaxial) distal rib/sternum, such that 7

total vertebral number (TVN) in WT and
Gdf11 +/− embryos, with or without the presence of Hoxd11 OE and Hoxd12 OE transgenes. Raw data is presented in the upper plot (vertical error bar = mean and standard deviation). Mean differences relative to WT are presented in the lower plot as bootstrap sampling distributions. Each mean difference is depicted as a dot and 95% confidence interval is indicated by the ends of the vertical error bar. n refers to the number of individual animals used for this analysis. Source data for b are provided as a Source data file. c Schematic summary of changes in axial formulae, relative to WT, observed in single Hoxd11 OE or Hoxd12 OE transgenic lines, and following cross-breeding of each transgenic with the Gdf11 mutant line. Numbers represent the rounded (to the next whole number) unpaired mean difference for a given genotype. C = cervical, T = thoracic, L = lumbar, S-Ca = sacral-caudal, question marks indicate dysmorphic and non-quantifiable elements. White drawn line = frequently observed reduction/malformation of ribs on the last rib-bearing thoracic element in mice carrying the Hoxd11 OE or Hoxd12 OE transgene. cervical vertebral bodies always form and axial diversity arises due to the gain or loss of ribs 8,9 . Here, in RA-depleted Gdf11 −/− embryos with 8 cervical elements, we observe serial anteriorising homeotic transformations for at least 2 elements on either side of the C-T boundary supporting a genuine increase in cervical number. However, in the majority of RA-depleted Gdf11 +/− embryos where 7 cervical elements formed, we see that the pattern of centra ossification of the first rib-bearing element is more indicative of a cervical element, and T2 is displaced posteriorly. This suggests that loss of even one Gdf11 allele in this context has the ability to preferentially repattern primaxial tissue, and supports the mis-alignment of primaxial and abaxial tissues as a mechanism to maintain constraint (here) or drive diversity (in sloths).
At a molecular level, the 3′-to-5′ sequential activation of Hox cluster expression within axial progenitors over development time -the Hox clock 51 -prefigures Hox spatial collinearity along the A-P axis 52,53 , though the consequences of clock manipulation are not easy to address. The mouse mutants presented here, all unified in their temporal control of global Hox transitions, offer a granular view as to the consequences of manipulating that Hox clock in vivo. We show that the speeding up, or the slowing down, of Hox cluster signatures result for the most part in serial homeotic transformations that associate with changes in total vertebral number. At a minimum, this implies a very tight association between Hox patterning and elongation mechanisms, though importantly, our data also revealed that posterior Hox genes have the capacity to positively influence axial elongation. This latter data, viewed along with similar results for trunk Hox genes in mouse 30 and Hox13 paralogs in zebrafish 41 , allows us to propose a model whereby multiple post-occipital expressing Hox paralogs participate in construction of the main body axis, not solely in its patterning. We do not suggest that Hox genes are the primary drivers of axial elongation, but that a minimum level is required. This model does not preclude a role for high levels of posterior Hox expression in terminating axial elongation 15,30,43 , and indeed our data on ectopic Hoxd11/Hoxd12 can be viewed as supporting this. However, caution should be taken when interpreting overexpression studies 41 while awaiting the genetic deletion of Hox12 and Hox13 paralog groups in the mouse. Why then has a role for Hox genes in shaping vertebral number not been apparent in the extensive Hox mouse mutant literature? Cumulative mutant analysis has shown that each individual vertebral element is patterned by at least two, if not more, Hox paralog groups 34 , thus only by removal of multiple adjacent paralog groups would this function be revealed. Experimentally, this would be a complex undertaking and of questionable relevance from an evolutionary perspective. However, by altering the timing of collective Hox code transitions (i.e. in vivo pacing of the Hox clock) rather than Hox function per se, through changes in widespread post-transcriptional regulation (miR-196) or alterations in higher-order signalling (Gdf11 haploinsufficiency), this unanticipated vertebrate Hox function has now been revealed.
From a cellular perspective, we show here the persistence of axial progenitors in the tailbud of Gdf11 −/− ;Hoxd12 OE embryos beyond what is seen in WT. Whether the restoration of posterior Hox expression in this scenario is acting cell-autonomously to maintain progenitor proliferation and support tail formation remains to be elucidated, though recent work in zebrafish argues against such a role 41 . Serial transplantation of axial progenitors in vivo 54 has shown that these cells do not intrinsically exhaust following a timer mechanism, highlighting the importance of extrinsic signals in the maintenance and exhaustion of a progenitor pool. In this light, the restoration of tailgut formation along the full extent of the Gdf11 −/− ;Hoxd12 OE ventral tail was of particular interest, since surgical removal of the caudal endoderm in chick has been shown to redirect the notochord ventrally into the hindgut 55 with striking similarity to what is seen in Gdf11 −/− mutants. Moreover, expression of a dominant-negative Hoxa13 protein within chick caudal endoderm, but not mesoderm, resulted in tail truncation 55 , further supporting the positive requirement for posterior Hox function in tail construction across multiple species, and highlighting important cellular targets for future investigation.
By focusing on how one mammalian species cannalises axial formulae, this work has reshaped the current understanding of body plan formation and the critical function of Hox networks within. This should help inform the basis for supernumerary vertebral elements that have been reported to appear in~3% of the human population 56 , constituting congenital anomalies such as an 8th cervical vertebra 57 , additions of up to 3 thoraco-lumbar vertebrae [58][59][60] and rare cases of humans with tail vertebrae 61 . Moreover, the comprehensive allelic deletion series performed here provides us with a window into how evolutionary changes may have occurred, not by complete ablation of one signal or regulatory mechanism, but by more subtle and/or tissuerestricted changes in multiple pathways that converge in their ability to shape axial formulae. Mice. miR-196a1 GFP , miR-196a2 GFP , miR-196a1 −/− , miR-196a2 −/− and miR-196b −/− mouse lines have been previously described 16 , and were maintained on a C57BL/6 (mixed substrain J/N) background. Gdf11 −/− mutant mice have been previously described 17 and were maintained on a C57BL/6J background. Cdx2P:Hoxd11 (Hoxd11 OE ) and Cdx2P:Hoxd12 (Hoxd12 OE ) transgenic mice were generated in-house by pronuclear injection into C57BL/6J zygotes according to standard protocols 62 . cDNA for mouse Hoxd11 or Hoxd12 was inserted downstream of the 9.4 kb Cdx2 regulatory region (Cdx2P) 42 . Founder animals were verified by standard PCR, copy number assessed by digit drop PCR and germline transmission confirmed (Cdx2P:Hoxd11 1 line established; Cdx2P:Hoxd12 2 lines established).

Methods
Oral gavage. The pan retinoic acid receptor inhibitor AGN193109 (Tocris) was administered by oral gavage to pregnant dams at a concentration of 0.8 mg/kg on 3 consecutive days of development (E7.5, E8.5 and E9.5). Stock AGN193109 was resuspended in DMSO (2 mg/mL), stored at −20°C and freshly diluted in corn oil (1:25) on the first day of administration.
Skeletal preparation and imaging. Skeletal preparation was performed on E18.5 embryos or p0 postnatal pups as previously described 63 . Embryos or pups were skinned and eviscerated, followed by sequential 2-day incubations in 100% ethanol, 100% acetone and stain solution (containing 0.15 mg/ml Alican blue, 0.5 mg/ml alizarin red, 70% ethanol and 5% glacial acetic acid), with constant rocking at room temperature (RT). Soft tissue matter was cleared from the skeleton by incubation in 1% potassium hydroxide (KOH) for 2-5 days. Once sufficiently cleared, samples were subjected to an increasing glycerol/1% KOH gradient (25%, 50%, 75%, 100% glycerol) for a minimum of one day each. Samples were stored in 100% glycerol. Skeletal images were taken with the Vision Dynamic BK Lab System at Monash University in the Paleontology Lab on a Canon 5d MkII with a 100 mm Macro lens (focus stop 1:3/1:1). To extend the focal depth on the skeletal images multiple images were taken and stacked in ZereneStacker using the PMax algorithm.
Skeletal phenotyping and statistical analysis. Minor differences in axial formulae have been reported between mouse strains, and thus it was of critical importance that all analyses were performed on an isogenic C57BL/6 background, decreasing inter-animal variation and allowing robust quantification. For each skeletal preparation, axial formulae were assessed independently by two individuals, blinded to genotype. Where unilateral identity changes were observed, the vertebral identity that deviates from WT was annotated, e.g. unilateral L1 rib would be annotated as a L1-to-T transformation. In skeletons that have been exposed to AGN193109, we commonly observed a dorsal bifurcation of C2 which fused laterally into one and therefore were counted as one vertebral element. Unpaired mean differences in total vertebral number counts are presented in Supplementary  Table 4. Animal sex has not previoulsy been shown to alter skeletal formulae and thus was not assessed. For statistical analysis and data visualisation of vertebral numbers, the R-package "Data Analysis using Bootstrap-Coupled ESTimation" (DABEST) 64 was used. To determine mean differences to the respective shared wild-type control 5000 bootstrap samples were taken and the confidence interval was bias-corrected and accelerated.
Droplet digital PCR assay and copy number detection. Each droplet digital (dd) PCR assay contained 20 ng of DNA template, 1×ddPCR Master Mix (Bio-Rad, USA), 0.5 μM primer-F, 0.5 μM primer-R and 0.25 μM probes. Primer and probe sequences are shown in Supplementary Information. A QX200 ddPCR droplet generator (Bio-Rad, USA) was used to divide each 20 μL ddPCR mixture into approximately 20,000 droplets. The ddPCR probe for the reference sequence (Rpp30) was labelled with hexachlorofluorescein (HEX) at the 5′ end and black hole quencher (ZEN/IowaBlack) at the 3′ end. The ddPCR probe for the target sequence (Hoxd12) was labelled with 6-carboxyfluorescein (FAM) at the 5′ end and black hole quencher (ZEN/IowaBlack) at the 3′ end. The thermal parameters were as follows: 10 min at 95˚C, followed by 40 cycles of 30 s at 94˚C and 1 min at 60˚C, followed by 98˚C for 10 min. The amplified products were analysed using a QX200 droplet reader (Bio-Rad, USA) and copy number analysed using QuantaSoft Analysis Pro Version 1.0 (Bio-Rad, USA).
In situ hybridisation. Whole-mount in situ hybridisation was performed as previously described 65 with some modifications. Embryos were dissected in ice-cold PBS (Gibco, 14190-144) and tissue overlying the hindbrain pierced. Embryos were fixed in 4% paraformaldehyde (PFA), rocking overnight at 4°C. All steps were performed with gentle rocking at RT unless otherwise specified. Embryos were washed twice with PBS-Tween (0.1%, PBT) for 5 min. Embryos were dehydrated using a methanol/PBT gradient (25%, 50%, 75%, 100% methanol) with each step ranging from 5-20 min depending on age. Embryos could be stored in 100% methanol, or to begin in situ hybridisation, embryos were rehydrated following a reversed methanol/PBT series as above, washed twice with PBT for 5 min and treated with 10 µg/ml of proteinase K in PBT (E9.  Hoxb13).
Cell lines. The Bruce4 mouse embryonic stem cell line was originally isolated from a C57BL/6 mouse 66 and kindly provided by Dr Jeff Mann. Two in-house generated iPSC lines were derived from the miR-196 triple knockout mouse strain 16 and wild-type mice of an isogenic background.
Generation of mouse induced pluripotent stem cells (iPSCs). iPSCs were generated from tail tip fibroblasts of adult WT and miR196a1 gfp/gfp ,a2 gfp/gfp ,b −/− mice as previously described 24  Gene expression analysis by quantitative PCR. Additional quantitative PCR analysis was performed on a Lightcycler 480 (Roche) using the SYBR Green I Master Mix (Roche). Two microlitres of cDNA was amplified per reaction using the following programme: 95°C for 10 s (1 cycle), 95°C for 10 s, 60°C for 15 s, 72°C for 10 s (45 cycles). Primers used in these studies were either previously published 69 or synthesised in-house and listed in Supplementary Table 6.
Dissection and tissue processing for sectioning. Whole mouse embryos from various stages were dissected and fixed overnight in 4% PFA/PBS at 4°C. Samples were washed twice for 10 min in PBS at 4°C. Samples were equilibrated through a series of sucrose solutions (5%, 20%, 30% sucrose in PBS) before being embedded in OCT medium (TissueTek) using an ethanol-bath and stored at −80°C. Frozen tissue was sectioned at 12-14-μm thickness on a Leica cryotome. Sections were mounted directly onto slides (Superfrost, Fisher Scientific).
Imaging and image processing. Whole-mount brightfield images of E12.5 and E13.5 embryos were acquired using a NSZ-405 Zoom Stereo Microscope at ×1.5 and ×4 magnifications. Images were adjusted for brightness and some images were reflected to have the tails oriented in the same direction for easy visualisation. Fluorescent images of slide-mounted embryonic sections were imaged with a Zeiss AxioImager Z1 with an AxioCam HRm camera at a ×20 objective magnification. All images were processed using AxioVisionRel.4.8 and Fiji softwares. Sectioned embryonic images of Gdf11 −/− with or without Cdx2P:Hoxd12 tail regions were automatically tiled on the Zeiss AxioImager Z1 and compiled into one image using AxioVisionRel.4.8 in-built image stitching function MosaiX. Tiled image and individual image size ratios were maintained to ensure direct comparisons between tiled images. Single-cell layer images of E12.5 split tail tips were acquired using Zeiss AxioImager Z1 in Apotome mode (Gitter 7.5 LP/mm gef; 423662-0000-000) and ×20 objective magnification. Linear adjustments, such as for "Brightness" or "Contrast", in which the same change is made to each pixel according to a linear function, were used on whole immunofluorescence images using Fiji (ImageJ2.0) or AxioVision Software. Fluorescence images of individual serial sections were uniformly processed using identical linear adjustments. For the cropped regions of apotome acquired split tail regions linear adjustments were made to all images and background subtraction was applied to the single-fluorescence images to extract/distinguish the signal above autofluorescence.
Images of E18.5 skeletal preparations and E10.5 embryos following wholemount RNA in situ hybridisations were acquired with a Vision Dynamic BK Lab System at the Monash University Paleontology Lab. Images were taken with a Canon 5d MkII with a 100 mm Macro lens (focus stop 1:3/1:1). Multiple images were taken to extend the focal depth and stacked in ZereneStacker using the PMax algorithm. Scales were determined via known pixel counts for each focus stop in Adobe Photoshop.
For whole-mount immunofluorescence-stained embryos, Z sections were acquired every 2 µm on a lightsheet microscope (Ultramicroscope 2, LaVision, Miltenyi Biotec GMbH). Image stacks were imported into Imaris (v9.6; BitPlane, Oxford Instrument) and visualised in 3D. Contrast was adjusted for each channel and segmentation was performed using the Surface tool. Pictures and animations were created in Imaris.
Adobe Illustrator was used to compose all multipanel figures in this study.
Statistics and reproducibility. All skeletal phenotyping was performed independently by two investigators, blinded to genotype. No statistical method was used to predetermine sample size. For gene expression analyses, Ct values greater than the manufacturer recommended threshold were excluded.