A. BA organization in mammals. BA3-6 are collectively indicated as PBA. The same colour code (BA2 red, PBA green) is used throughout the manuscript. B. Heatmap of Hox expression in E10.5 mouse BA1, BA2 and PBA, based on the normalized expression values count per million (CPM) [30]. C. Overlap of HOXA3 binding in PBA and HOXA2 binding in BA2 (200 nt summits, overlap at least 1 nt). Only peaks with FE≥10 are considered. D. UCSC tracks (mm9) of HOXA3 (green) and HOXA2 (red) specific and shared peaks. E. Overlap (expressed as %) of increasing numbers of top HOXA2 and HOXA3 peaks (ranked by FE). High-confidence peaks show the smallest overlap. FG. GREAT analysis of HOXA3- (F) and HOXA2- (G) specific peaks (non-overlapping, green and red bars respectively) shows association with genes involved in different biological processes and whose mutations generate different phenotypes in mouse. The length of the bars corresponds to the binomial raw (uncorrected) P-values (x-axis values). H. HOXA2 binding in PBA. Overlap of HOXA2 summit regions in PBA (FE ≥10, green) with HOXA2 summit regions in the BA2 (red) and HOXA3 summit regions in the PBA (green); same rule as in C. HOXA2 binding locations are similar in BA2 and PBA.

A. Homer detects different variants of the HOX-PBX motif in top 250 HOXA2 and HOXA3 peaks, with a G/C (HOXA3) or mainly a G (HOXA2) in the second variable position. B. Occurrence of HOX-PBX motif variants (all permutations of the variable nucleotides in TGATNNAT) in top 250 HOXA2 and HOXA3 peaks (ordered into 50 region bins by decreasing FE). The TGATTCAT motif (red arrows) is among the most enriched variants in HOXA3 peaks but does not virtually occur in HOXA2 peaks. C. Distribution of differential H3K27 acetylation (PBA/BA2 ratio) at top 250 HOXA3 peaks containing the four most frequent TGATNNAT variants. Data is ordered based on fold change of acetylation (fold change of normalized read counts) and divided into quartiles. The central quartiles (50% of the data) are shown in the box, while top and bottom whiskers represent the top 25% and bottom 25% (the top and the bottom quartiles) of the data, respectively. HOXA3 peaks containing the TGATTCAT variant are associated with increased enhancer activity in PBA (red line). D. UCSC tracks with HOXA3, HOXA2, PBX and MEIS binding profiles in BA2 (red) and PBA (green) at the Sulf2 locus, containing TGATTCAT. No HOXA2 binding is detected in BA2 or PBA. E. Sequence of HOXA3 peak summit in D, corresponding to the probe used in F. The TGATTCAT motif (underlined) is flanked by two MEIS motifs (also underlined); the C➔G substitution tested in G is indicated in red. F. HOXA3 can selectively bind the Sulf2 probe in complex with PBX/MEIS. Incubation of the Sulf2 probe with TNT reticulocytes expressing HOXA2, HOXA3, MEIS/PBX, HOXA2/MEIS/PBX or HOXA3/MEIS/PBX. MEIS/PBX bind the Sulf2 probe in combination (arrow). Addition of HOXA3 to the probe results in the formation of a complex only in the presence of PBX/MEIS (arrow). No complex is formed when PBX/MEIS are co-translated with HOXA2. G. Same experiment as in F, using a mutant Sulf2 probe (the nucleotide substitution is shown in E). HOXA2 can bind the mutant probe in combination with MEIS/PBX (asterisk), similar to HOXA3 (arrow).

A. UCSC tracks of HOXA2, HOXA3, PBX, MEIS binding and H3K27 acetylation profiles in BA2 (red) and PBA (green) at the Meis2 locus. Strong HOX and TALE binding is observed in both tissues, with higher acetylation levels in BA2. B. Heatmap shows Meis2 and Zfp703 expression in E11.5 mouse BA1, BA2 and PBA, based on the normalized expression values CPM [30]. C. Meis2 enhancer is active in the hindbrain (h) and the BAs (ba, arrow) of developing zebrafish (72 hours post fertilization), which correspond to Meis2 expression domains in mouse [18]. The enhancer sequence spans the 200nt summit of HOXA2 peak in A. D. Luciferase activity driven by Meis2 enhancer co-transfected with Hoxa2 (red bar) or Hoxa3 (green bar) in combination with Meis2 and Pbx1a expression vectors in NIH3T3 cells. The combination of Hoxa2 with Meis2 and Pbx1a results in the highest activation. Changing the HOX-PBX site (empty bars, mutant sequence in F) reduces HOX-TALE activation. E. Luciferase activity driven by Meis2 enhancer co-transfected with Hoxa2-a3HD (red empty bar) or Hoxa3-a2HD (green empty bar) and Meis2 and Pbx1a. Values shown in DE represent fold activation over basal enhancer activity and are presented as the average of at least two independent experiments, each performed in triplicate. Error bars represent the standard error of the mean (SEM). One-way ANOVA with post-hoc Tukey HSD (Honestly Significant Difference) test: ** = p<0.005; *** = p<0.0005. F-G. Incubation of the Meis2 probe with TNT reticulocyte expressing HOXA2, HOXA3, MEIS/PBX, HOXA2/MEIS/PBX or HOXA3/MEIS/PBX as indicated. F. MEIS and PBX bind the Meis2 probe together (black arrow, see also S3E and S3F Fig). Addition of HOXA2 and HOXA3 results in a trimeric protein complex (red and green arrows respectively). Co-translation of HOXA2 with PBX/MEIS gives rise to a stronger trimeric complex relative to HOXA3 and, accordingly, to a decrease in PBX-MEIS binding as a dimer (black arrow; see also S3G Fig for quantification). G. HOX binding to Meis2 probe, alone and in combination. The total HOX dose is unvaried and corresponds to 100%; HOXA2 and HOXA3 are translated alone and co-translated with each other at different ratios, as follows: 100% HOXA2; 75% HOXA2 + 25% HOXA3; 50% HOXA2 + 50% HOXA3; 25% HOXA2 + 75% HOXA3; 100% HOXA3. Colour coded % are indicated below the gel; HOXA2 is red and HOXA3 is green. HOXA2 forms a stronger trimeric complex (red arrow) than HOXA3 (green arrow) at comparable concentrations. PBX/MEIS dimeric complex, black arrow (see also S3H and S3I Fig for quantification). H. Swapping HOXA3-HD with HOXA2-HD does not improve the ability of HOXA3 to form a ternary complex with PBX and MEIS (green arrows), and does not decrease HOXA2 binding with MEIS and PBX (red arrows). Adding HOXA2 (or HOXA2-A3HD) results in higher intensity of the trimeric complex (red arrows) and lower intensity of TALE dimeric complex relative to HOXA3 (or HOXA3-A2HD), as observed in F. I. Quantification of MEIS/PBX binding to Meis2 mutant probes. The sequence of Meis2 wild-type and mutant probes is shown below: HOX-PBX and MEIS motifs are underlined and nucleotide substitutions are shown in red. Mutants are as follows: m1 and m2 contain mutations in each half of the HOX-PBX site, m1+m2 in both halves; m3 contains mutations in MEIS site closest to HOX-PBX motif and m4 bears mutations in both MEIS sites. MEIS/PBX binding to each mutant probe is expressed as relative percentage to PBX/MEIS binding to Meis2 wild-type in the absence of HOX. Binding of PBX/MEIS dimer was quantified in the absence (black boxes) and in the presence of HOXA2 (red boxes) and HOXA3 (green boxes). Values are the result of three independent experiments; for accurate quantification, each gel was run with the reference condition (PBX/MEIS bound to Meis2 wild-type in the absence of HOX). All mutations significantly decrease PBX/MEIS binding (p<0.0005, one-way ANOVA with post-hoc Tukey HSD test). J. Top HOXA2 and HOXA3 overlapping peaks (total of 60/250 intersecting HOXA2 and HOXA3 top peaks) are more frequently associated with genes with higher expression in BA2 (red) relative to PBA (green). The white portion of the pie chart refers to genes that are not differentially expressed between BA2 and PBA (no DE). Gene association is based on GREAT standard association rules; expression levels are extracted from E11.5 RNA-seq [30].

AB. In situ hybridization on E9.5 embryos, using Hoxa2 (A) and Hoxa3 (B) probes. A. Hoxa2 is highly expressed in the neural crest migrating from rhombomere 4 (asterisk) to the BA2 (arrow). The portion of neural crest migrating just below the otic vesicle (OV) into the BA3 (arrowhead) is also Hoxa2-positive. B. Hoxa3 is expressed in the BA3 (arrowhead). C. Distribution of HOXA2 peaks FE in BA2 and PBA. HOXA2 peaks in BA2 and PBA are ordered based on log2 FE and divided into quartiles. The central quartiles (50% of the data) are shown in the box, while top and bottom whiskers represent the top 25% and bottom 25% (the top and the bottom quartiles) of the data, respectively. HOXA2 peaks display a significantly higher FE in BA2 compared to PBA (p < 2.2e-16; one sided Welsh two sample t-test). D. Comparison of HOXA2 binding in BA2 (red bars) and PBA (green bars) by ChIP-qPCR. Enrichment of each region following immunoprecipitation with HOXA2 and IgG negative control antibody (Neg Ab) is calculated as percentage input; numbers indicate the corresponding FE values in HOXA2 ChIP-seq (BA2 and PBA). Peaks are labelled by their closest genes. Itih4 is a negative control (unbound region). Values represent the average of duplicate samples, and error bars indicate the SEM. E. HOX binding to Meis2 probe, alone and in combination. Co-translation of HOXA2 with PBX/MEIS and with or without HOXA3; the same applies to HOXA3. The highest HOX dose corresponds to 100%; HOXA2 (red) is translated at progressively decreasing levels, either alone, or co-translated with progressively increasing levels of HOXA3 (green) (e.g. 75% HOXA2; 75% HOXA2 + 25% HOXA3; 50% HOXA2; 50%HOXA2 + 50% HOXA3; 50% HOXA3 as indicated below the gel). At 50% and 25% levels, HOXA2 forms a significantly stronger trimeric complex (red arrow) on its own than when HOXA3 is added. PBX/MEIS dimeric complex, black arrow; HOXA3/PBX/MEIS trimeric, green arrow (see also S4C Fig). F. Luciferase activity driven by Meis2 and Zfp703 enhancers co-transfected with expression vector for Hoxa2 or Hoxa3, alone, or at diverse ratio of Hoxa2 to Hoxa3 (3:1; 2:2; 1:3) as indicated. All samples contain Meis2 and Pbx1a expression vectors; Hoxa2 and Hoxa3 are added as indicated. For both enhancers, luciferase activity decreases as Hoxa2 is progressively replaced by Hoxa3 (asterisks indicate significant difference relative to full dose of Hoxa2); for Meis2, enhancer activation is significantly higher with half dose of Hoxa2 alone, relative to Hoxa2 and Hoxa3 together. Values represent fold activation over basal enhancer activity and are presented as the average of at least two independent experiments, each performed in triplicate. Error bars represent the SEM. One-way ANOVA with post-hoc Tukey HSD test: * = p<0.01; ** = p<0.001.

A. Overlap of HOXA3 with MEIS and PBX peaks in the same tissue (PBA) and at the same embryonic stage (E11.5) (200nt summit regions, overlap at least 1nt). The proportional Venn diagram is cropped to focus on HOXA3 peaks. B. Barplots of fold change in MEIS binding levels in PBA versus BA1. Regions co-occupied by MEIS with HOXA3 in PBA generally display higher MEIS binding levels in PBA (HOX-positive) relative to the HOX-negative BA1. In contrast, MEIS binding not overlapping HOXA3 can be higher in BA1 or in PBA. Fold changes in MEIS peaks were calculated using EdgeR (see also S7 Fig). C. Distribution of FE across different groups of MEIS peaks (PBA). MEIS peaks are sorted into: peaks not overlapping HOX (light green), MEIS peaks overlapping HOXA3 only (‘exclusive’ peaks, i.e. not overlapping HOXA2 in PBA, darker green) and MEIS peaks overlapping HOXA2 and HOXA3 (darkest green). Kernel density plot shows that MEIS peaks overlapping HOX display higher FE relative to MEIS only peaks. D. Distance of HOXA3 peaks relative to MEIS peaks (PBA). HOXA3 peaks are binned according to their log₁₀ distance to the nearest MEIS peak and labelled according to FE (high FE, dark red bars; low FE, dark blue bars). High HOXA3 peaks mainly occur within 1kb of a MEIS peak. E. Co-precipitation assays. HEK293 cells were co-transfected with expression vectors for FLAG-tagged HOXA2 or HOXA3 and GST-tagged MEIS1, GST-tagged MEIS2 or GST alone. Protein interactions were assayed by co-precipitation on glutathione beads directed toward the GST tag and eluted proteins analysed by western blotting to detect the presence of FLAG-HOXA2 or FLAG-HOXA3 (red box, Co-P). Cell lysates were analysed by western blotting prior to co-precipitation to detect protein expression (input).

A. Proportional Venn diagram shows highly overlapping binding of MEIS in BA1, BA2 and PBA. MEIS peaks were combined and re-centered using DiffBind; out of 215830 MEIS peaks, 101055 are in common between the three tissues. B. Top MEIS peaks are different in BA2 and PBA. At high FE, there is a lower overlap between MEIS binding in BA2 and PBA. The ratio of MEIS peaks, which are common to BA2 and PBA, increases as FE decreases. C. UCSC tracks illustrate the different MEIS binding levels in BA1-BA2-PBA at the Zfp496 and Zfpm1 loci. Instances of common MEIS peaks higher in one tissue (PBA) are shaded. D. Correlation plot of differential MEIS binding and differential acetylation (enhancer activity) at intergenic regions (PBA versus BA2). Each point corresponds to a region with MEIS log₂ fold change >1 (FC>2); the corresponding H3K27ac value is plotted. Changes in MEIS binding levels are positively correlated with increased enhancer activity in the same tissue (correlation = 0.73). E. Boxplots of MEIS binding and H3K27Ac levels (log₂ fold change) associated with differentially expressed genes. Peaks are associated with intragenic and intergenic (+/-100kb from transcription start and end sites) regions of differentially expressed genes (higher expression in PBA, green; higher expression in BA2, red; non-differentially expressed, grey). Irrespective of location, MEIS binding and acetylation levels are significantly higher (p<2.2e-16, one-sided t-test) when associated with genes with higher expression in the same tissue, relative to MEIS binding and acetylation associated with genes expressed at the same levels in BA2 and PBA. Data is ordered based on fold change of acetylation and MEIS binding (fold change of normalized read counts) and divided into quartiles. The central quartiles (50% of the data) are shown in the box, while top and bottom whiskers represent the top 25% and bottom 25% of the data, respectively. F. Changes in H3K27ac (log₂RPKM) in BA2 and PBA for all HOX peaks and for HOX peaks overlapping MEIS differential binding higher in BA2 (HOXA2 peaks) and higher in PBA (HOXA3 peaks). Data is ordered based on fold change of acetylation (fold change of normalized read counts), as described above. HOX binding generally increases H3K27Ac; peaks associated with increased MEIS binding display a significantly higher increment of H3K27Ac in the same tissue (One-way ANOVA with post-hoc Tukey HSD; ** = p<0.001). G. CNN models of MEIS differential peaks uncover enrichment of tissue-specific sequence motifs as described [35]. MEIS binding was classified in six categories (i.e. peaks with higher/lower binding in BA1, BA2, PBA). CNN analysis identifies tissue-specific sequence features in each class of MEIS peaks. Predicted GATA binding in a MEIS PBA up-binding region is visualised as in the example (a feature matching GATA TF recognition motif on chr5:104257972-104258015 is shown). GATA6 ChIP-seq in PBA verifies this prediction. HOMER was used to cluster and annotate tissue-specific sequence features; features enriched in the different classes of MEIS peaks are matched to TF families with known tissue-specificity (see also Fig 6H). H. Heatmap of the expression of selected TF families, corresponding to cognate recognition motifs identified in MEIS PBA-up, in E11.5 mouse BA2 and PBA. Members of the GATA and TBX families, and the majority of expressed Forkhead TFs are enriched in PBA relative to BA2. Only TFs with expression values > 10 cpm in at least one tissue are shown. I. HOMER de novo motif discovery in HOXA3-specific and HOXA2-specific peaks. HOXA3-specific are HOXA3 peaks excluding peaks overlapping with HOXA2 BA2; similarly, HOXA2-specific are HOXA2 peaks excluding peaks overlapping with HOXA3 PBA. HOMER identifies enrichment of the same motifs enriched in BA-specific MEIS differential binding, Forkhead motif in HOXA3-specific (shaded in green) and BHLH motif in HOXA2-specific subsets (shaded in red). Variations of HD recognition motifs potentially recognized by HOX and attributed by HOMER to PBA-specific TFs NKX and ISL1 in PBA and LHX/DLX in BA2 are also enriched. J. Effects of Forkhead and GATA motifs on HOX and MEIS binding, assessed by in silico knockout. CNN MEIS PBA ‘up-binding’ features (Fig 6G) were annotated as HOX, GATA, and Forkhead (see methods). Co-occurring HOX- Forkhead motifs (distance between 1 nt to 100 nt) were selected for in silico mutagenesis. Forkhead mutagenesis results in a significant drop in HOXA3 binding in PBA, but shows no average significant effect on HOXA2 in BA2. Similarly, Forkhead mutagenesis significantly decreases Meis PBA binding across most tested sites. In comparison, much weaker effects are predicted on BA1 and BA2 MEIS differential binding. As a negative control, the same procedure was applied to co-occurring HOX-GATA motifs. GATA motif mutagenesis does not show significant average effects on HOX, or MEIS in HOX-bound regions. K. Luciferase activity driven by Sfrp2 enhancer co-transfected with Meis alone (grey bar) and with Hoxa2 (red empty bars), Hoxa3 (green empty bars); Meis and Pbx with Hoxa2 (red bars) and Hoxa3 (green bars) in the absence (left half of the graph) and in the presence of Foxc1 (right half of the graph, as indicated) in NIH3T3 cells. The last two samples contain Sfrp2 enhancer co-transfected with Hoxa2 (red empty bars), Hoxa3 (green empty bars). Adding Foxc1 to Hoxa2 or Hoxa3 with Meis2 and Pbx1a results in the highest activation.

Low affinity, widespread binding of MEIS (blue square) defines a large subset of accessible chromatin (grey line) for activation (PBX is not shown as PBX and MEIS binding almost entirely overlap). Direct cooperativity with HOX (A2 and A3, red and green circles respectively) and/or indirect cooperativity with tissue-specific TFs (triangle) increase MEIS binding affinity and residence time; prolonged residence time of MEIS at enhancers promotes recruitment of general co-activators (yellow) and activation of transcription. HOX paralogs preferentially bind different subsets of MEIS occupied regions, resulting in differential transcription. Three examples of BA-specific transcription are shown. In a, the red site is bound with higher affinity by HOXA2 than HOXA3, resulting in the formation of a more stable HOX-TALE complex on DNA and a (higher) transcriptional output in BA2. Conversely, in c, the green site is only recognized by HOXA3, leading to high affinity MEIS binding only in PBA, and to PBA-specific transcription. In b, the effect of HOXA3 is potentiated by a PBA-specific TF binding in the vicinity. Co-binding with tissue-specific TFs may positively contribute to HOX-MEIS cooperativity by competing with nucleosome for DNA binding, especially at HOX and/or MEIS low affinity sites. These mechanisms result in BA-specific transcription.