Axial skeleton anterior-posterior patterning is regulated through feedback regulation between Meis transcription factors and retinoic acid

Anterior-posterior (AP) patterning is an essential feature of the bilaterian body plan and its mechanisms have been extensively studied. Segmental epithelial sacs known as somites emerge from the paraxial mesoderm as it is produced and these sacs progressively incorporate to the AP axis. The initially homogeneous somites later subdivide into compartments, including the sclerotome, which is the precursor of the vertebrae and ribs, and the myotome, which is the precursor of the skeletal muscles (Musumeci et al., 2015). Crosstalk from the myotome to the sclerotome is essential for sclerotome patterning and in particular for rib specification and patterning (Vinagre et al., 2010).

An important breakthrough in understanding AP axis patterning was the identification of Hox mutants in Drosophila, which cause the transformation of one part of the body into another: a phenomenon known as homeotic transformation (Lewis, 1978). Hox genes are conserved in evolution and appear to be organized in genetic complexes in most animals (Duboule and Dolle, 1989; Sanchez-Herrero et al., 1985). Mutations in Hox genes in different species produce AP homeotic transformations, which in mammals is best exemplified in the hindbrain and in the axial skeleton (Krumlauf, 1994).

Hox gene transcription is activated sequentially in axial precursors during gastrulation (Deschamps et al., 1999). As cells derived from these precursors colonize the different AP segments, they carry the successive Hox expression combinations to the progressively forming body axis, resulting in an AP nested patterns (Alexander et al., 2009). Thus, temporal information is translated into spatial domains during axial elongation (Deschamps and Duboule, 2017).

HOX proteins bind DNA through a highly conserved 60 amino acid region called the homeodomain (McGinnis et al., 1984). HOX proteins alone show limited DNA-binding ability, but they gain specificity and affinity for target sequences through interactions with co-factors of the PBC and MEINOX families, both belonging to the three amino acid loop extension (TALE) class of homeodomains (Mann and Affolter, 1998). PBC and MEINOX proteins form heterodimers and heterotrimers with HOX proteins, conferring them with increased target sequence selectivity and affinity (Merabet and Mann, 2016). Mutation of single members of the PBC and MEINOX families in Drosophila show AP phenotypes compatible with a generalized Hox gene loss of function, without affecting Hox AP expression (Chan et al., 1994; Rieckhof et al., 1997). In mammals, redundancy of the PBC (four members) and MEINOX (five members) families has hampered the study of their roles in axial skeletal patterning. Although knowledge has been obtained from Pbx mutants in zebrafish and mouse, indicating essential roles in axial skeleton patterning (Capellini et al., 2008; Popperl et al., 2000; Selleri et al., 2001), the role of Meis genes in this context remains unexplored.

MEIS proteins directly bind HOX proteins encoded by paralogs 9-13 (Shen et al., 1997) and form DNA-bound heterotrimeric complexes with PBX and HOX proteins encoded by paralogs 1-10 (Chang et al., 1996). The repertoire of MEIS-, PREP- and PBX-binding sites revealed by ChIP-seq analysis in E11.5 mouse embryos identified HOX and HOX-PBC binding sites as the preferred sites for MEIS binding, above the MEIS-only binding sites, suggesting that MEIS factors are strongly dedicated to interactions with HOX and PBX proteins (Penkov et al., 2013). In addition, a large number of MEIS-binding sites was found within the Hox gene complexes, which suggested that, in addition to their HOX co-factor role, they may regulate Hox gene transcription. Studies in zebrafish (Choe et al., 2014) and mouse (Amin et al., 2015) embryos indeed show that some of these binding sites represent Hox auto-regulatory elements that are active in the neural tube, and Hox transcription regulation by MEIS factors has been demonstrated in neural tube patterning (Dibner et al., 2001; Vlachakis et al., 2001; Waskiewicz et al., 2001). MEIS elimination or overexpression also affects Hox gene expression during limb skeletal patterning (Delgado et al., 2020; Mercader et al., 1999,, 2009; Rosello-Diez et al., 2014); however, this aspect has not been studied in axial skeleton patterning.

An important pathway connecting MEIS, HOX and axial patterning is that of vitamin A. The active form of vitamin A, retinoic acid (RA), regulates gene expression during embryonic development by binding to nuclear receptors RARα, RARβ or RARγ (Rhinn and Dolle, 2012). Meis genes have been identified in screens for RA targets (Berenguer et al., 2020; Oulad-Abdelghani et al., 1997) and respond to RA fluctuations in vivo (Mercader et al., 2000; Yashiro et al., 2004). RA excess produces axial skeleton alterations and modifies the Hox AP expression domains (Kessel and Gruss, 1991) and mutations in RA-receptor genes result in homeotic transformations; however, the mechanism by which this takes place is not clear. Although RAR binding sites have been described in Hox complexes (Marshall et al., 1996), and RA administration in vitro regulates Hox gene transcription (Deschamps et al., 1987), RA administration in vivo can lead to axial skeleton homeotic transformations without changes in Hox expression (Kessel, 1992) and changes in Hox expression in Rar-deficient mice have not been reported.

Here, we study the role of MEIS factors in axial skeleton formation and its relationship to the RA pathway by characterizing mouse genetic models of Meis1, Meis2 and Raldh2 (the synonym for Aldh1a2). We dissect the regulatory and functional relationships between Meis genes, Hox genes and Raldh2, and formulate a new model that explains the ability of RA to produce homeotic transformations without modifying Hox expression.

We studied the mRNA expression pattern of Meis1 and Meis2: the two Meis genes extensively expressed in paraxial and lateral mesoderm (Fig. 1). We detected the earliest expression of Meis2 in early-streak stage embryos in a posterior region of the embryo close to the boundary with the extra-embryonic region (Fig. 1H). This expression extends distally and anteriorly as development progresses (Fig. 1I); at the early-headfold stage, an anterior stripe of Meis2 transcripts was found bilaterally close to the extra-embryonic region, and continuous with its posterior expression (Fig. 1J). At late-headfold stage, Meis2 started to disappear from the posterior region (Fig. 1K) and at E8, the posterior embryonic bud was devoid of Meis2 transcripts (Fig. 1L). Meis1 expression started slightly later than Meis2, being first detected at the late-streak stage, bilaterally in the mesoderm close to the extra-embryonic region (Fig. 1B), and at the early-headfold stage, forming a stripe of expression similar to the Meis2 anterior stripe at this stage (Fig. 1C). Both Meis1 and Meis2 expression domains extend posteriorly into the lateral plate mesoderm at the late-headfold stage (Fig. 1D,K), but high levels of Meis1 transcripts were never observed in the posterior embryonic bud. Finally, at E8 Meis1 and Meis2 expression patterns converge to a similar domain, being strongly expressed in paraxial and lateral plate mesoderm up to the pharyngeal region (Fig. 1E,L). At this stage, expression of both genes is excluded from the posterior embryonic bud, whereas it appears in the presomitic mesoderm and adjacent regions precursor to the lateral plate mesoderm. This expression pattern is maintained at later stages, indicating that as new precursors from the posterior bud incorporate to the presomitic area, they activate Meis1 and Meis2, and this activity persists as they differentiate into paraxial and lateral plate mesoderm. To determine the early activation pattern of Meis1 and Meis2 in the embryonic germ layers, we studied Meis mRNA and protein distribution in sections (Fig. 1). Detection of MEIS proteins in sections with an antibody that recognizes the majority of embryonic isoforms, but does not discriminate between MEIS1 and MEIS2, shows early expression in all three germ layers at early allantoic bud stage (Fig. 1G,N). Sections of the RNA in situ hybridization of both genes, showed that Meis1 expression was not detected in the epiblast/ectoderm (Fig. 1F), while Meis2 expression affected the three germ layers (Fig. 1M). This result suggests Meis2 is activated in epiblast cells and its expression persists as they gastrulate to contribute to mesoderm. The early Meis2 expression pattern thus resembles the activation pattern of Hox genes.

We used conditional deletion of Meis1 and Meis2 and studied the mutant skeletal pattern. We first studied the consequences of Meis2 deletion using different Cre alleles that allow dissecting the putative specific functions of Meis2 early expression. Deletion of a Meis2^flox allele with Sox2^Cre leads to Meis2 elimination in the epiblast. Lethality of Sox2^Cre;Meis2^flox/flox embryos around E14.5-E15.5 due to cardiac defects did not allow us to study the pattern at later stages; however, the general vertebral formula could be determined at E14.5. We observed defects at the occipital-cervical transition, where the first cervical vertebra (C1 or atlas) was fused to the exoccipital bone (n=14/14) and in its ventral part showed a position and shape that resemble the exoccipital bone, while its dorsal part was not formed (Fig. 2A,B; Fig. S1; Table S1A). These changes correlated with a change in the shape of the second vertebra (C2 or axis), which acquired a C1-like morphology (n=13/14) (Fig. 2A,B; Fig. S1). With low penetrance, the C3 vertebra presented a morphology that resembles C2 (n=2/14). These observations are compatible with an anterior homeotic transformation of the cervical vertebrae. In addition, disconnected isolated elements often appeared (arrowhead in Fig. 2B; Fig. S1), suggesting segmentation problems are also present in this region. Outside the axial skeleton, we observed a vestigial otic capsule in the mutants.

In the thoracic region, the most prominent defect was rib, rib-sternum attachment and sternum mispatterning (Fig. 2D,E). We observed failures in sternum fusion, rib bifurcations, fusions and alteration of the sternal/floating rib formula. The gain of a rib in the first lumbar vertebra (L1) in some specimens (n=4/14) and the tendency to reduction of the first rib (R1) suggests the anterior transformations observed in the cervical region may also affect the thoracic region (Fig. 2B). More caudal regions did not show any defects.

To investigate whether Meis2 activity in the epiblast is involved in the observed defects, we combined the Meis2^flox allele with Mesp1^Cre to eliminate Meis2 from the nascent mesoderm. Although Mesp1 activates in the early embryo in a similar pattern to Meis2, because of the time lag between Cre expression and effective recombination, the recombination pattern of Mesp1^Cre affects only the mesoderm and it does so down to the forelimb level (Saga et al., 1999). As occurred in Sox2^Cre;Meis2^flox/flox mice, lethality due to cardiac defects allowed us to study the phenotype only until E14.5. In the Mesp1^Cre model, we observed lower penetrance, but the same type and distribution of defects found in the Sox2^Cre model, except for the reduction of R1 and the otic capsule defects (Fig. 2C,F; Fig. S1; Table S1A).

To further dissect the specific tissues in which Meis2 activity is required during early embryogenesis, we studied a third model in which we deleted Meis2^flox using Dll1^Cre, a line that recombines the mesoderm in the presomitic region (Wehn et al., 2009), i.e. at a later step of mesodermal allocation than does Mesp1^Cre. Dll1^Cre;Meis2^flox/flox mice survive to adulthood, allowing a full assessment of the skeletal pattern at the end of gestation. In this model, we observed similar defects to those previously observed in the Sox2^Cre and Mesp1^Cre models in the occipital, cervical and thoracic regions (Figs 2G-N and 3M; Fig. S1; Table S1A). In addition, we observed a defect in supraoccipital ossification (Figs 2K,M and 3M) and fusions between the basioccipital and the anterior arch of the atlas (aaa) (Fig. 2G,I; Fig. S1 and Fig. 3M), which could not be determined at earlier stages because these bones form late in gestation. Again, we did not detect the formation of a rib in L1, suggesting this phenotype requires an early deletion of Meis2.

The irrelevance of early Meis2 expression for most aspects of axial patterning is not due to compensatory activation of Meis1, as we detected no ectopic Meis1 mRNA expression in early Sox2^Cre;Meis2^flox/flox embryos (Fig. S2). These results indicate that the expression of Meis2 in the epiblast and early nascent mesoderm is to a large extent dispensable for its functions in axial skeletal patterning, although it might be needed for a proper specification of the thoracic-to-lumbar transition.

Next, to determine whether Meis1 and Meis2 cooperate in axial patterning, we combined Meis1 and Meis2 mutant alleles. Combining Meis1 and Meis2 deletion is not possible using the Sox2^Cre or the Mesp1^Cre deleters, owing to lethality of double heterozygous mice. We therefore used the Dll1^Cre line for these experiments. The defects observed in the allelic series generated affected the same skeletal elements that were altered in the Meis2 mutant models (Fig. 3) and the type of defects were similar, with anterior transformations of C1-C3 (Fig. 3G-I,M; Fig. S3; Table S1B,C) and defects in the occipital bones that either did not form or appeared fused to the C1 or C2 element (Fig. 3A-I,M; Fig. S3; Table S1B,C). In the thoracic region, we also detected rib fusions and defects in rib-sternum attachment (Fig. 3J-M; Table S1B,C). Although we found two cases of extra ribs on L1, one case was also found in controls, suggesting this observation was nonspecific. In general, skeletal defects become more severe as the number of Meis alleles deleted increases, with the absence of Meis2 being more detrimental than absence of Meis1 (Fig. 3M; Fig. S3; Table S1B,C). However, for some aspects of the phenotype, E18.5 Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox specimens appeared less affected in comparison with Dll1^Cre;Meis1^+/flox;Meis2^flox/flox specimens, which was paradoxical. However, we observed, that the viability of Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox mice at E18.5 was 37%, which suggested that specimens of this genotype at E18.5 represent escapers; thus, missing specimens could be more affected than appreciated. We then studied the phenotype of Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox fetuses at E14.5, when viability of double mutants was 67%, and observed a fraction of embryos with defects compatible with those observed at E18.5 and, in addition, we found very strongly affected fetuses in which all the cervical vertebrae fused, there was no apparent development of occipital condensations, and there were widespread rib fusions and truncations (Fig. 3N,O). In addition to the axial skeleton, defects in the limb skeleton were obvious in Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox fetuses (Fig. 3N,O) and have been independently reported (Delgado et al., 2020).

Although defects in vertebral segmentation or overt rib deformations have not been described in Hox mutants, the AP specification defects observed in the occipital and cervical regions are very similar to those observed in the mutants of Hox paralog groups 3-5, and some of the rib cage defects are similar to those found in paralog groups 5-9 (Horan et al., 1995; Jeannotte et al., 1993; Manley and Capecchi, 1997; McIntyre et al., 2007). These coincidences and the previous reports of MEIS proteins binding to Hox gene clusters (Amin et al., 2015; Penkov et al., 2013) prompted us to study the Hox mRNA expression pattern in Meis mutants. We did not detect alterations of Hox expression initiation or definitive anterior expression borders in the paraxial mesoderm in either Sox2^Cre;Meis2^flox/flox or in Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox embryos (Fig. 4A,B). These results indicate that eliminating Meis2 function with Sox2^Cre or Meis1 and Meis2 Dll1^Cre does not modify Hox gene expression patterns and therefore, the phenotypes observed in these models do not relate to a role for Meis genes in regulating Hox gene transcription in the paraxial mesoderm. To study the generality of these observations, we combined simultaneous maternal and paternal deletion of Meis1^flox/flox and Meis2^flox/flox alleles, using the maternal deleter Zp3^Cre and the paternal deleter Stra8^Cre (Fig. S4). Using this approach, we were able to completely eliminate Meis1 and Meis2 zygotic expression. Such embryos die around E9 with profound alterations of cardiac development; however, this allowed us to study early Hox expression patterns. Although previous reports in embryos at E9.5 or later stages have described paralog group Hox3 gene expression starting at somite 5 (Alexander et al., 2009), we found that in control embryos of up to 10 somites, expression of the Hox3 paralog group extended from somite 2-3 into more posterior somites (Fig. 4C, Fig. S5). In embryos of 12 somites, the most anterior expressing somite is somite 3-4, while in embryos of 15 or more somites, expression starts at somite 5. These observations show transient Hox3 gene expression in occipital somites and a later progressive posteriorization towards their definitive expression domain. This expression pattern is consistent with the observation that the occipital region, which mostly originates from somites 1-4, is strongly affected in Hox3 mutants. The defects present in Hox3 mutants in fact strongly affect the supraoccipital bone, which is exclusively contributed to by somites 1 and 2 (Huang et al., 2001; Muller and O'Rahilly, 1994). Meis mutant embryos showed a normal Hox3 group gene expression in the paraxial mesoderm of embryos of 4-10 somites (Fig. 4C, Fig. S5). In more advanced stages, although counting somites was very difficult in mutant embryos of 15-20 somites, due to the developmental abnormalities, we concluded that the expression patterns in the paraxial mesoderm were either normal or anteriorized by 1-2 somites (Fig. 4C, Fig. S5). In contrast, the anterior border of expression in the neural tube appeared clearly posteriorized (Fig. 4C, Fig. S5). The study of the expression of Hoxd4 showed similar results, with a transient early expression starting at somite 4 and later becoming restricted to its definitive anterior border at somite 6. In mutants, Hoxd4 expression was similar at early stages and appeared anteriorized to somite 4-5 at later stages. A posteriorization of Hox mRNA expression in the neural tube was again evident (Fig. 4C, Fig. S5). Most likely, the failure of the mutants to relocate Hox expression to more posterior somites at late stages does not indicate a direct role for MEIS proteins in regulating Hox gene expression, but a general blockade in development of Meis double knockout (DKO) embryos beyond the somite 7 stage. In fact, Meis DKO embryos do not undergo turning, body wall folding or neural tube closure, morphologically resembling E8.5 embryos at E9. The fact that no alterations were observed upon deletion with Dll1^Cre support this conclusion. We therefore conclude that transcriptional regulation of Hox genes is not involved in MEIS protein-mediated regulation of axial skeleton patterning.

To identify the molecular mechanism underlying the skeletal phenotypes observed, we performed a transcriptomic analysis of E9 Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox embryos. To discriminate between alterations during the somite differentiation phase and early patterning defects, we separately analyzed the anterior region containing the first 10-12 somites and the posterior region, including the rest of the somites and the posterior embryonic bud. We identified 9 upregulated genes and 25 downregulated genes in the analysis of the anterior region, whereas in the posterior region there were 58 upregulated and 58 downregulated genes (Fig. S6A,B). Regulatory pathway analysis showed that ‘skeletal and muscular system development’ appears as the top tissue-specific altered class (Fig. S6C). Differences in other processes, such as cell death, cell-to-cell interactions, cell assembly and organization, were also found in this analysis (Fig. S6C). No alterations were found in Hox gene expression, which confirmed the results observed in the Hox mRNA in situ analysis.

We then focused on the in situ analysis of genes involved in somite development found altered in the RNAseq analysis and in additional genes relevant to somite patterning. When comparing the expression pattern of this set of genes between control and Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox embryos, we found that a set of genes expressed and/or involved in hypaxial myotomal development were downregulated in the hypaxial region of the mutants (Fig. 5), including Eya1 (Grifone et al., 2007) (Fig. 5A,B), Sim1 (Ikeya and Takada, 1998) (Fig. 5C,D), Shisa2 (Nagano et al., 2006) (Fig. 5E,F) and Pax3 (Tremblay et al., 1998) (Fig. 5G,H). Regarding sclerotome markers, we found no alteration of Pax1 expression (Fig. 5I,J); however, an abnormal expression pattern of Pax9 was observed in sclerotomes of the cervical region, which appeared incorrectly segmented (Fig. 5K,L).

Crosstalk between the myotome and sclerotome is essential for proper rib patterning and mice deficient for the myogenic factors MYF5 (Braun et al., 1992), MRF4 (Zhang et al., 1995) or myogenin (Vivian et al., 2000) show rib defects that resemble those described here in Meis mutants. Therefore, we next studied the main myogenic factors. Expression of Myf5 appears first in the epaxial somite at E8, followed by Mrf4 and Myog at E9, later extending hypaxially caudal to the forelimb at E10.5 (Fig. 5M-R). In mutant mice, the early epaxial expression of Myf5 shows incomplete segmentation, whereas at E10.5, expression in myotomes anterior to the forelimb extends ventrally and appears as a continuous band between adjacent somites in a pattern that is not detected in control embryos (Fig. 5M,N). Both Mrf4 and Myog show mis-segmented and bifurcating patterns in mutant embryos (Fig. 5O-R). In addition, the ventral hypaxial extension of the expression domain was reduced, as observed before for other hypaxial markers. In contrast, defects in the early expression of Myog at E9 are not as evident as for Myf5. Myod1 shows as well a disorganized and spread expression in cervical myotomes of mutants, whereas hypaxial extension of the expression is also defective in more caudal myotomes (Fig. 5S-T).

We finally studied Fgf4 and Fgf6, which are involved in myogenesis through their expression in the medial myotome (Grass et al., 1996). We found that expression of Fgf4 and Fgf6 appeared highly reduced in mutant embryos (Fig. 5U-X).

In summary, re-segmentation of the paraxial mesoderm appears impaired in Meis mutants, with defects in the separation of adjacent sclerotomal/myotomal domains and bifurcated myogenic domains. These defects affect mainly the cervical region, although defects were also seen sometimes in the interlimb region. During further myotome development, a defect in myogenic Fgf gene expression was found and the hypaxial developmental program seemed especially affected, with a failure in hypaxial myotomal migration in correlation with an inability to properly activate Pax3 expression. To study whether any of the skeletal or myogenic defects found correlate with alterations of the segmentation clock, we analyzed the expression of Hes5 and Lnfg – factors involved in cyclic expression in the presomitic mesoderm; however, we did not find any alterations in the expression of the mRNAs for these factors (Fig. S7).

In the transcriptomic analysis of Meis mutants, Raldh2 – the gene encoding the main enzyme responsible for embryonic RA synthesis – and Cyp26b1 – the gene encoding the main enzyme responsible for RA degradation in the embryo – appeared to be downregulated in the anterior trunk region by RNA-seq (Fig. S6). In situ hybridizations for both genes were consistent with the transcriptomic analysis. Raldh2 expression appeared reduced at E9.5 in the differentiating derivatives of anterior somites but not in the presomitic area or in newly produced somites (Fig. 6A-B″). A similar pattern is present at E10.5, where tail regions with newly produced somites do not show alterations, but more anterior regions do show a reduction in Raldh2 transcripts (Fig. 6C-D″).

In the somitic region of E9 embryos, Cyp26b1 is expressed exclusively in the endothelium of the dorsal aortae and inter-somitic vessels, whereas it is strongly expressed in areas of the hindbrain. The hindbrain signal was preserved in mutants; however, the endothelial signal in the trunk region was lost (Fig. 6E-F′). Cyp26b1 is a direct target of the RA pathway that is activated in response to RA. The concomitant downregulation of Raldh2 and Cyp26b1 thus suggest that Meis mutant embryos have defective RA. We then studied the expression of the gene encoding the RA β receptor (Rarb), which has been described as a RA-responsive gene. Contrary to expectations, no change in the pattern of Rarb mRNA was detected between controls and Meis mutants (Fig. 6G-H′), which is consistent with the RNAseq analysis, which identified no differences in Rara, Rarb or Rarg. The observations on the expression levels of Raldh2, Cyp26b1 and Rarb mRNAs were also confirmed by quantitative PCR (Fig. S8).

Unexpectedly, several of the embryos studied showed Raldh2 reduction in a mosaic fashion. To understand why the reductions in Raldh2 appeared mosaic, we combined Meis1^flox and Meis2^flox alleles with Dll1^Cre and a Rosa26R^tdTomato reporter. In these embryos, Dll1^Cre recombines the Rosa26R^tdTomato reporter, allowing the Cre recombination pattern to be determined. At E10.5, E9.5 and E8.5, we observed a mosaic pattern of Tomato⁺ cell distribution in both control and Meis mutant embryos, with variability in the proportion of Tomato⁺ cells found in the somites of different embryos (Fig. S9). This mosaicism had not been described for this line before (Wehn et al., 2009), and therefore it might depend on the genetic background. To determine whether the observed mosaicism results from inefficient recombination in all cells or from mosaic activation of Dll1^Cre expression, we studied the correlation between Meis immunodetection and Tomato expression in Dll1^Cre;Rosa26R^tdTomato;Meis1^flox/flox;Meis2^flox/flox embryos. We found that Tomato⁺ cells were devoid of Meis, whereas their neighboring Tomato⁻ cells showed Meis expression (Fig. 6I-J‴). Image profiling shows anti-correlation between Tomato and Meis detection in mutants (Fig. 6J‴), whereas this was not found in control embryos (Fig. 6I‴). These observations indicate that the pattern observed results from mosaic inactivation of Dll1^Cre and therefore the Tomato⁺ cell distribution reports the distribution of Meis-deficient cells. In mutants, we found a tendency of knockout and wild-type cells to segregate from each other, resulting in large aggregates of Tomato⁺ cells that were not found in controls (Fig. 6I″,J″). We did not find any reproducible difference between mutant and control embryos in the distribution of Tomato⁺ cell patches by tissues. In addition, the anterior-most border of Tomato⁺ cell distribution was established at the occipital level and this did not differ between control and mutant embryos.

We therefore used the mosaic inactivation of Meis alleles to study the regulation of Raldh2 by MEIS. We performed Raldh2 immunostaining and correlated this signal with that of Tomato. We found that Tomato⁺ cells lacking MEIS function did not present detectable RALDH2 expression, while their Tomato⁻, Meis-expressing neighboring cells showed normal Raldh2 expression (Fig. 6K-L‴). The result was similar to that observed for MEIS immunostaining, with the signals for RALDH2 and Tomato being mutually exclusive in mutant embryos but not in controls (Fig. 6K″,L″). These results indicate a strict and cell-autonomous requirement for Meis function in Raldh2 expression in the differentiating trunk mesoderm.

We then analyzed Raldh2 expression in embryos with double maternal/zygotic inactivation of Meis1 and Meis2 (Fig. 6M-R). Raldh2 mRNA distribution in the early embryo resembles Meis expression pattern; however, it starts slightly later and affects only the mesoderm (Fig. S10). In mutant embryos, we observed no alteration of the expression pattern in the axial and paraxial mesoderm; however, the lateral plate domain close to the extra-embryonic region was abolished (arrowheads in Fig. 6M,N). Up to E8.75, no alteration of Raldh2 expression in the paraxial mesoderm was observed (Fig. 6O,P); however, at E9, all trunk Raldh2 expression was strongly decreased in mutants (Fig. 6Q,R). The reduction in Raldh2 and Cyp26b1 was confirmed by qPCR in four out of five embryos with complete absence of Meis1 and Meis2 (Fig. S8); however, no reduction was observable for Rarb by either RNA in situ hybridization or qRT-PCR (Fig. 6G-H′,S,T; Fig. S8). Although RALDH2 is the main enzyme for RA synthesis in the embryo, RALDH3 is also present in embryos and contributes to RA generation. To determine whether Raldh3 activation could contribute to RA synthesis, we studied transcript abundance by qRT-PCR and found instead reduced expression in three out of six embryos in the Dll1Cre model and in all five embryos with complete absence of Meis1 and Meis2 (Fig. S8). These results show that Meis controls RA synthesis by maintaining transcription of genes encoding essential RA-synthesizing enzymes. Contrary to expectations, this does not result in reduced transcription of the RA direct target Rarb.

Given that retinoic acid has been shown to regulate Meis expression in different settings (Berenguer et al., 2020; Mercader et al., 2000; Oulad-Abdelghani et al., 1997; Yashiro et al., 2004), we studied whether the elimination of Raldh2-mediated RA synthesis affects axial Meis expression. We studied Meis1 and Meis2 mRNA expression in Sox2^Cre;Raldh2^flox/flox embryos, and found that both genes presented lower levels of transcripts along the trunk region of E8 embryos (Fig. 6U-X).

These results indicate that Meis is required for maintenance of Raldh2 expression in the differentiating paraxial mesoderm but not for its initial expression before somite differentiation. These conclusions correlate with the observed downregulation of Raldh2/Cyp26b1 in the transcriptome of the anterior trunk but not the posterior trunk of E9.5 embryos. In contrast, the early lateral plate mesoderm – likely fated to the cardiogenic area – requires Meis activity for Raldh2 expression from the earliest stages. Reciprocally, Raldh2 expression is required to maintain proper Meis expression levels, but not for initiating Meis expression, given that Meis expression starts before Raldh2 expression. These results indicate that there is a positive regulatory loop between Meis and Raldh2 that is relevant to mutually maintain but not initiate their expression.

Although retinoic acid has long been postulated as a regulator of axial skeleton, there is no direct study of the consequences of eliminating RA on anterior-posterior axial identities. Here, we conditionally deleted Raldh2 using Dll1^Cre to investigate whether this affects the axial skeleton and the extent to which RA might be related to Meis roles in axial patterning. In the occipital region, the basioccipital presented similar alterations to those observed in Meis mutants (n=14/43) (Fig. 7A,F; Table S1D,E), including its fusion with the aaa (arrowhead in Fig. 7F). Strikingly, similar modifications of the basioccipital were also found in some control embryos, although in a lower proportion (n=5/47) (Fig. 7F,P), suggesting a genetic background prone to these particular defects. In mutants, C1 appeared fused to, and/or adopting a shape and position similar to, the exoccipital (n=8/49) (Fig. 7H,P; Fig. S11). In the cases in which C1 showed transformation to exoccipital, C2 adopted a C1 morphology (n=8/41), whereas some cases in which C1 retained its morphology, C2 adopted a C1 morphology and partially fused to C1 (n=9/41) (Fig. 7H,P; Fig. S11). C3 to C2 transformations/fusions were also observed (n=11/41) (Fig. 7P). At the cervical thoracic transition, tuberculi anterior were found in C7 instead of C6 (n=5/21) (Fig. 7C,H, arrow in 7H; Fig. S11), suggesting that anterior transformations also take place at this axial level. Altogether, the alterations found in Raldh2 mutants in the occipital/cervical regions were similar to those observed in Meis mutants but displayed lower penetrance (Fig. 7P).

In the thoracic region, shortening or fusion of the first rib with the second rib and generalized rib fusions and bifurcations were observed, similar to the defects found in Meis mutants (Fig. 7D,I). In the most affected mutant embryos, we observed defects in the inter-sternal cartilage and the sternebrae, although we did not observe a split sternum (Fig. 7E,J). Some incidences of an extra sternal rib and an extra rib on L1 were also observed (Fig. 7E,J,P), suggesting AP transformations were extensive down to the thoracic/lumbar transition.

The compared analysis of Meis and Raldh2 mutants supports the idea that MEIS and the retinoic acid pathway act in a positive-feedback loop that is relevant in patterning the axial skeleton. To obtain evidence for the functional relevance of this regulatory loop and determine its relevant output in axial patterning, we used a Rosa26R^Meis2-EYFP allele (Rosello-Diez et al., 2014) that provides Meis2 overexpression coupled to IRES-EYFP expression upon Cre recombination. We then simultaneously eliminated Raldh2 and activated Meis2 with Dll1^Cre. Co-staining for detection of GFP and RALDH2 confirmed deletion of Raldh2 concomitant with Meis2 activation (Fig. S12). Interestingly, in this mouse model, all defects produced by Raldh2 mutation in the axial skeleton were rescued (Fig. 7K-P; Table S1D,E), indicating that Meis suppresses the effect of RA deficiency on axial skeleton patterning.

Although Meis1 and Meis2 expression starts at gastrulation, their early patterns are different in time and expression domains. Meis2 is activated earlier than Meis1 in a pattern that coincides spatially and temporally with that of Hox gene activation in the posterior epiblast. We have not observed major alterations in Hox gene expression patterns or transcript abundance in Meis mutants. These results indicate that, despite the profuse binding of Meis proteins to the Hox complexes (Penkov et al., 2013), Meis is not essential for Hox gene transcriptional regulation during axial skeleton patterning. This does not apply to other embryonic regions, given that we have observed clear alterations of Hox mRNA expression domains in limb buds with altered Meis function (Delgado et al., 2020; Mercader et al., 1999, 2009; Rosello-Diez et al., 2014) and similar results have been reported during neural patterning in other vertebrate species (Dibner et al., 2001; Vlachakis et al., 2001; Waskiewicz et al., 2001) and are reported here in the mouse (Fig. 4). The observed binding of MEIS to the Hox complexes might therefore be involved in regulating Hox transcription in several tissues but not in the paraxial mesoderm, at least at the stages studied here.

Despite the absence of changes in Hox transcription in the paraxial mesoderm, Meis mutants produce anterior homeotic transformations and defects similar to those previously described for Hox mutants affecting the occipital and cervical regions (Horan et al., 1995; Jeannotte et al., 1993; Manley and Capecchi, 1997; McIntyre et al., 2007). This is consistent with studies in flies in which the elimination of the Meis ortholog homothorax produces homeotic phenotypes by modifying the DNA affinity and target selectivity of HOX proteins without altering Hox gene transcription (Merabet and Mann, 2016; Rieckhof et al., 1997). We thus propose that, during axial patterning in vertebrates, Meis regulates Hox function by directly binding HOX proteins and regulating their affinity and selectivity for target DNA sequences, through their described HOX-co-factor activity.

We deleted Meis2 using Cre lines that recombine at different stages of epiblast cell recruitment to the paraxial mesoderm; however, we did not find any substantial influence of the timing of Meis2 removal on the phenotypes obtained. The early expression of Meis2 in the posterior epiblast thus seems not to play any role in axial skeletal patterning, whereas both Meis1 and Meis2 cooperate in the presomitic mesoderm, or at later stages of somite development, in axial patterning.

Apart from its function in segmental identity, the transcriptional analysis of the mutants indicates an important function in hypaxial myotome development, with profound alterations of the myogenic pathways and phenotypes that have not been reported in Hox mutants. Interestingly, Myf5, Mrf4 and myogenin-deficient mice show rib defects similar to those described here (Braun et al., 1992; Hasty et al., 1993; Patapoutian et al., 1995), and, therefore, the failure in proper activation of the hypaxial myogenic program is sufficient to explain rib mispatterning in Meis1/Meis2 DKOs. Moreover, hypaxial myotomal Fgf4 and Fgf6 expression, which is required for rib patterning downstream of the myogenic factors (Vinagre et al., 2010), is strongly impaired in Meis mutants, indicating a function of MEIS in the crosstalk between myotome and sclerotome. In addition, the activation of the myogenic program involved in rib pattering is under direct control of a specific set of HOX proteins involved in the specification of thoracic segments (Vinagre et al., 2010). The rib mispatterning phenotypes may therefore also partly involve the impairment of Hox function in the absence of MEIS.

RNA-seq analysis, quantitative PCR and/or in situ hybridization revealed a reduction in Raldh2, Raldh3 and Cyp26b1 in Meis mutants. RALDH enzymes are essential for RA synthesis and, therefore, Meis mutant embryos presumably contain reduced levels of RA. As the activation of Cyp26b1 is RA dependent, its downregulation in Meis mutants could be a secondary event, due to the reduction in RA synthesis by RALDH2. A contradictory result is the observation that Rarb, a well described target of RA, does not show reduced expression in Meis-deficient embryos. We do not have a clear explanation for this result; it could indicate a lower RA threshold for Rarb activation than for Cyp26b1 activation or a more complex situation in which its transcriptional dependence on RA is modified by the absence of Meis.

Cyp26b1 mutants show posterior homeotic transformations in the occipital/cervical region associated with increased RA levels (Sakai et al., 2001). The posterior transformations in this model are the opposite of those observed in Meis mutants, which concurs with the idea that Meis is a positive regulator of RA synthesis. In addition, in vivo treatments with RA during mouse gestation caused either anterior or posterior homeotic transformations, depending on the stage of the treatment (Kessel and Gruss, 1991). In the cervical region, anterior transformations were observed following treatments at E7, while posterior transformations were found following RA treatment from E8 (Kessel, 1992). Here, we generated a Raldh2 conditional knockout using the Dll1^Cre driver and found homeotic transformations affecting the occipital/cervical region, and additional patterning defects in the thoracic region that significantly overlap with those observed in Meis mutants. In agreement with this, mutations in retinoic acid receptor genes produce homeotic transformations (Lohnes et al., 1993, 1994), similar to those observed in Meis mutants. In particular, Rarg and Rarb loss-of-function mutants show anterior transformations without showing any changes in Hox expression patterns (Folberg et al., 1999a,b).

We described a positive regulatory loop between MEIS and RALDH2 that we confirmed using functional genetic analysis. The similarities in skeletal transformations between Raldh2 and Meis mutants and the cross-regulation between Raldh2 and Meis suggest that the positive regulatory loop between Raldh2 and Meis is involved in axial patterning. Although Meis genes are RA targets in various contexts (Mercader et al., 2000; Oulad-Abdelghani et al., 1997; Yashiro et al., 2004), Raldh2 is a direct Meis target in the hindbrain (Vitobello et al., 2011), and ChIPseq analysis in E10.5 limb buds identifies MEIS-binding sites in the Raldh2 locus (Delgado et al., 2020). In fact, MEIS could promote RA accumulation at various levels, as it also represses Cyp26b1 during limb development in a cell-autonomous manner (Rosello-Diez et al., 2014). The requirement of Meis activity for Raldh2 transcription in the paraxial mesoderm is restricted to the differentiation stages and does not take place in the nascent or segmenting mesoderm. In concurrence with our findings, Pbx1/Pbx2-null embryos show normal Raldh2 expression at early embryonic stages, but strong downregulation in the paraxial mesoderm at E9.0 and beyond (Vitobello et al., 2011). Finally, we studied the functional output of the MEIS-RA regulatory loop by genetic rescue. The complete rescue of Raldh2 mutants by Meis overexpression suggest that MEIS is the main functional output of the positive regulatory loop between Meis and RA in the paraxial mesoderm.

We propose a model for the RA-Meis-Hox network in the paraxial mesoderm in which Meis is involved in a positive-feedback loop with RA through Raldh2 regulation (Fig. 8). Meis is the main output of this regulatory loop and is required for the specification of axial skeletal identities, likely through regulating HOX protein activity. The proposed model would explain the ability of RA and RARs to phenocopy Hox mutants without affecting their transcriptional expression. In addition, as deduced from the segmentation and myogenic phenotypes, Meis factors play roles unrelated to Hox function.

Experiments were performed using mice (Mus musculus). Mice were handled in accordance with CNIC Ethics Committee, Spanish laws and the EU Directive 2010/63/EU for the use of animals in research. All mouse experiments were approved by the CNIC and Universidad Autónoma de Madrid Committees for ‘Ética y Bienestar Animal’ and the area of ‘Protección Animal’ of the Community of Madrid with reference PROEX 220/15.

Meis conditional knockouts were generated mating Meis1^flox (Unnisa et al., 2012) and Meis2^flox (Delgado et al., 2020) with different Cre lines: Sox2^Cre (Hayashi et al., 2002), Mesp1^Cre (Saga et al., 1999), Dll1^Cre (Wehn et al., 2009), Stra8^Cre (Sadate-Ngatchou et al., 2008) and Zp3^Cre (de Vries et al., 2000). Raldh2 conditional knockouts were obtained by mating Raldh2^flox mice (Vermot et al., 2006) with Dll1^Cre and Sox2^Cre. For Cre⁺ cell-lineage tracing, we used Gt(ROSA)^{26Sortm14(CAG-tdTomato)Hze} (Madisen et al., 2010). For conditional Meis overexpression we used the Rosa26R^Meis2-EYFP line (Rosello-Diez et al., 2014).

To obtain embryos at different gestational stages, mice were mated in the afternoon and females were checked every morning for the presence of a vaginal plug; noon on the day the plug was observed was considered as gestational day 0.5 (E0.5). Embryos at somitogenic stages were staged according to age and somite number. Embryos that had not started somitogenesis were staged according to Downs and Davies (1993).

Embryos were fixed in 4% PFA overnight at 4°C. Embryos were dehydrated and rehydrated washing them with increasing and decreasing, respectively, concentrations of methanol in PBT (25%, 50%, 75% and 100%). Bleaching was carried out by incubation in 6% H₂O₂ in PBT for 1 h. Proteinase K (Sigma) digestion was performed at 10 µg/ml with different incubation times depending on the stage. After permeabilization, embryos were washed with PBT for 5 min and fixed with 0.05% glutaraldehyde in 4% PFA for 20 min. Embryos were incubated in hybridization buffer [50% formamide, 4× SSC (pH 4.5), 1% SDS, 50 µg/ml heparin (Sigma), 10 µg/ml tRNA from baker's yeast (Sigma) and 1% w/v Blocking reagent (Sigma)] for 2 h and hybridized with the probe overnight at 65°C. Post-hybridization washes were performed with 0.1% w/v CHAPS (Sigma) and 2× SSC (pH 5.5), followed by a second round of post-hybridization washes with 0.1% w/v CHAPS and 0.2× SSC for 3 h at 65°C. Embryos were incubated overnight at 4°C with 1:2000 anti-digoxigenin AP antibody (Roche) in 20% goat serum, 1% blocking reagent in TBST [5 mM Tris-HCl (pH 7.5), 15 mM NaCl and 0.1% Triton X-100 (Sigma)]. After several washes in TBST, embryos were washed with 125 mM Tris-HCl (pH 9.5), 125 mM NaCl, 62.5 mM MgCl₂ and 0.5% Triton X-100, and stained with BMPurple (Roche) at room temperature until the signal was optimal. After the staining, embryos were washed with TBST, fixed in 4% PFA and stored at 4°C. Occasionally, after in situ hybridization, embryos were gelatin embedded and cryosectioned.

RNA antisense probes were synthesized by transcription of linearized DNA from plasmids or from cDNA amplified with specific primers (Table S2). Transcription was carried out with digoxigenin-labeled nucleotides (Roche) and T7 RNA polymerase (Roche). Synthesized RNA was precipitated with 0.8 M ammonium acetate in 75% ethanol or 0.1 M LiCl in 75% ethanol and finally resuspended in 50% formamide-50% RNase free water.

Embryos at E14.5 were eviscerated and fixed in 10% formaldehyde overnight and then washed in acid alcohol (3% HCl in 70% ethanol) several times. Embryos were stained for 3 h with 0.5% w/v Victoria Blue (Sigma) in acid alcohol and after staining embryos were washed in acid alcohol until the embryos were white, then they were washed in 70% ethanol and 95% ethanol. Finally, embryos were clarified with increasing concentrations of methyl salicylate in ethanol (30% and 50%) and stored in 100% methyl salicylate.

Embryos at E18.5 were eviscerated and the skin and soft tissues were removed as much as possible. Embryos were fixed overnight with 95% ethanol and after fixation were submerged in Alcian Blue solution [0.03% w/v Alcian Blue (Sigma), 80% ethanol and 20% glacial acetic acid] overnight. Alcian Blue solution was removed and several washes with 70% ethanol were made during the day with incubation in 95% ethanol overnight. Once the tissue becomes whiter, embryos were cleared with 1% KOH for 3-6 h depending on the stage and the amount of soft tissue that the embryos had. Once cleared, Alizarin Red solution [0.005% Alizarin Red (Sigma) and 1% w/v KOH] was added until the bones were stained. Another clarification step with 1% KOH could be carried out if necessary after staining with Alizarin Red solution, otherwise embryos were transferred to increasing concentrations of glycerol (20% and 50%) and finally placed in 100% glycerol for long-term storage.

Embryos were fixed in 2% PFA, gelatin embedded and cryosectioned. Sections were permeabilized with 0.5% Triton X-100 in PBS for 20 min and blocking was performed with 20% goat serum in PBS for 1 h. The primary antibodies used overnight at 4°C were anti-GFP (Aves laboratories, GFP-1020, 1:200) rabbit anti-RALDH2 (Abcam, ab96060, 1:200) and an anti-MEISa, recognizing C-terminal short isoform of MEIS1 and MEIS2 (Mercader et al., 2005). Secondary antibodies were incubated for 45 min at room temperature. Secondary antibodies used were Alexa-488 (Life Technologies, A11034, 1:500) for anti-MEISa and anti-RALDH2; anti-rabbit-HRP (Dako, P0448, 1:200) for anti-RALDH2 and anti-chicken Alexa-594 (Abcam, ab15172, 1:500) for anti-GFP. Following anti-rabbit-HRP incubation, amplification with Tyr-FITC (Perkin Elmer, NEL74100KT, 1:100) for 3 min at room temperature was performed. Sections were incubated with DAPI and mounted in Vectashield or Dako fluorescent mounting media for acquisition. Images were acquired using a Nikon A1R confocal microscope using 405, 488 and 561 nm wavelengths and Plan Apo 10× DIC L or Plan Apo VC 20× DIC N2 dry objectives.

E9 embryos were eviscerated and the embryonic tissue posterior to the otic vesicle, containing the paraxial mesoderm and the neural tube, was selected and frozen. RNA purification was carried out using the Qiagen miRNeasy Mini Kit and retrotranscribed using Invitrogen SuperScript II Reverse Transcriptase. qPCR was then performed using Sybr Green PCR Master Mix and run in an ABI Prism 7000 from Applied Biosystems. Primers used were as follows: Raldh2 (forward), GAAGGATGGATGCGTCTGAAA; Raldh2 (reverse), TTCCACCAAGTCTGCAAGCTT; Rarb (forward), TATGAGATGACAGCGGAGCTAGAC; Rarb (reverse), GGCTTTCCGGATCTTCTCAGT; Cyp26b1 (forward), AATGAGATTCTGCCGGAGACA; Cyp26b1 (reverse), GAGGCTACACCGTAGCACTCAA; Raldh3 (forward), AGAGTGCGAACCAGTTATGGC; Raldh3 (reverse), ATCTCCTTCTTCCACCTCACATA.

Differential gene expression analysis was carried out between Dll1^Cre;Meis1^flox/flox;Meis2^flox/flox and control embryos at E9. Four embryos were used for each condition and were staged by somite number, choosing the embryos with 20-24 somites. Total RNA was isolated using RNeasy Micro Kit (Qiagen) separating the anterior region, containing the first 10-12 somites (the head was excluded), from the posterior region, with the rest of the somites and the posterior embryonic bud. Total RNA (20 ng) was used to generate barcoded RNA-seq libraries using the NEBNext Ultra RNA Library preparation kit (New England Biolabs). The size and the concentration of the libraries were checked using the TapeStation 2200 DNA 1000 chip. Libraries were sequenced on a HiSeq2500 (Illumina) to generate 60-base single reads. FastQ files for each sample were obtained using bcltofastQ software 2.20.

Sequencing reads were pre-processed by means of a pipeline that used FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to assess read quality, and Cutadapt (Martin, 2011) to trim sequencing reads (eliminating Illumina adaptor remains) and to discard reads that were shorter than 30 bp. The resulting reads were mapped against the mouse transcriptome (GRCm38, release 91; dec2017 archive) and quantified using RSEM v1.2.20 (Li and Dewey, 2011). Data were then processed with a pipeline that used Bioconductor package Limma (Ritchie et al., 2015) for normalization and differential expression analysis, using a blocking strategy to consider gender and developmental stage (number of somites). Genes with at least one count per million in at least four samples (14,731 genes) were considered for further analysis. We considered as differentially expressed those genes with Benjamini-Hochberg adjusted P value <0.05. Fold change and log(ratio) values were calculated to represent gene expression differences between conditions. Ingenuity Pathway Analysis software from Qiagen was used to identify regulatory pathway enrichment.

We thank members of the Torres group for stimulating discussions and suggestions. We thank members of the Genomics, Bioinformatics, Pluripotent Cell Technology, Transgenesis and Animal Facility CNIC units for excellent support. Meis1 floxed mice were generated by Keith Humphries and kindly provided by Hesham Sadek. We thank Miguel Manzanares for providing mouse strains, and Marian Ros, Aimee Zuniga, Paola Bovolenta, Robb Krumlauf and Jose Luis de la Pompa for kindly providing in situ probes.

The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.193813.reviewer-comments.pdf