Discovery of genes required for body axis and limb formation by global identification of retinoic acid–regulated epigenetic marks

Identification of target genes that mediate required functions downstream of transcription factors is hampered by the large number of genes whose expression changes when the factor is removed from a specific tissue and the numerous binding sites for the factor in the genome. Retinoic acid (RA) regulates transcription via RA receptors bound to RA response elements (RAREs) of which there are thousands in vertebrate genomes. Here, we combined chromatin immunoprecipitation sequencing (ChIP-seq) for epigenetic marks and RNA-seq on trunk tissue from wild-type and Aldh1a2-/- embryos lacking RA synthesis that exhibit body axis and forelimb defects. We identified a relatively small number of genes with altered expression when RA is missing that also have nearby RA-regulated deposition of histone H3 K27 acetylation (H3K27ac) (gene activation mark) or histone H3 K27 trimethylation (H3K27me3) (gene repression mark) associated with conserved RAREs, suggesting these genes function downstream of RA. RA-regulated epigenetic marks were identified near RA target genes already known to be required for body axis and limb formation, thus validating our approach; plus, many other candidate RA target genes were found. Nuclear receptor 2f1 (Nr2f1) and nuclear receptor 2f2 (Nr2f2) in addition to Meis homeobox 1 (Meis1) and Meis homeobox 2 (Meis2) gene family members were identified by our approach, and double knockouts of each family demonstrated previously unknown requirements for body axis and/or limb formation. A similar epigenetic approach can be used to determine the target genes for any transcriptional regulator for which a knockout is available.


Introduction
Retinoic acid (RA) is generated from retinol by the sequential activities of retinol dehydrogenase 10 (RDH10) [1] and aldehyde dehydrogenase 1A2 (ALDH1A2) [2,3]. Knockout studies of these enzymes revealed an essential role for RA in many early developmental programs, including those controlling hindbrain anteroposterior patterning, neuromesodermal progenitor (NMP) differentiation, spinal cord neurogenesis, somitogenesis, forelimb bud initiation, and heart anteroposterior patterning [4,5]. RA functions as a ligand for nuclear RA receptors (RARs) that bind DNA sequences known as RA response elements (RAREs) as a heterodimer complex with retinoid X receptors (RXRs) [6]. Binding of RA to RAR alters the ability of RAREs to recruit nuclear receptor coactivators (NCOAs) that activate transcription or nuclear receptor corepressors (NCORs) that repress transcription [7]. Thus, RA functions are mediated by transcriptional activation or repression of key genes via RAREs. Identification of RA-regulated genes that are required for development has been difficult, as loss or gain of RA activity alters the mRNA levels of thousands of genes in various cell lines or animals, perhaps most being indirect targets of RA or regulated posttranscriptionally. As RA target genes are dependent upon RAREs, identification of RAREs by RAR-binding studies, cell line transfection assays, and enhancer reporter transgenes in mouse or zebrafish have been used to identify RA target genes that may be required for development, but progress is slow, as each gene is analyzed separately [5]. Genomic RAR chromatin immunoprecipitation sequencing (ChIP-seq) studies on mouse embryoid bodies and F9 embryonal carcinoma cells reported approximately 14,000 potential RAREs in the mouse genome [8,9], but it is unclear how many of these RAREs are required to regulate genes in any specific tissue, and many may not function in any tissue at any stage of development. Only a few RAREs have been shown to result in gene expression and developmental defects when subjected to deletion analysis in mouse, i.e., a RARE enhancer that activates Hoxa1 in the hindbrain [10], a RARE enhancer that activates Cdx1 in the spinal cord [11], and a RARE that functions as a silencer to repress caudal Fgf8 in the developing trunk [7]. In 1 additional case, a RARE described within intron 2 of Tbx5 that was suggested to be required for activation of Tbx5 in the forelimb field based on a mouse enhancer reporter transgene [12] was found to be unnecessary for Tbx5 activation and forelimb budding when subjected to CRISPR deletion analysis, suggesting Tbx5 is not an RA target gene [13]. Many DNA control elements (including RAREs) that exhibit appropriate tissue-specific expression in enhancer reporter transgene assays have been shown to not be required as an enhancer in vivo when deleted; this may be due to enhancer redundancy or because the control element is really not an enhancer but appeared to be when inserted as a transgene at a random location in the genome near a heterologous promoter [14]. Thus, additional methods are needed (preferably genome-wide) to locate functional RAREs in a particular tissue that can be used to identify new candidate RA target genes that are required for development.
Epigenetic studies have found that histone H3 K27 acetylation (H3K27ac) associates with gene activation and histone H3 K27 trimethylation (H3K27me3) associates with gene repression [15,16]. We suggest that genes possessing nearby H3K27ac and H3K27me3 marks that are altered by loss of RA may point to direct transcriptional targets of RA (either activated or repressed) that are excellent candidates for performing functions downstream of RA. Here, we performed genomic ChIP-seq (H3K27ac and H3K27me3) and RNA-seq studies on embryonic day (E)8.5 mouse embryonic trunks from wild-type and Aldh1a2-/-mouse embryos lacking RA synthesis to globally identify RA target genes for embryonic trunk. Candidate targets are defined as genes whose mRNA levels are decreased or increased by genetic loss of RA that also have nearby RA-regulated epigenetic marks associated with conserved RAREs, suggesting they have important downstream functions. Our approach was able to identify many previously reported RA target genes known to control embryonic trunk development (including all 3 known RA target genes from RARE knockout studies: Hoxa1, Cdx1, and Fgf8); plus, we identified numerous new candidate RA target genes that may control trunk development. CRISPR knockout studies on several of these new candidate RA target genes validated them as being required for body axis and/or limb formation. Our approach is generally applicable to determine tissue-specific target genes for any transcriptional regulator that has a knockout available.

Comparison of RNA-seq and H3K27ac/H3K27me3 ChIP-seq for Aldh1a2-/-trunk tissue
Embryonic trunks were obtained from E8.5 embryos dissected to remove the head (including the pharyngeal region and anterior hindbrain), heart, and caudal tissue below the most recently formed somite as previously described [17]. We performed RNA-seq analysis comparing E8.5 trunk tissue from wild-type embryos and Aldh1a2-/-embryos that lack the ability to produce RA [3]. This analysis identified 4,298 genes whose mRNA levels in trunk tissue are significantly decreased or increased when RA is absent (fragments per kilobase of transcript per million mapped reads [FPKMs] > 0.5; a cutoff of log2 <−0.85 or >0.85 was employed to include Sox2 known to be activated by RA; data available at GEO under accession number GSE131584).
We performed ChIP-seq analysis for H3K27ac and H3K27me3 epigenetic marks comparing E8.5 trunk tissue from wild-type and Aldh1a2-/-embryos isolated as described previously for RNA-seq. This analysis identified 314 RA-regulated ChIP-seq peaks for H3K27ac located within or near 214 genes (i.e., the genes with the nearest annotated promoters) using a log2 cutoff of <0.51 or >0.51 to include an RA-regulated peak near Sox2 known to be activated by RA [18,19]. We identified 262 RA-regulated peaks for H3K27me3 located within or near 141 genes (i.e., the genes with nearest annotated promoters) using a log2 cutoff of <−0.47 or >0.47 to include an RA-regulated peak near Fst known to be repressed by RA [20]; all ChIPseq data are available at GEO under accession number GSE131624. Thus, we found a much smaller number of RA-regulated ChIP-seq peaks for H3K27ac/H3K27me3 compared with the very large number of genes found to have altered mRNA levels with RNA-seq.
In order to identify genes that are good candidates for being transcriptionally activated or repressed by RA (RA target genes), we compared our ChIP-seq and RNA-seq results to identify RA-regulated ChIP-seq peaks in which nearby genes have significant changes in expression in wild type versus Aldh1a2-/-based on RNA-seq. We found 73 RA-regulated peaks for H3K27ac near 63 genes with significant changes in expression when RA is lost (S1 Table), plus 46 RA-regulated peaks for H3K27me3 near 41 genes with significant changes in expression when RA is lost (S2 Table). As some genes have more than 1 nearby RA-regulated peak for H3K27ac or H3K27me3, and some genes have nearby RA-regulated peaks for both H3K27ac and H3K27me3 (Rarb, Dhrs3, Fgf8, Cdx2, Fst, Meis homeobox 1 [Meis1], Meis homeobox 2 [Meis2], nuclear receptor 2f2 [Nr2f2], Foxp4, Ptprs, and Zfhx4), a total of 93 RA-regulated genes have nearby RA-regulated peaks for H3K27ac and/or H3K27me3 when RA is lost, thus identifying them as candidate RA target genes for trunk development (S1 and S2 Tables, S1A  Fig).
Among the 93 candidate RA target genes for trunk development identified with our approach are many examples of genes previously reported to be regulated by RA in the trunk based on studies of Aldh1a2-/-embryos [5,21] or RA-treated NMPs [20]; this includes Hoxa1, Cdx1, Rarb, Crabp2, Sox2, Dhrs3, and Pax6, whose expression is increased by RA, plus Fgf8, Cdx2, and Fst, whose expression is decreased by RA (Table 1). H3K27ac peaks near Cdx1, Rarb, Crabp2, Sox2, Dhrs3, and Pax6 are reduced in Aldh1a2-/-trunk, consistent with these being RA-activated genes, whereas H3K27ac peaks near Fgf8, Cdx2, and Fst are increased in Aldh1a2-/-, consistent with these being genes repressed by RA. Conversely, H3K27me3 peaks near Fgf8, Cdx2, and Fst are decreased in Aldh1a2-/-, whereas H3K27me3 peaks near Rarb, Hoxa1, and Dhrs3 are increased in Aldh1a2-/-, consistent with the former being genes repressed by RA and the latter being genes activated by RA (Table 1). In addition to these 10 well-established RA target genes that are required for trunk development, we identified 83 Identification of RAREs associated with RA-regulated deposition of H3K27ac or H3K27me epigenetic marks As RA target genes need to be associated with a RARE, the DNA sequences within the RA-regulated H3K27ac/H3K27me3 ChIP-seq peaks we found near our list of 93 RA-regulated genes were searched for RARE sequences using the Homer transcription factor binding site program for the mm10 genome; we searched for 3 types of RAREs including those with a 6-bp direct repeat separated by either 5 bp (DR5), 2 bp (DR2), or 1 bp (DR1) [5], and the presence or absence of RAREs is summarized (S1 and S2 Tables). We found that 46 of these 93 genes contained at least 1 RARE in their nearby RA-regulated H3K27ac and/or H3K27me3 ChIP-seq peaks, thus narrowing down our list of candidate RA target genes to 49% of the genes originally identified. Our approach identified the 3 RAREs previously shown to have required functions during trunk development in vivo by knockout studies (RAREs for Hoxa1, Cdx1, Fgf8) plus several RAREs associated with known RA-regulated genes in the E8.5 trunk from Aldh1a2-/-studies (Rarb, Crabp2, Sox2, Dhrs3, Cdx2, Fst), thus validating our approach for identifying RA-regulated genes required for trunk development. The sequences of all the RAREs near these 46 RA target genes here are summarized; included are 65 RAREs near 34 RA-activated genes (we refer to these as RARE enhancers associated with increased H3K27ac and/or decreased H3K27me3 in the presence of RA) and 20 RAREs near 12 RA-repressed genes (we refer to these as RARE silencers associated with increased H3K27me3 and/or decreased H3K27ac in the presence of RA) (S3 Table).
The results here provide evidence that many of the RA-regulated H3K27ac and H3K27me3 marks are associated with regulation of the nearest genes. However, it is possible that some H3K27ac and H3K27me3 RA-regulated peaks may be related to RA-regulated genes located further away in the same topologically associated domain (TAD). In order to address this issue, we assigned each RA-regulated H3K27ac and H3K27me3 peak to a TAD using the 3D Genome Browser (http://promoter.bx.psu.edu/hi-c/view.php); TAD analysis has not been performed on mouse E8.5 trunk tissue, but as TAD domains are similar between different tissues [22], we used the TAD database for mouse embryonic stem (ES) cells, which is the closest biologically relevant database in the 3D Genome Browser. Then, the genes in each TAD containing an RA-regulated peak were searched in our RNA-seq database to identify genes whose mRNA levels are decreased or increased when RA is lost, and if at least 1 gene was found, we determined whether a RARE is present in the ChIP-seq peak. This analysis resulted in the identification of 82 additional RARE enhancers near RA-activated genes and 40 additional RARE silencers near RA-repressed genes, where the gene is not the gene nearest to the RARE in the TAD; in some cases, more than 1 RA-regulated gene was identified in a TAD (S3 Table).
Up to now, Fgf8 represents the only example of a gene that is directly repressed by RA at the transcriptional level as shown by developmental defects upon knockout of the RARE at −4.1 kb and by the ability of this RARE to stimulate binding of NCOR and polycomb repressive complex 2 (PRC2) plus deposition of H3K27me3 in an RA-dependent manner [7,17].
Here, in addition to Fgf8, we found many more candidates for genes repressed by RA in the trunk based on identification of nearby RARE silencers (S3 Table).

Analysis of known RA target genes for trunk validates our approach for finding new targets
The RA-regulated H3K27ac and/or H3K27me3 peaks we identified near Rarb, Crabp2, Hoxa1, and Cdx1 all overlap previously reported RAREs for these genes (Fig 1). In the case of Rarb, the DR5 RARE in the 5 0 -untranslated region [23] overlaps RA-regulated peaks for both H3K27ac and H3K27me3, suggesting that this RARE in the presence of RA stimulates deposition of H3K27ac and removal of H3K27me3 during activation of Rarb; we identified a DR1 RARE in the 5 0 -noncoding region of Rarb within another RA-regulated H3K27me3 ChIP-seq peak ( Fig 1A). For Crabp2, 2 closely spaced RAREs previously reported in the 5 0 -noncoding  [23]; here, H3K27ac is decreased and H3K27me3 increased near the native RARE when RA is lost in trunk tissue, supporting its function as a RARE enhancer in vivo. We also found a RARE in the 5 0 -noncoding region of Rarb within an H3K27me3 ChIP-seq peak that is increased when RA is lost. (B) RA-regulated peaks for H3K27ac and RAREs are shown for Crabp2. The 2 RAREs in the 5 0 -noncoding region were previously shown to function as RA-dependent enhancers in cell line studies [24]. Our epigenetic studies also identified another RARE enhancer in the 3 0 -noncoding region. (C) RA-regulated peaks for H3K27ac and/or H3K27me3 and RAREs are shown for Hoxa1. KO studies in mouse embryos have shown that the RARE in the 3 0 -noncoding region is essential for hindbrain Hoxa1 expression and development [10]. (D) RA-regulated peaks for H3K27ac and H3K27me3 and RAREs are shown for Cdx1. KO studies in mouse embryos have shown that the RARE in the 5 0 -noncoding region is essential for Cdx1 expression and body axis development [11]. RA-regulated peaks in the genome browser view shown here and elsewhere are for 1 replicate, with the other replicate showing a similar result. ChIP-seq, chromatin immunoprecipitation sequencing; E, embryonic day; H3K27ac, histone H3 K27 acetylation; H3K27me3, histone H3 K27 trimethylation; KO, knockout; RA, retinoic acid; RARE, RA response element; WT, wild-type.
https://doi.org/10.1371/journal.pbio.3000719.g001 region [24] associate with RA-regulated peaks for H3K27ac; plus, another RARE we identified in the 3 0 -noncoding region associates with changes in H3K27ac ( Fig 1B). For Hoxa1, the RARE located in the 3 0 -noncoding region is associated with RA-regulated peaks for both H3K27ac and H3K27me3; plus, another RARE we identified in the 3 0 -untranslated region is associated with RA-regulated peaks for H3K27me3 ( Fig 1C); importantly, knockout studies on the Hoxa1 RARE in the 3 0 -noncoding region demonstrated that it is required in vivo for Hoxa1 expression and normal development [10]. For Cdx1, 2 RAREs have been reported, 1 in the 5 0 -noncoding region that was shown by knockout studies to be required for Cdx1 expression and body axis development [11] plus another RARE in intron 1 [25]. Both of these Cdx1 RAREs are overlapped by RA-regulated peaks for both H3K27ac and H3K27me3 (Fig 1D). These findings demonstrate that our approach can identify genes that are already known to be transcriptionally activated by RA via a RARE and required for development.

Identification of RA-regulated epigenetic marks and RAREs near RAregulated genes known to control NMPs
Ingenuity Pathway Analysis (IPA) of our list of 93 RA target genes shows enrichment for the pathway "development of body trunk", including Sox2, Cdx2, and Fgf8 known to be required for NMP function during trunk development (S1B Fig). NMPs are bipotential progenitor cells in the caudal region coexpressing Sox2 and T/Bra that undergo balanced differentiation to either spinal cord neuroectoderm or presomitic mesoderm to generate the postcranial body axis [26][27][28][29][30][31][32][33]. NMPs are first observed in mouse embryos at about E8.0 near the node and caudal lateral epiblast lying on each side of the primitive streak [34][35][36]. Caudal Wnt and fibroblast growth factor (FGF) signals are required to establish and maintain NMPs [34,[36][37][38][39][40][41]. Cdx2 is required for establishment of NMPs [33]. During development, RA is first produced at E7.5 in presomitic mesoderm expressing Aldh1a2 to generate an anteroposterior gradient of RA with high activity in the trunk and low activity caudally [5]. Loss of RA does not prevent establishment or maintenance of NMPs but does result in unbalanced differentiation of NMPs, with decreased caudal Sox2 expression and decreased appearance of neural progenitors, plus increased caudal Fgf8 expression and increased appearance of mesodermal progenitors and small somites because of the encroachment of caudal Fgf8 expression into the trunk where it reduces epithelial condensation of presomitic mesoderm needed to form somites [19,20,42,43]. Cdx2 expression is increased when RA is lost in Aldh1a2-/-embryos [44].
Here, when RA is lost, we observed RA-regulated H3K27ac and/or H3K27me3 peaks near several genes required for NMP function that show decreased (Sox2) or increased (Fgf8 and Cdx2) expression (Fig 2A-2C). Most of these RA-regulated peaks contain RAREs, providing evidence that Sox2, Fgf8, and Cdx2 are direct RA target genes (S3 Table). For Sox2, we observed 2 RA-regulated H3K27ac ChIP-seq peaks, but only the one in the 3 0 -noncoding region was found to have a RARE (Fig 2A). In the case of Fgf8, previous studies reporting knockout of the RARE located in the 5 0 -noncoding region at −4.1 kb resulted in increased caudal Fgf8 expression and a small somite phenotype (although the defect is not as severe as for Aldh1a2-/-embryos), demonstrating that this RARE functions in vivo as a silencer by RAdependent recruitment of NCORs [7]; RARE redundancy may explain the milder phenotype, as our approach suggests that Fgf8 has 2 additional candidate RARE silencers ( Fig 2B). RARE redundancy may be common, as we observe that Cdx2 has 2 candidate RARE silencers ( Fig  2C), and our overall analysis shows that many genes have more than 1 nearby RARE (S3 Table). These findings provide evidence that RA controls NMP differentiation directly at the transcriptional level by activating Sox2 and repressing Fgf8 and Cdx2 as progenitor cells progress from a caudal to a trunk location. . A RARE (green) was found in the 3 0 -noncoding peak (but not the 5 0 -noncoding peak), suggesting it may function as a RARE enhancer as the H3K27ac peak is decreased when RA is lost. (B) Shown are RA-regulated ChIP-seq peaks for H3K27me3 and H3K27ac near Fgf8. In the 5 0 -noncoding region of Fgf8, we found 2 RAREs on either end of the peak for H3K27me3 (repressive mark) that is decreased in KO, indicating they are candidate RARE silencers; the RARE furthest upstream in the 5 0 -noncoding region at −4.1 kb was shown by knockout studies to function as an RA-dependent RARE silencer required for caudal Fgf8 repression and somitogenesis [7]. We also found another RARE in the 3 0 -noncoding region of Fgf8 that is another candidate for a RARE silencer, as it is contained within an RAregulated peak for H3K27ac (activating mark) that is increased when RA is lost. (C) Cdx2 has a peak for H3K27ac that is increased

Evidence for genes regulated indirectly by RA at the transcriptional level
Our studies show that many genes that are down-regulated or up-regulated following loss of RA are associated with RA-regulated peaks for H3K27ac or H3K27me3 (either nearby or in the same TAD) that do not contain RAREs (S1 and S2 Tables). Such genes may be indirectly activated or repressed by RA at the transcriptional level. In the case of Pax6, our results indicate that RA stimulates H3K27ac deposition in Pax6 introns 2 and 6 that do not contain RAREs, with no other RA-regulated peaks in the same TAD (Fig 3A). Previous studies identified an enhancer in Pax6 intron 6 containing a SOXB1 binding site that is important for activation in the spinal cord [45]. Activation of Pax6 in the spinal cord requires CDX proteins in the posterior-most neural tube, and CDX binding sites have been identified in Pax6 intron 2 [46]; in addition to expression in the caudal progenitor zone, mouse Cdx1 is expressed in the posterior neural plate where Pax6 is activated, and this expression domain requires RA [44]. Activation of Pax6 requires that caudal FGF signaling be down-regulated [43]. Thus, although it is possible that our H3K27ac/H3K27me3 studies failed to identify an unknown RARE near Pax6, our findings suggest that the RA requirement for Pax6 activation may operate through several indirect mechanisms because of the ability of RA to activate Sox2 and Cdx1 and repress Fgf8 (Figs 1 and 2).
We observed that Spry4 (shown here to be down-regulated by RA) does not have a RARE associated with its RA-regulated ChIP-seq peak for H3K27me3; no other RA-regulated peaks were found in its TAD (Fig 3B). Many of the RA-regulated ChIP-seq peaks observed with our approach that do not contain RAREs may be indirect RA-regulated peaks that contain DNA binding sites for transcription factors other than RARs whose expression or activity is altered by loss of RA, thus resulting in changes for H3K27ac/H3K27me3 marks that are caused by the other transcription factors.

Conservation of RAREs identified with our approach identifies candidate RA target genes
Although it possible that some RAREs that are conserved only in mammals perform mammalspecific functions, RAREs that are conserved from mammals to birds or lower may play fundamental roles in regulation of target genes. The candidate RARE enhancers and RARE silencers we identified here that are associated with RA-regulated epigenetic marks were searched for evolutionary conservation using the University of California Santa Cruz (UCSC) genome browser. Among the RAREs in which the nearest gene is RA-regulated, we found 6 RAREs that are conserved from mouse to zebrafish, 11 conserved to frog (Xenopus tropicalis), 18 conserved to reptile (lizard, painted turtle), 20 conserved to bird (chicken, turkey), 39 conserved to human, 65 conserved to rodent (rat), and 20 that are not conserved with rat (S3 Table). The large number of RAREs (i.e., 20) conserved beyond mammals to bird, lizard, frog, or fish demonstrate that our approach is able to identify highly conserved RAREs that point to excellent candidate genes required for development. Among the additional RAREs we found located farther away in the TAD from an RA-regulated gene, we identified only 4 more RAREs conserved beyond mammals to bird, lizard, frog, or fish, thus bringing the total to 24 highly conserved RAREs (S3 Table). Thus, most of the highly conserved RAREs we identified are located and an overlapping peak for H3K27me3 that is decreased, along with 3 RAREs included within both peaks, indicating that all these RAREs are candidates for RARE silencers. ChIP-seq, chromatin immunoprecipitation sequencing; E, embryonic day; H3K27ac, histone H3 K27 acetylation; H3K27me3, histone H3 K27 trimethylation; KO,; NMP, neuromesodermal progenitor; RA, retinoic acid; RARE, RA response element; WT, wild-type.
https://doi.org/10.1371/journal.pbio.3000719.g002 very close to an RA-regulated gene rather than further distant in the TAD. In addition, all these highly conserved RAREs are either identical to the RARE consensus or have only 1 mismatch. Here we summarize the 24 most highly conserved RAREs that point to 38 RA-regulated genes that may be required for development ( Table 2).
As RAREs need to be bound by an RAR in order to function, we examined previously reported RAR ChIP-seq databases for mouse embryoid bodies [8] and mouse F9 embryonal carcinoma cells [9] to determine whether the highly conserved RAREs we identified are included in RAR-binding regions. We found that 19 of our 24 highly conserved RAREs are included in the RAR ChIP-seq peaks from at least 1 of those studies (S4 Table).
Our list of best candidate RA target genes ( Table 2) includes several for which knockout studies have already demonstrated required functions during trunk development, i.e., in RA signaling (Rarb, Dhrs3), body axis formation (Hoxa1, Hoxa4, Hoxa9, Sox2, Fgf8, Pbx1, Tshz1, Zbtb16), and foregut formation (Foxp4); mouse knockout data are summarized by Mouse Genome Informatics (http://www.informatics.jax.org). This list also includes many genes for which knockout studies have either not been performed or knockouts resulted in no reported early developmental defects. This list of genes thus contains excellent new candidates that can be tested for function during trunk development by generating knockouts or double knockouts in the case of gene families.

Nr2f and Meis gene families have nearby RA-regulated epigenetic marks associated with highly conserved RARE enhancers
We identified 2 gene families (Nr2f and Meis) in which 2 family members have decreased expression when RA is lost and nearby RA-regulated peaks for H3K27ac or H3K27me3 containing RAREs. Nr2f1 and Nr2f2 encode orphan nuclear receptors NR2F1 and NR2F2 (previously known as COUP-TFI and COUP-TFII, respectively) that regulate transcription, although they have not been found to have endogenous ligands that control their activity [47].  Meis1 and Meis2 encode transcription factors belonging to the three amino acid loop extension (TALE) family of homeodomain-containing proteins [48,49]. Previous studies suggested that Nr2f genes are activated by RA in Ciona, zebrafish, and mouse cell lines [50][51][52][53]. Here, Nr2f1 and Nr2f2 were both found to have a single RARE in the 5 0 -noncoding region close to exon 1 that is overlapped by or close to the edge of RA-regulated H3K27ac and H3K27me3 peaks (Fig 4A and 4B). Recent studies in zebrafish identified RAREs in similar locations in the nr2f1a and nr2f2 genes [52], and this conservation to mouse was detected by our analysis (Table 2).
Meis1 and Meis2 were previously shown to be up-regulated by RA in chick limbs treated with RA [54], and loss of RA in Aldh1a2-/-embryos results in reduced expression of Meis2 in paraxial mesoderm [55]. Other studies have shown that Meis1 and Meis2 are activated by RA in embryonic stem cells and other cell lines, and RAREs were identified in their 5 0 -noncoding regions [56,57]. Here, Meis1 was found to have 4 RAREs in introns 1, 6, and 7 that are overlapped by RA-regulated peaks for H3K27ac and/or H3K27me3; plus, we identified the previously reported RARE in the 5 0 -noncoding region that is located at the edge of a small RAregulated H3K27ac peak (Fig 4C). Meis2 was found to have 2 RAREs that are overlapped by RA-regulated peaks for H3K27ac and/or H3K27me3, 1 in the 5 0 -noncoding region (previously identified) and another in intron 7 (Fig 4D). Our analysis shows that Meis1 and Meis2 each have a highly conserved DR5 RARE enhancer (Table 2). Together, these studies identify Nf2f1, Nr2f2, Meis1, and Meis2 as candidate RA target genes in the developing trunk.

Nr2f1 and Nr2f2 function redundantly to control body axis formation
In order to be a biologically important direct RA target gene, the gene must be controlled by a RARE and must also perform a function downstream of RA during trunk development, which can be determined by gene knockout studies. Here, we sought to validate our approach by performing knockout studies on some of the candidate RA target genes, particularly those that have nearby highly conserved RAREs. One could also undertake deletion studies of the RAREs, but this is only relevant after a knockout of the gene itself shows a defect. In addition, as genes are often controlled by redundant enhancers (which we observed here for many genes that have 2 or more RAREs associated with RA-regulated epigenetic marks; S3 Table), studies in which predicted enhancers are deleted often have no effect on development [13,14,[58][59][60]; this includes knockout studies we performed on a RARE that was predicted by enhancer transgene studies to be needed for Tbx5 expression in forelimb bud that had no effect on Tbx5 or development [13]. Next, we describe gene knockout studies on candidate RA target genes with nearby highly conserved RAREs to determine whether these genes have a required function in trunk development.
Nr2f1 and Nr2f2 were selected for gene knockout because they both have nearby candidate RARE enhancers (identified by our H3K27ac/H3K27me3 ChIP-seq analysis) that are conserved from mouse to zebrafish (Table 2). Nr2f1 (formerly known as COUP-TFI) and Nr2f2 (formerly known as COUP-TFII) are both expressed at E8.5 in somites and presomitic mesoderm but not the spinal cord, suggesting they may function in mesoderm formation during body axis formation [61,62]. Here, in situ hybridization analysis shows that Nr2f1 and Nr2f2 have reduced expression in the trunk of E8.5 Aldh1a2-/-embryos compared with wild type (S3 Fig).
The Nr2f1 knockout is lethal at birth with brain defects, but no somite, spinal cord, or body axis defects are observed [63]. The Nr2f2 knockout is lethal at E10.5 with defects in heart development but not body axis formation [64]. As redundancy may have masked a body axis defect, we generated Nr2f1/Nr2f2 double mutants. As it would be quite time-consuming and  expensive to obtain (if possible) the previously described Nr2f1 and Nr2f2 single-knockout mouse lines, then generate a double heterozygote mouse line, and then generate double homozygote embryos at a ratio of 1:16, we employed CRISPR/Cas9 gene editing to examine F0 embryos as we previously described for Ncor1/Ncor2 double mutants [7]. Fertilized mouse oocytes were injected with single-guide RNAs (sgRNAs) designed to generate frameshift knockout deletions in the second exons of both Nr2f1 and Nr2f2. After dissecting F0 embryos at E9.0, we obtained Nr2f1/Nr2f2 double knockouts that exhibited a body axis growth defect, more similar in size to that of wild-type E8.25 embryos (Fig 5). Genotyping showed that embryos carrying 1 or 2 knockout alleles were normal in size compared with E9.0 wild type (Fig 5A), whereas embryos carrying either 3 or 4 knockout alleles exhibited a defect in body axis extension and are similar in size to E8.25 wild type; n = 7 (Fig 5B and 5C). Staining for Uncx somite expression demonstrated that embryos with 1-2 knockout alleles all have a normal number of somites with normal size (Fig 5A), whereas embryos with 3-4 knockout alleles all have fewer somites that are smaller in size; embryos with 3 knockout alleles (either Nr2f1het/Nr2f2-hom or Nr2f1-hom/Nr2f2-het) or 4 knockout alleles (Nr2f1-hom/Nr2f2-hom) have a similar small somite defect (Fig 5B and 5C). As E9.0 Nr2f1/Nr2f2 mutants carrying 3-4 knockout alleles are more similar in size to E8.25 wild type, in order to estimate somite size along the anteroposterior axis, we compared them to Uncx-stained E8.25 wild-type embryos (Fig 5D), thus revealing that the E9.0 mutants have somites about 57% the size of somites in E8.25 wild-type embryos, showing they have a specific defect in trunk development rather than a global body growth defect (Fig 5E).
Overall, our findings show that loss of 3 or 4 alleles of Nr2f1 and Nr2f2 hinders body axis formation and results in smaller somites, thus validating our approach for finding new genes regulated by RA that are required for body axis formation. Nr2f1/Nr2f2 double mutants exhibit a small somite phenotype (57% of normal) that is similar to RA-deficient Aldh1a2-/-embryos (50% of normal); however, they cease body axis extension sooner, as E9.0 Nr2f1/Nr2f2 double mutants have fewer than 10 somites (Fig 5B and 5C), whereas most E9.0 Aldh1a2-/-embryos reach 15 somites [65]. Thus, the more severe growth defect we observe for Nr2f1/Nr2f2 double mutants may be caused at least in part by a cardiac defect [64]. However, as our RNA-seq and ChIP-seq studies were performed on trunk tissue in which the heart was excluded, our findings show that RA-regulated epigenetic marks are observed near Nr2f1/Nr2f1 genes located outside the heart in the trunk. Also, our results show that double mutants have smaller somites than wild-type embryos of a comparable size, revealing a specific effect on somitogenesis (body axis extension) rather than just an overall effect on embryonic growth. In the future, more detailed studies of Nr2f1/Nr2f2 double mutants can be performed to determine how these genes control body axis extension. In addition, future studies can be performed to determine how RARE enhancers function along with other factors to control Nr2f1 and Nr2f2 expression during body axis formation.

Meis1 and Meis2 function redundantly to control both body axis and limb formation
Meis1 and Meis2 were selected for gene knockout, as Meis1 has a nearby candidate RARE enhancer conserved from mouse to frog, and Meis2 has a nearby candidate RARE enhancer conserved from mouse to bird (Table 2). Meis1 and Meis2 are both expressed throughout the trunk and in the proximal regions of limb buds [54]. Here, in situ hybridization analysis shows that Meis1 and Meis2 have reduced expression in the trunk of E8.5 Aldh1a2-/-embryos compared with wild type (S3 Fig). The Meis1 knockout is lethal at E11.5 with hematopoietic defects, but no body axis or limb defects are observed [66]. The Meis2 knockout is lethal at E14.5 with defects in cranial and cardiac neural crest, but no defects in body axis or limb formation were observed [67]. As redundancy may have masked a body axis or limb defect, we generated Meis1/Meis2 double mutants via CRISPR/Cas9 gene editing of fertilized mouse oocytes employing sgRNAs designed to generate frameshift knockout deletions in the second exons of both Meis1 and Meis2. Embryos were dissected at E10.5 and stained for somite Uncx expression. Genotyping showed that E10.5 embryos carrying 1 or 2 knockout alleles for Meis1/Meis2 were normal in size with normal-size somites compared with E10.5 wild type (Fig 6A). However, E10.5 embryos carrying 3 or 4 knockout alleles for Meis1/Meis2 exhibited a body axis extension defect and were either similar in size to Uncx-stained E9.5 wild-type embryos (n = 3) or smaller (n = 4); comparison of somite size along the anteroposterior axis for 5 of these E10.5 mutants shows that somite sizes range from that seen in E9.5 wild type to about 60% of normal (Fig 6B-6D). We observed that E10.5 Meis1/Meis2 mutants carrying 3-4 knockout alleles that grew similar in size and somite number to E9.5 embryos exhibited a lack of forelimb bud outgrowth; n = 3 ( Fig 6E).
Overall, our findings show that loss of 3 or 4 alleles of Meis1 and Meis2 hinders body axis and forelimb formation, thus providing further evidence that our method of identifying candidate RA target genes can identify genes essential for development. Meis1/Meis2 double mutants all exhibit a growth defect but display a variable range of somite sizes (from normal to 60% of normal), which is not as severe as for Aldh1a2-/-embryos (consistent 50% reduction) [65]; however, Meis1/Meis2 double mutants that develop to E9.5 (20-25 somites) fail to develop forelimbs (Fig 6E), similar to E9.5 Aldh1a2-/-embryos [68]. In the future, more detailed studies of Meis1/Meis2 double mutants can be performed to determine how these genes control body axis and limb formation; plus, additional studies can be performed to determine how RARE enhancers function along with other factors to control Meis1 and Meis2 expression during early development.

Discussion
Our epigenetic ChIP-seq studies combined with RNA-seq on wild-type versus Aldh1a2-/-RAdeficient trunk tissue provides a means for identifying new candidate RA target genes that may be required for development. By focusing on RA-regulated genes that also have changes in nearby RA-regulated H3K27ac and/or H3K27me3 epigenetic marks associated with highly conserved RARE enhancers or silencers, our approach can be used to identify excellent candidates for gene knockout studies to learn more about gene function.
Here, in our studies on Aldh1a2-/-trunk tissue, we were able to narrow down 4,298 genes identified with RNA-seq that have significant changes in gene expression following loss of RA to 38 excellent candidate RA target genes in E8.5 trunk that also have significant changes in H3K27ac and/or H3K27me3 marks (located nearby or further away in the same TAD) associated with highly conserved RAREs. Our method allows one to identify genes that are most likely to be transcriptional targets of the RA signaling pathway as opposed to those whose expression is changed by effects downstream of RARs and RA signaling such as changes in expression or activity of other transcription factors or posttranscriptional changes in mRNA abundance. Our findings allow us to predict that some genes are likely to be indirect transcriptional targets of RA, as they have nearby RA-regulated peaks for H3K27ac or H3K27me3 but no RAREs, i.e., Pax6 that is transcriptionally regulated by factors whose expression is altered by loss of RA including Sox2 [45], Cdx [46], and Fgf8 [43].
Our findings provide evidence for additional RARE silencers. Previous methods designed to identify RAREs favored discovery of RARE enhancers, as studies were designed to find DNA elements that, when fused to a heterologous promoter and marker gene, would stimulate expression of the marker gene in the presence of RA. Also, when NCOAs and NCORs that control RA signaling were originally discovered, the model proposed for their function suggested that binding of RA to RAR favored binding of NCOA to activate transcription, with unliganded RAR favoring release of NCOA and binding of NCOR to repress transcription [69]. However, analysis of the Fgf8 RARE silencer at −4.1 kb demonstrated that RARs bound to RAREs can recruit NCOR in an RA-dependent manner; plus, this RARE is required for normal body axis extension [7]. The Fgf8 RARE silencer was also found to recruit PRC2 and histone deacetylase 1 (HDAC1) in an RA-dependent manner, providing further evidence that RA can directly control gene silencing [17]. Here, we identified additional RARE silencers near Fgf8 and Cdx2 plus several additional genes. Our studies indicate that RARE silencers are less common than RARE enhancers, and we found that Fgf8 is the only gene associated with a RARE silencer conserved beyond mammals. These additional RARE silencers can be further examined in comparison to the Fgf8 RARE silencer to determine the mechanism through which RA directly represses transcription. It will be important to determine how RAREs can function as RA-dependent enhancers for some genes but RA-dependent silencers for other genes.
RA has been shown to be required for balanced NMP differentiation during body axis formation by favoring a neural fate over a mesodermal fate [19,29,32]. Our studies provide evidence that RA directly regulates several genes at the trunk/caudal border needed for NMP differentiation-i.e., activation of Sox2 in the neural plate that favors neural differentiation, repression of Fgf8 that favors mesodermal differentiation, and repression of Cdx2 that helps define the location of NMPs. We now provide evidence for a candidate RARE enhancer that activates Sox2, 3 candidate RARE silencers that repress Cdx2, and 2 additional candidate RARE silencers for Fgf8. As the knockout of the original Fgf8 RARE silencer at −4.1 kb exhibited a body axis phenotype less severe than loss of RA in Aldh1a2-/-embryos [7], it is possible that the additional 2 candidate RARE silencers found here provide redundant functions for Fgf8 repression.
Our observation of highly conserved candidate RARE enhancers near 2 members of 2 different gene families (Nr2f and Meis) was intriguing, as it suggested that these gene family members may play redundant roles in body axis formation downstream of RA. As we were not sure whether the lack of previous studies on Nf2f1;Nf2f2 double-knockout and Meis1;Meis2 double-knockout mouse lines may be due to defects in double heterozygote adults that prevent generation of double homozygote embryos by conventional genetic approaches, we employed CRISPR gene editing to directly generate F0 double knockouts. Our Nr2f1/Nr2f2 doubleknockout studies indeed revealed a defect in body axis formation and small somites that is not observed in each single knockout. Interestingly, zebrafish nr2f1a/nr2f2 double-knockout embryos reported recently exhibit a heart defect more severe than each single knockout but not a body axis defect or body growth defect [52]. This observation is consistent with studies showing that RA is not required for NMP differentiation or body axis formation in zebrafish [70,71]. Thus, it appears that the ancestral function of Nr2f genes in fish was to control heart formation but that during evolution, another function to control body axis formation was added. Future studies can be directed at understanding the mechanism through which Nr2f1 and Nr2f2 control body axis formation.
The Meis1/Meis2 double knockouts we describe here revealed an unexpected function for Meis genes in body axis extension and forelimb initiation. Meis1 and Meis2 are markers of the proximal limb during forelimb and hindlimb development and were proposed to be activated by RA in the proximal limb as part of the proximodistal limb patterning mechanism in chick embryos [54,72,73]. However, knockout of Rdh10 required to generate RA demonstrated that complete loss of RA in the limb fields prior to and during limb development did not affect hindlimb initiation or patterning, whereas forelimbs were stunted but with Meis1 and Meis2 expression still maintained in a proximal position in both stunted forelimbs and hindlimbs [74,75] (reviewed in [5]). Our epigenetic results here support the previous proposal that RA can up-regulate Meis1 and Meis2 (but in the body axis prior to limb formation as opposed to the limb itself), and we provide evidence that Meis1 and Meis2 are transcriptional targets of RA in the body axis. Future studies can be directed at understanding the mechanism through which Meis1 and Meis2 control body axis and limb formation.
Our studies demonstrate the power of combining gene knockouts, ChIP-seq on epigenetic marks, and RNA-seq to identify transcription factor target genes required for a particular developmental process. In addition to H3K27ac and H3K27me3 epigenetic marks that are quite commonly observed near genes during activation or repression, respectively, it is likely that further ChIP-seq studies that identify RA-regulated binding sites for coactivators and corepressors will provide additional insight into RA target genes and transcriptional pathways. Such knowledge is essential for determining the mechanisms through which RA controls developmental pathways and should be useful to address RA function in adult organs. A similar epigenetic approach can be used to determine the target genes for any transcriptional regulator for which a knockout is available, thus accelerating the ability to understand gene regulatory networks in general.

Ethics statement
All mouse studies conformed to the regulatory standards adopted by the Institutional Animal Care and Use Committee at the SBP Medical Discovery Institute, which approved this study under Animal Welfare Assurance Number A3053-01 (approval #18-092). Animal care and use protocols adhered to the guidelines established by the National Institutes of Health (USA).

Generation of Aldh1a2-/-mouse embryos and isolation of trunk tissue
Aldh1a2-/-mice have been previously described [3]. E8.5 Aldh1a2-/-embryos were generated via timed matings of heterozygous parents; genotyping was performed by PCR analysis of yolk sac DNA. E8.5 trunk tissue was released from the rest of the embryo by dissecting across the posterior hindbrain (to remove the head, anterior hindbrain, pharyngeal region, and heart) and just posterior to the most recently formed somite (to remove the caudal progenitor zone) as previously described [17].

RNA-seq analysis
Total RNA was extracted from E8.5 trunk tissue (2 wild-type trunks and 2 Aldh1a2-/-trunks), and DNA sequencing libraries were prepared using the SMARTer Stranded Total RNA-Seq Kit v2 Pico Input Mammalian (Takara). Sequencing was performed on Illumina NextSeq 500, generating 40 million reads per sample with single read lengths of 75 bp. Sequences were aligned to the mouse mm10 reference genome using TopHat splice-aware aligner; transcript abundance was calculated using expectation-maximization approach; FPKM was used for sample normalization; generalized linear model likelihood ratio test in edgeR software was used as a differential test. High-throughput DNA sequencing was performed in the Sanford Burnham Prebys Genomics Core.

qRT-PCR analysis
Total RNA was extracted from 20 trunks of either E8.5 wild-type or Aldh1a2-/-embryos with the RNeasy Micro Kit (Qiagen #74004). Reverse transcription was performed with the High-Capacity cDNA RT Kit (Thermo Fisher Scientific #4368814). Quantitative PCR (qPCR) was performed using Power SYBR Green PCR Master Mix (Life Tech Supply #4367659). Relative quantitation was performed using the ddCt method with the control being expression of Rpl10a. Primers used for PCR (5 0 -3 0 ):

ChIP sample preparation for ChIP-seq
For ChIP-seq, we used trunk tissue from E8.5 wild-type or Aldh1a2-/-embryos dissected in modified PBS, i.e., phosphate-buffered saline containing 1X complete protease inhibitors (concentration recommended by use of soluble EDTA-free tablets sold by Roche #11873580001) and 10 mM sodium butyrate as a histone deacetylase inhibitor (Sigma # B5887). Samples were processed similar to previous methods [76]. Dissected trunks were briefly centrifuged in 1.5-ml tubes, and excess PBS dissection buffer was removed. For cross-linking of chromatin DNA and proteins, 500 μl 1% formaldehyde was added, the trunk samples were minced by pipetting up and down with a 200-μl pipette tip and then incubated at room temperature for 15 minutes. To stop the cross-linking reaction, 55 μl of 1.25 M glycine was added, and samples were rocked at room temperature for 5 minutes. Samples were centrifuged at 5,000 rpm for 5 minutes, and the supernatant was carefully removed and discarded. A wash was performed in which 1,000 μl of ice-cold modified PBS was added and mixed by vortex, followed by centrifugation at 5,000 rpm for 5 minutes and careful removal of supernatant that was discarded. This wash was repeated. Cross-linked trunk samples were stored at −80˚C until enough were collected to proceed, i.e., 100 wild-type trunks and 100 Aldh1a2-/-trunks to perform ChIP-seq with 2 antibodies in duplicate.
Chromatin was fragmented by sonication. Cross-linked trunk samples were pooled, briefly centrifuged, and excess PBS removed. A 490-μl lysis buffer (modified PBS containing 1% SDS, 10 mM EDTA, 50 mM Tris-HCl [pH 8.0]) was added, mixed by vortexing, then samples were incubated on ice for 10 minutes. Samples were divided into 4 sonication microtubes (Covaris AFA Fiber Pre-Slit Snap-Cap 6x16 mm, #520045) with 120 μl per tube. Sonication was performed with a Covaris Sonicator with the following settings: Duty, 5%; Cycle, 200; Intensity, 4; #Cycles, 10; 60 seconds each for a total sonication time of 14 minutes. The contents of the 4 tubes were recombined by transfer to a single 1.5-ml microtube, which was then centrifuged for 10 minutes at 13,000 rpm, and the supernatants transferred to a fresh 1.5-ml microtube. These conditions were found to shear trunk DNA to an average size of 300 bp using a 5-μl sample for Bioanalyzer analysis. At this point, 20 μl was removed for each sample (wild-type trunks and Aldh1a2-/-trunks) and stored at −20˚C to serve as input DNA for ChIP-seq.
Each sample was divided into four 100-μl aliquots to perform immunoprecipitation with 2 antibodies in duplicate. Immunoprecipitation was performed using the Pierce Magnetic ChIP Kit (Thermo Scientific, #26157) following the manufacturer's instructions and ChIP-grade antibodies for H3K27ac (Active Motif, Cat#39133) or H3K27me3 (Active Motif, Cat#39155). The immunoprecipitated samples and input samples were subjected to reversal of cross-linking by adding water to 500 μl and 20 μl 5 M NaCl, vortexing and incubation at 65˚C for 4 hours; then addition of 2.6 μl RNase (10 mg/ml), vortexing and incubation at 37˚C for 30 min; then addition of 10 μl 0.5 M EDTA, 20 μl 1 M Tris-HCl (pH 8.0), 2 μl proteinase K (10 mg/ml), vortexing and incubation at 45˚C for 1 hour. DNA was extracted using ChIP DNA Clean & Concentrator (Zymo, # D5201 & D5205). After elution from the column in 50 μl of elution buffer, the DNA concentration was determine using 2-μl samples for Bioanalyzer analysis. The 2 input samples ranged from 16-20 ng/l and the 8 immunoprecipitated samples ranged from 0.1 to 0.2 ng/μl (5-10 ng per 100 trunks). For ChIP-seq, 2 ng was used per sample to prepare libraries for DNA sequencing.

ChIP-seq genomic sequencing and bioinformatic analysis
Libraries for DNA sequencing were prepared according to the instructions accompanying the NEBNext DNA Ultra II kit (catalog # E7645S; New England Biolabs). Libraries were sequenced on the NextSeq 500 following the manufacturer's protocols, generating 40 million reads per sample with single read lengths of 75 bp. Adapter remnants of sequencing reads were removed using cutadapt v1.18 [77]. ChIP-Seq sequencing reads were aligned using STAR aligner version 2.7 to Mouse genome version 38 [78]. Homer v4.10 [79] was used to call peaks from ChIP-seq samples by comparing the ChIP samples with matching input samples. Homer v4.10 was used to annotate peaks to mouse genes and quantify reads count to peaks. The raw reads count for different peaks were compared using DESeq2 [80]. p Values from DESeq2 were corrected using the Benjamini and Hochberg (BH) method for multiple testing errors [81]. Peaks with BH-corrected p-value < 0.05 (BHP < 0.05) were selected as significantly differentially marked peaks. Transcription factor binding sites motif enrichment analyses were performed using Homer v4.10 [79] to analyze the significant RA-regulated ChIP-seq peaks; DR1 RAREs were found by searching for TR4(NR),DR1; DR2 RAREs by Reverb(NR),DR2; and DR5 RAREs by RAR:RXR(NR),DR5. Evolutionary conservation of RAREs was performed via DNA sequence homology searches using the UCSC genome browser software. TAD analysis was performed using the 3D Genome Browser (http://promoter.bx.psu.edu/hi-c/view.php); we used the TAD database from Hi-C data for mouse ES cells reported for the mouse mm10 genome. IPA was used to identify pathways for our list of target genes; from IPA results, heatmaps were designed with Prism software, and associated networks were created using STRING software. High-throughput DNA sequencing was performed in the Sanford Burnham Prebys Genomics Core, and bioinformatics analysis was performed in the Sanford Burnham Prebys Bioinformatics Core.

Generation of mutant embryos by CRISPR/Cas9 mutagenesis
CRISPR/Cas9 gene editing was performed using methods similar to those previously described by others [82,83] and by our laboratory [7]. sgRNAs were generated that target exons to generate frameshift null mutations, with 2 sgRNAs used together for each gene. sgRNAs were designed with maximum specificity using the tool at crispr.mit.edu to ensure that each sgRNA had no more than 17 out of 20 matches with any other site in the mouse genome and that those sites are not located within exons of other genes. DNA templates for sgRNAs were generated by PCR amplification (Phusion DNA Polymerase; New England Biolabs) of ssDNA oligonucleotides (purchased from Integrated DNA Technologies) containing on the 5 0 end a minimal T7 promoter, then a 20-nucleotide sgRNA target sequence (underlined below), and finally the tracrRNA sequence utilized by Cas9 on the 3 0 end, shown as follows: The 20-nucleotide target sequences used were as follows: Meis2 exon 2 (#2) CGACGCCTTGAAAAGAGACA sgRNAs were then transcribed from templates using HiScribe T7 High Yield RNA Synthesis Kit (New England Biolabs) and purified using Megaclear Kit (Life Technologies). sgRNAs were tested in vitro for their cleavage ability in combination with Cas9 nuclease (New England Biolabs); briefly, genomic regions flanking the target sites were PCR amplified, then 100 ng was incubated with 30 nM Cas9 nuclease and 30 ng sgRNA in 30 μl for 1 hour at 37˚C, followed by analysis for cleavage by gel electrophoresis.
For injection into mouse embryos, a solution containing 50 ng/μl Cas9 mRNA (Life Technologies) and 20 ng/μl for each sgRNA used was prepared in nuclease free water. Fertilized oocytes were collected from 3-4-week-old superovulated C57Bl6 females prepared by injecting 5 IU each of pregnant mare serum gonadotrophin (PMSG) (Sigma Aldrich) and human chorionic gonadotropin (hCG) (Sigma Aldrich). Fertilized oocytes were then transferred into M2 medium (Millipore) and injected with the Cas9 mRNA/sgRNA solution into the cytoplasm. Injected embryos were cultured in KSOMaa medium (Zenith) in a humidified atmosphere with 5% CO 2 at 37˚C overnight to maximize the time for CRISPR/Cas9 gene editing to occur at the 1-cell stage, then reimplanted at the 2-cell stage into recipient pseudopregnant ICR female mice. Implanted females were sacrificed to obtain F0 E9.0 embryos (Nr2f1/Nr2f2) or F0 E10.5 embryos (Meis1/Meis2). As fertilized mouse oocytes spend a long time at the 1-cell and 2-cell stages, this facilitates CRISPR/Cas9 gene editing at early stages and allows many F0 embryos to be examined for mutant phenotypes [7]. For genotyping, yolk sac DNA was collected, and PCR products were generated using primers flanking the sgRNA target sites; PCR products were subjected to DNA sequence analysis from both directions using either upstream or downstream primers. For each gene analyzed, embryos were classified as heterozygous (het) if the DNA sequence contained both a wild-type allele and a frameshift allele; embryos were classified as homozygous (hom) if only frameshift alleles were detected but no wild-type sequence.

Body axis length analysis of embryos
ImageJ software (https://imagej.net) [84] was used to measure body axis length along a 3-somite region (Nr2f1/Nr2f2 double mutants) or 4-somite region (Meis1/Meis2 double mutants) compared with wild type, with each specimen photographed at the same magnification. Statistical analysis was performed using 1-way ANOVA (nonparametric test) with data presented as mean ± standard deviation (SD) and with p > 0.05 indicating significance.

In situ gene expression analysis
Embryos were fixed in paraformaldehyde at 4˚C overnight, dehydrated into methanol, and stored at −20˚C. Detection of mRNA was performed by whole-mount in situ hybridization as previously described [85].