Abstract The integument plays a critical role in functional adaptation, with macro-regional specification forming structures like beaks, combs, feathers, and scales, while micro-regional specification modifies skin appendage shapes. However, the molecular mechanisms remain largely unknown. Craniofacial integument displays dramatic diversity, exemplified by the Polish chicken (PC) with a homeotic transformation of comb-to-crest feathers, caused by a 195–base pair (bp) duplication in HoxC10 intron. Micro-C analyses show that HoxC-containing topologically associating domain (TAD) is normally closed in the scalp but open in the dorsal and tail regions, allowing multiple long-distance contacts. In the PC scalp, the TAD is open, resulting in high HoxC expression. CRISPR-Cas9 deletion of the 195-bp duplication reduces crest feather formation, and HoxC misexpression alters feather shapes. The 195-bp sequence is found only in Archelosauria (crocodilians and birds) and not in mammals. These findings suggest that higher-order regulation of the HoxC cluster modulates gene expression, driving the evolution of adaptive integumentary appendages in birds. __________________________________________________________________ Higher-order HoxC cluster regulation controls region-specific epigenetic profiles and generates diverse skin appendages. INTRODUCTION Skin serves as the interface between an organism and its environment, providing a layer of protection. However, skin and its appendages vary in different body regions (e.g., scalp, face, trunk, paws, etc.) enabling regionally diverse functions to optimize adaptation to new ecospaces. The avian integument shows remarkable multilevel appendage diversity. Macro-level regional specificity refers to different appendage types, such as scales, feathers, beaks, combs, etc., and signaling pathways regulating scale versus feather phenotypes have been studied ([48]1, [49]2). In early ornithischian dinosaurs, filamentous appendages coexist with scales, preceding the origin of feathers ([50]3). Micro-level regional specificity refers to distinct shapes within a particular type of integumentary appendage such as beaks or feathers ([51]4). For feathers, birds exhibit radially symmetric down feathers, bilaterally symmetric contour feathers, and bilaterally asymmetric flight feathers in the trunk, tail, and wing, respectively, to optimize their functions in thermoregulation, communication display, and flight ([52]5, [53]6). These diverse appendages emerged gradually over the course of evolution, eventually converging into the region-specific integumentary arrangement seen in modern bird ([54]3, [55]7). For example, early theropods started to display diverse integumentary appendages as seen in the Beipiaosaurus, which displays regionalized feathers ([56]8). The theropod, Sinosauropteryx, exhibits downy-like protofeathers all over the body ([57]9). In Microraptor, all four limbs were shown to have asymmetric feathers. Through evolution, skin appendages in the hindlimbs became scales. Apparent regional differences between flight and tail feathers in Caudipteryx imply a novel molecular mechanism regulating the evolution of skin regional specificity ([58]10). In Longirostravis, unique crown feathers were observed on the cranium, likely serving a display (fig. S1A) ([59]11). Furthermore, with stem cells that allow feather cycling, old feathers are replaced with new ones in each cycle ([60]12–[61]14). Feather transitions from natal down to juvenile feathers ([62]15), offered birds the opportunity to molt and alter feather phenotypes in response to hormonal and seasonal changes ([63]16, [64]17), making the integument a highly interactive biointerface that could adapt to varied environments. Thus, regionalizing the avian integument for flight or display was and continues to be a major driving force for successful avian evolution. In Discussion, we will further discuss integumentary fossil records demonstrating regional specification in feathered dinosaurs and Mesozoic birds. How does skin regional specification form? This can be considered from the tissue interaction level, morphogen level, and epigenetic level. At the tissue interaction level, classical epithelial-mesenchymal recombination experiments demonstrated that the dermis controls the skin appendage phenotype ([65]18, [66]19). At the morphogen level, Wnt-β-catenin, fibroblast growth factor (FGF), and bone morphogenetic protein (BMP) pathways have been shown to determine appendage distribution and phenotypes ([67]20–[68]29). Diffusible morphogens controlling different aspects of feather morphology have been identified: BMP/noggin regulate barb branch formation ([69]30), Wnt3a positions the rachis ([70]24), Wnt2b controls the formation of vanes versus fluffy barbs ([71]6), and retinoic acid mediates mediolateral asymmetry formation ([72]5). In mice, Wnt is involved in hair growth ([73]23). Dkk2, a Wnt inhibitor, is expressed in hairless plantar skin ([74]25). Normal laboratory mice show short hair with low Wnt and Hox expression in their ear skin. The Koala mouse mutant exhibits long hairs in the ear with high Wnt activity and HoxC expression levels in their ear skin ([75]31). While we have established a general logic on how diffusible morphogens (e.g., Wnt, FGF, BMP, etc.) work together to form specific feather phenotypes, we need to learn how different skin regions express different combinations of these morphogens at the epigenetic level. We speculate that coupling region-specific epigenetic profiles with different morphogenetic modules helps establish region-specific skin phenotypes. A comparative transcriptome study revealed multiple regulatory modules involved in the scale-feather transition ([76]2). Hox genes have been considered to provide positional information during vertebrate development ([77]32–[78]35). For the skin, Hox codes were postulated to define chicken appendages in different body regions ([79]36, [80]37) and also were shown to be differentially expressed in fibroblasts derived from different human skin regions ([81]38). Hox gene expression was also linked to region-specific differences in skin coloration in crows ([82]39) and finches ([83]40). Mutations in coding or noncoding regions can cause homeotic transformations in which an appendage forms within a different region of the body ([84]41). In chickens, recent studies have pinpointed that genomic structural variations are related to the integumental variation in chicken breeds ([85]42–[86]45). Obvious scalp skin appendage differences exist between White leghorn (WL) chickens and Polish chickens (PC). The scalp of WL and most other chicken breeds has a large comb adjacent to short cranial feathers ([87]Fig. 1A), whereas the PC scalp bears elongated crest feathers surrounded by short cranial (non-crest) feathers ([88]Fig. 1B). WL scalp skin consists of three domains: the comb domain, the comb-base domain (an apteric region separating the comb and feather domains), and the feather domain, whereas PC scalp was occupied by crest and non-crest feathers ([89]Fig. 1C). The PC crest locus was mapped to the HoxC cluster ([90]46). Later, a 195–base pair (bp) duplication in the intron of HoxC10 was identified as the causal mutation for the crest phenotype in PC and other crested chickens ([91]47). This prompted us to study in depth the regulation of HoxC genes in different avian skin regions. Fig. 1. Colinear HoxC expression patterns and chromatin accessibility in chicken skin, along the AP axis. [92]Fig. 1. [93]Open in a new tab (A) The WL chicken displays different skin appendage types (i.e., comb and feather) and morphologies (feather shapes) along the body AP axis. (B) The PC head grows back-like contour feathers in the crest region. Note the morphology and size differences between the crest and non-crest feathers in the PC head. (C) Schematic drawing shows the top view of the developing WL and PC head. Note different comb, comb-base, and feather domains in WL and the crest or non-crest feathers in PC. (D) Transcriptome expression heatmap of HoxC gene cluster expression in E9 WL. Five regions from the scalp to the tail are used. Colinear expression is shown in the HoxC gene cluster. (E) SISH of HoxC8, HoxC10, and HoxC12 in E7 chicken samples show no expression in the scalp and the strongest signal in sequential anterior domains. Scale bars, 200 μm. (F) Summary of the HoxC gene expression pattern in embryonic WL chicken skin. (G) Transcriptome expression heatmap of HoxC gene cluster expression in the adult WL feather follicles along the body AP axis. (H and I) Gene expression heatmaps of significantly up- (H; FC >10) and down-regulated (I; FC ≤10) scalp/body/tail-specific marker genes. (J) Integrated Genomics Viewer showing the E9 scalp, back, and tail ATAC-seq and RNA-seq data within the HoxC cluster. Region-specific DAR peaks that on average obtain greater than twofold accessibility compared to the other two regions are indicated with arrows. Purple arrows point to promoters, while red arrows indicate other noncoding regulatory regions. (K) Illustration showing our model depicting how different open chromatin profiles direct regional specificity. ante. back, anterior back; DEGs, differentially expressed genes; FC, fold change; post. back, posterior back. Large-scale regulation of gene clusters often involves higher-order chromatin organization ([94]48, [95]49). The elementary units of genome architecture are topologically associating domains (TADs) ([96]50). A TAD contains sites of internal chromosomal interactions and limits interactions with sequences residing outside the TAD ([97]51). We hypothesize that regional differences in TADs can create region-specific gene expression altering tissue morphologies and their functions. Here, we study the three-dimensional (3D) configurations of chromatin containing HoxC clusters in different body regions using Micro-C. We find that the HoxC cluster TAD is closed in the WL scalp but is open in the PC scalp. We integrate RNA sequencing (RNA-seq) and assay for transposase-accessible chromatin using sequencing (ATAC-seq) with Micro-C data to identify connections with external sites that may regulate HoxC expression in different body regions. We further demonstrate that in vivo CRISPR-mediated deletion of the 195-bp regulatory sequence and the misexpression of HoxC members can perturb skin region specificity. Thus, this study presents an example of a platform for understanding how different regulatory landscapes, with spatiotemporally controlled gene expression profiles, can lead to region-specific skin appendage phenotypes. RESULTS The distribution of HoxC gene expression along the AP axis of WL skin We examined the HoxC gene expression profile in chicken skin regions along the anterior-posterior body axis (AP axis) of WL by performing bulk RNA-seq analysis from E9 dermis in five regions along the AP axis (scalp, neck, anterior back, posterior back, and tail). An RNA-seq heatmap demonstrates spatial colinear HoxC cluster gene expression ([98]Fig. 1D). Few HoxC genes are expressed in the scalp. HoxC4 is moderately expressed starting from the neck, whereas HoxC5-C8, HoxC9-C10, and HoxC11-C13 display high expression starting from the anterior back, posterior back, and tail, respectively. Whole-mount in situ hybridization (WISH) confirmed the colinear HoxC distribution (fig. S2, A and A′). Section in situ hybridization (SISH) shows that the spatial distribution of HoxC8, HoxC10, and HoxC12 begin in the anterior back, posterior back, and tail skin, respectively, but are absent in the scalp skin ([99]Fig. 1E, and see different regions in fig. S1B). Hence, along the AP axis, HoxC gene expression appears to be colinear with their genomic locations ([100]Fig. 1F). We next examined whether adult WL feather follicles retain the embryonic HoxC expression pattern. In adult WL, feathers grown from different regions are morphologically distinct ([101]Fig. 1A). We performed RNA-seq on adult feather follicles from the same regions used for embryonic skin (fig. S1, C to E). The HoxC RNA expression adult profile parallels those seen in embryonic skin and demonstrates colinear expression patterns ([102]Fig. 1G). To find the molecular signatures of feathers from different adult body regions, we focused on the scalp, anterior back, and tail. We identified differentially expressed genes (DEGs) among these regions [defined as fold change (FC) >2 or <−2 with a false discovery rate (FDR) < 0.05] (fig. S1F). More stringent criteria (FC >10 or <−10; FDR < 0.05) identified 54 up-regulated and 65 down-regulated region-specific transcripts ([103]Fig. 1, H and I). Hox genes on this list indicate that they may be involved in establishing regional specificity. The colinear HoxC expression corresponds to chromatin accessibility How chromatin accessibility differs in various skin regions resulting in differential HoxC cluster expression profiles was explored by conducting ATAC-seq on E9 scalp, anterior back, and tail dermis (fig. S3, A and B). ATAC-seq peaks (542,955) were identified. To understand region-specific chromatin architecture, we compared genome-wide differential accessible regions (DARs) among these samples (FC >2 or <−2; P value < 0.05) and found specific ATAC-seq peaks for each region (fig. S3, C and D). The region-specific open chromatin segments ([104]Fig. 1J) are classified into seven gene regions based on their genomic locations [promoter, intergenic, 5′ untranslated region (5′UTR), 3′UTR, intron, transcription termination site (TTS), and coding sequence (CDS); fig. S3E]. The scalp- and tail-specific peaks are enriched in intronic and intergenic sections, suggestive of enhancers, while the back-specific peaks are enriched in promoter regions. To unveil DAR functions, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The top 10 GO terms of biological processes and KEGG pathways (fig. S3F) show that Wnt signaling is one of the enriched morphogen pathways. We found that scalp-specific ATAC peaks were absent in the HoxC cluster. Back-specific peaks were present toward the 3′, and tail-specific peaks were toward the 5′ region of the HoxC cluster ([105]Fig. 1J). These peaks belong to either promoter (purple arrows) or other noncoding regulatory regions (red arrows). One back-specific peak is in the HoxC5 promoter region, and seven peaks are in the noncoding regulatory regions from HoxC4-C8. Tail-specific peaks are found in the HoxC12 and HoxC13 promoters as well as several peaks in the noncoding regulatory regions from HoxC6-C13. Integrating the RNA-seq data with chromatin status shows that scalp dermis without open chromatin segments lacks HoxC gene expression. In contrast, back and tail samples with open chromatin regions display HoxC6-C10 and HoxC10-C13 expression, respectively. Therefore, the DARs generally follow spatial Hox colinearity, suggesting that chromatin accessibility could contribute to region-specific HoxC gene expression ([106]Fig. 1K). The higher-order chromatin organization of the HoxC gene cluster is skin region specific in the WL To understand how chromatin architecture relates to region-specific HoxC gene expression, we performed Micro-C experiments using E9 [Hamburger and Hamilton stage 35 (HH35)] WL scalp, back, and tail skin dermis. Micro-C is a genome-wide chromosome conformation capture 3C-based method that provides chromatin interactome data at single nucleosome resolution ([107]52). For comparison purposes, we also included RNA-seq and ATAC-seq data for these samples ([108]Fig. 2, A to C). The HoxC gene cluster is located on chromosome 33 according to the chicken galGal6 assembly. By identifying significant TADs and chromatin contact regions (CRs, or anchors that form the loops), we observed a scalp-specific TAD that covers the whole HoxC cluster ([109]Fig. 2A, green triangle and horizontal black bar). In contrast, this TAD is absent in both back and tail dermis ([110]Fig. 2, B and C, dotted green triangles). Moreover, in scalp dermis, most chromatin contacts were observed within the scalp-specific TAD ([111]Fig. 2A, red arches), whereas in back and tail dermis, numerous loops were detected outside of the HoxC cluster ([112]Fig. 2, B and C, gray arches). These findings led us to hypothesize that the TAD in the WL scalp may set up barriers around the HoxC cluster and inhibit gene expression by restricting promoter and distal enhancer interactions. Fig. 2. Different skin regions display distinct higher-order chromatin configurations in the WL chicken HoxC gene cluster. [113]Fig. 2. [114]Open in a new tab (A to C) The Micro-C contact maps, ATAC-seq, and gene expression profiles of WL E9 scalp (A), back (B), and tail (C) skin dermis at the HoxC gene cluster (chr33:7,435,000–7,715,000) demonstrate a significant scalp-specific TAD [chr 33: 7,490,000–7,670,000, (A), black bar] and intracluster (red arching lines), and off-cluster (gray arching lines) chromatin contacts around the cluster. Purple-shaded boxes indicate the whole HoxC cluster. Vertical pink line, region-specific DARs at intracluster loop anchors; vertical green line, ATAC-seq peak regions at intracluster loop anchors; vertical blue line, intracluster loop anchors contain no ATAC-seq peaks and DARs; black arrows and red asterisks indicate the 195-bp region. (A′ to C′) Schematic drawing of accessible chromatin, intracluster, and adjacent chromatin contacts surrounding the HoxC gene cluster. (D) A close look of the dashed box in (B) with additional histone H3K27ac, H3K4me1, H3K4me3, H3K27me3 chromatin immunoprecipitation sequencing, and nucleosome profiles. The green horizontal bar indicates a clustered enhancer (Super-E) analyzed by the rank order of super-enhancer (ROSE) software. Vertical bars indicate a E9 back dermis-selective DAR (red) and accessible chromatin (green), and chromatin contacts without ATAC-seq peaks (blue). Red asterisks indicate the 195-bp region. Brown bars, chromatin CRs or so-called loop anchors. (E) Schematic drawing summarizes off-cluster chromatin contacts to the HoxC cluster in E9 WL chicken skin dermis. To better understand how chromatin at the HoxC locus is organized, we further categorized chromatin loops into two types—(I) intracluster contacts ([115]Fig. 2, A to C, red arches, and [116]Fig. 2, A′ to C′, light brown squares) and (II) off-cluster contacts ([117]Fig. 2, A to C, gray arches, and [118]Fig. 2, B′ and C′, dark brown squares). In type I, both loop anchors are within the HoxC clusters, whereas in type II, only one loop anchor is within the HoxC cluster. On the basis of ATAC-seq, we found that most of the type I contacts contain either region-specific DARs or accessible chromatin. Moreover, we observed that there are four conserved contacts among different skin regions. These contracts are located in a 5-kb region that lies upstream of the HoxC13, HoxC10, HoxC6, and HoxC4 genes. Different skin regions may use these conserved contacts, along with additional region-specific contacts, to create unique chromatin conformations. The WL has only one copy of the 195-bp sequence in the HoxC10 intron compared to the two copies that are present in the PC ([119]47). We wonder whether this single 195-bp copy in WL locates within any loop anchor among different skin regions. We found that there is one loop anchor containing this 195-bp sequence in the WL back and tail samples ([120]Fig. 2, B and C, red asterisks). The duplicated 195-bp in PC correlated to high HoxC10 expression in its scalp ([121]47), implying a possibility that the duplicated 195-bp segment may alter chromatin topology of the HoxC cluster in PC. To gain insight into potential intracluster regulation of the HoxC chromatin structure, we mapped histone modification profiles using histone chromatin immunoprecipitation sequencing (H3K27ac, H3K4me1, H3K4me3, and H3K27me3) ([122]53), focusing on the HoxC cluster in a back sample. As shown in [123]Fig. 2D, we found a super-enhancer (Super-E, 80 kbp, horizontal green bar) covers the region from HoxC11 to HoxC5. In addition, 2 intracluster and 10 off-cluster contacts ([124]Fig. 2B) located within the super-E were identified. Since the HoxC cluster has more off-cluster contacts when the TAD is open in WL back and tail regions ([125]Fig. 2, B and C), we further examined the 5-kb CRs and the genes embedded therein (fig. S4, A and B, and table S1). We identified 12 and 6 significant off-cluster CRs in the WL back and tail samples, respectively, but no CR was found in the scalp ([126]Fig. 2E). We categorized these CRs as proximal contacts (CR10-12 and CR18), which flank the HoxC cluster, and distal contacts, which are at least 400 kb away from the HoxC cluster. In table S1, we further identified the TWIST3, NACA, TUBA1A, PRKAG1, SP1, and HNRNPA1 genes located in the back sample CRs. Another group of genes, NCKAP1L, TUBA1A, PRKAG1, and HNRNP1 genes were located in the tail sample CRs. All these genes had medium (11 to 1000 transcripts per million, TPM) to high (TPM > 1000) expression levels based on the RNA-seq results. These results suggest that the open TAD HoxC gene loci might serve as a regional morphogenetic hub, which contacts multiple morphogenetic modules to establish skin region–specific expression profiles. Overall, the WL Micro-C and histone map suggest that the scalp has a TAD and the regulatory regions fall within the TAD. However, the WL back and tail do not have a TAD covering the HoxC cluster. The lack of TAD status in the back and tail regions may enable the accessibility of morphogenetic modules located outside of the HoxC locus. This may account for the lack of HoxC gene expression in the scalp and high HoxC gene expression in the back and tail. The 195-bp duplication in HoxC10 regulatory regions allows HoxC genes to make long-distance contacts outside of the TAD in the PC scalp Genetic data demonstrated that the PC crest phenotype is caused by a 195-bp duplication in the HoxC10 intron ([127]47). To identify upstream factors that regulate differential HoxC expression in the PC scalp, we performed ATAC-seq on the crest feathers (crest) and the neighboring non-crest feathers ([128]Fig. 1C). Peaks (5265 and 23,980) were more accessible in the crest and non-crest feathers, respectively (fig. S5, A and B). Compared with the non-crest samples, the crest samples are enriched for reads at transcription start sites (TSSs) ([129]Fig. 3A). In addition, the PC crest has more accessible promoter peaks ([130]Fig. 3B), which is similar to WL back (fig. S3E). These results suggest that the PC crest obtains higher promoter activity to produce long crest feathers. Fig. 3. The TAD is lost within the HoxC cluster in PC scalp. [131]Fig. 3. [132]Open in a new tab (A) TSS profile plot and heatmap showing the read coverage near the TSSs within the ATAC-seq peaks. TSS, transcription start site. Y axis is the total read counts. (B) Distribution of genomic features among the non-crest versus crest-specific ATAC-seq peaks. UTR, untranslated region; TTS, transcription termination site; CDS, coding sequence. (C) Profiles of ATAC-seq and RNA-seq show DARs (green vertical bars) within the HoxC gene cluster. A crest-specific ATAC-seq peak is shown (blue vertical bar). (D) H&E and SISH analyses showing the distributions of DLX3, LMX1B, and SREBF2 in crest and non-crest feathers. (E) The Micro-C contact maps of PC scalp demonstrate 3D chromatin structures and significant chromatin loops surrounding the HoxC gene cluster. Black arrowheads show a loop between the 195-bp duplicate (red asterisk) and a putative enhancer (green arrowhead). ATAC-seq tracks show two PC scalp–specific DARs (vertical red bars) and two ATAC-seq peaks (vertical green bars). (F) Hypothetical model illustrating that DLX3 and MSX1 target their binding motifs. DLX3 binds to a skin region–specific DAR, while MSX1 may interact with an ATF7 intron. This interaction hub lies beside the 195-bp duplicate and may mediate de novo loop formation. (G) Schematic drawing showing off-cluster chromatin contacts of the PC scalp HoxC cluster. Purple square, HoxC gene cluster. Red square, CRs contact the HoxC10 gene. Black square, CRs contact the HoxC gene cluster. (H) Hypothetical model illustrating how skin region–specific 3D chromatin architectures modulate HoxC gene expression. A single 195-bp sequence allows the formation of a TAD that covers the HoxC gene cluster, which restricts off-TAD chromatin contacts. The duplicated 195-bp sequence prevents the formation of the TAD and allows HoxC proximal and distal chromatin loops to form, resulting in Hox and downstream gene activation. We identified two DARs at the HoxC locus in the crest sample (fig. S5C). These DARs are in the HoxC10 and HoxC9 promoters ([133]Fig. 3C, green vertical bars), supporting a correlation between increased chromatin accessibility and up-regulated gene expression in the crest. We next examined the enriched transcription factor (TF)–binding motifs in crest-specific DARs. We found five significantly enriched binding motifs (Dlx3, Dlx6, MSX1, En1, and Lmx1b) in both DARs (fig. S5D). We also found an Srebf2 motif in DAR 2 (fig. S5D). SISH shows that these TFs are up-regulated in the crest samples ([134]Fig. 3D). These findings suggest that these TFs might contribute to HoxC gene cluster regulation. Furthermore, we conducted Micro-C experiment to analyze the 3D chromatin organization of PC crest feathers. Unlike the WL scalp sample, no TAD was detected in the PC crest HoxC cluster ([135]Fig. 3E, compared to [136]Fig. 2A). The lack of a TAD in PC scalp resembles the WL back and tail samples ([137]Fig. 2, B and C). Next, we identified significant PC crest HoxC locus chromatin loops ([138]Fig. 3E, gray arches). All these loops are off-cluster contacts. Our previous work suggests that the 195-bp duplication in the HoxC10 gene intron is associated with the PC crest phenotype ([139]47). Excitingly, we found that the duplicated 195-bp sequence coincided with one of the chromatin loop anchors (red asterisk). Furthermore, this anchor also covers DAR 1 at the HoxC10 promoter ([140]Fig. 3F). Another loop anchor is detected in an ATF7 gene intron that lies 130 kb upstream of the HoxC cluster ([141]Fig. 3F, named CR21 in fig. S5E). An open chromatin region at the ATF7 anchor could function as an enhancer and interact with the open HoxC10 promoter, resulting in transcriptional activation of both HoxC10 and ATF7 ([142]Fig. 3F). Thus, our data suggest that the 195-bp duplication in PC scalp disrupts the formation of the HoxC TAD, resulting in formation of multiple off-cluster chromatin contacts, which may be critical to establish the crest feather-forming chromatin interactome. To understand more about these off-cluster HoxC chromatin contacts in the PC scalp, we also examined the 5-kb CRs and genes embedded in them. We identified six significant CRs and named them CR19-CR24 ([143]Fig. 3G, fig. S5E, and table S1). We summarize the HoxC gene cluster chromatin contacts in fig. S5F. These results show that chromatin contacts are restricted within the TAD in the WL scalp sample. The absence of TADs in WL back/tail regions and the PC scalp enabled the formation of off-cluster proximal and distal chromatin contacts. Together, our findings suggest a hypothetical model ([144]Fig. 3H) illustrating how skin region–specific 3D chromatin architectures modulate expression of key HoxC genes and contribute to the skin appendage phenotype. A single 195-bp sequence within the HoxC10 gene intron produces a TAD that covers the whole HoxC gene cluster in the WL scalp, restricting off-TAD chromatin contacts. The duplicated 195-bp sequence blocks TAD formation, allowing HoxC genes to form proximal and distal chromatin loops. These loops may result in HoxC gene activation and/or activation of genes associated with their downstream CRs. CRISPR-Cas9–mediated deletion of the 195-bp region inhibits long feather growth in the PC scalp To further evaluate the function of the 195-bp repeat region, we used CRISPR-Cas9–directed excision to remove the duplicated copy from the PC genome. To achieve this, we need an effective tool to remove one copy of 195 bp. To increase CRISPR efficiency in chicken embryos, we modified a single CRISPR plasmid targeting chicken tissues ([145]54) (fig. S6A) to generate a Tol2-CRISPR plasmid also encoding Citrine, a fluorescent marker of cells containing the CRISPR plasmid. DNA transposons are powerful tools for analyzing the regulatory genome ([146]55). We first injected a plasmid (Tol2-CRISPR-CTNNB1) targeting chicken β-catenin (CTNNB1) to the HH16 somite. Embryos collected at HH35 show reduced CTNNB1 expression and severe skin abnormalities, including thinner skin and a lack of feather bud formation (fig. S6, B and C). In DF1 chicken embryo fibroblasts, Tol2-CRISPR-CTNNB1 produced 57.3% CRISPR-mediated deletion efficiency (fig. S6D). We then targeted the PC duplicate 195-bp sequence by creating a Tol2-CRISPR–195-bp plasmid ([147]Fig. 4A). Its efficiency was tested in cultured PC embryonic day 10 (E10) crest fibroblasts by flow cytometry and extracting DNA from Citrine-positive cells. Polymerase chain reactions (PCRs) covering the duplicated 195-bp region revealed that there is a shorter DNA band in the CRISPR–195 bp–treated Citrine-positive cells ([148]Fig. 4B, red arrow). We cloned this shorter band into the p-drive plasmid and sequenced 10 clones, which all showed only one 195-bp copy ([149]Fig. 4C), demonstrating that we successfully deleted the duplicated copy in some cultured PC fibroblast cells. Fig. 4. CRISPR-Cas9 targeting of the 195-bp duplication in PC inhibits crest feather growth. [150]Fig. 4. [151]Open in a new tab (A) CRISPR-cas9 design strategy. We aim to delete one copy of the duplicated 195-bp region in PC using CRISPR-Cas9–directed excision. (B) CRISPR-Cas9 experiment using cultured PC scalp fibroblast cells. Cells are sorted using green fluorescent protein after 2 days of culture. The control sample only shows a 579-bp band representing the duplicated 195-bp sequence. There is a 384-bp band in CRISPR-Cas9–treated Citrine-positive cells (red arrow). (C) Sequencing shows the CRISPR control and CRISPR-195–bp shorter PCR band. Note that the shorter band contains only a single 195-bp copy. (D) CRISPR-cas9 deletion of the duplicated 195-bp region in PC reduced the crest feather phenotype. Region 1 is close to the beak (red box). Region 2 is at the top of head (green box). Deletion of the duplicated 195-bp sequence in PC generated three phenotypes. Type I, ectopic apteric zone (no feather growth); type II, small feather buds; type III, short feather buds with random orientation. (E) RNAscope shows the reduction of HoxC6/8/10 expression in region 1 of CRISPR-treated samples. Results for individual channels of region 1 and region 2 are shown in fig. S6. (F) RNAscope of LGR4/5/6 gene expression in region 1 of CRISPR-treated samples. Results for individual channels of region 1 and region 2 are shown in fig. S6. (G) qPCR analysis shows the down-regulation of HoxC6, HoxC8, and HoxC10 in CRISPR-treated PC samples. The data are analyzed using GraphPad Prism software. The CRISPR–195-bp plasmid was delivered to the HH11 PC neural crest region. We targeted migrating neural crest cells since classical neural crest transplantation work showed that they give rise to the craniofacial dermis and determine head feather patterning ([152]56, [153]57). Samples were collected at HH39 (N = 14). These PC samples show three phenotypes ([154]Fig. 4D and fig. S6E). For analysis, we highlighted two regions of the head. Region 1 is close to the beak at the anterior boundary of crest feathers (red box). Region 2 is in the scalp (green box). All three phenotypes appeared in region 1 but not in region 2. The type I phenotype has an ectopic apteric zone. The type II phenotype shows inhibited feather bud growth. The type III phenotype displayed short, disoriented feather buds. CRISPR control samples (N = 6) did not produce any phenotypes. We next used RNAscope to localize changes in Hox gene expression. RNAscope analysis revealed reduced expression of HoxC6/8/10 in region 1 for each of the three phenotypes ([155]Fig. 4E and fig. S6F). The expression of LGR4/5/6 in region 1 was decreased in type I and II but not in type III phenotypes ([156]Fig. 4F and fig. S6G). The reduction of HoxC6/C8/10 expression then was confirmed by real-time quantitative PCR (qPCR) ([157]Fig. 4G). This result suggests that the CRISPR-mediated deletion of the duplicated 195-bp sequence altered the expression of HoxC genes and suppressed the PC crest phenotype. However, we realize that the in vivo CRISPR efficiency is not sufficient to generate more dramatic phenotypes, but the results support the predicted role of the 195-bp duplication in crest feather formation. A parallel CRISPR–195-bp experiment conducted in WL did not induce any phenotype changes (N = 6) (fig. S7). HoxC modulates the specification of feather and comb regions via MSX2 and Wnt, among other signaling molecules It has been shown that crested chickens express higher HoxC8 levels than non-crested chickens ([158]46, [159]47), suggesting that HoxC8 may promote crest phenotype formation. We examined the endogenous HoxC8 expression on the WL scalp during development from E7 to E13 ([160]Fig. 5A). HoxC8 transcripts were not detected before E10 when they are expressed in developing feather buds. HoxC8 is also expressed in the tip of the comb primordia at moderate levels but is negative in the comb-base domain. These results suggest that HoxC8 may play a role in specifying the domains of different scalp skin appendages starting from E10. Fig. 5. HoxC genes regulate skin appendage phenotype in the WL chicken scalp. [161]Fig. 5. [162]Open in a new tab (A) WISH shows HoxC8 expression patterns on the E7–E13 chicken scalp. HoxC8 expression can be detected starting at E10 in the feather and comb domains. (B) Schematic drawing of RCAS-mediated misexpression. (C) Ectopic HoxC8 expression expands the feather domain and decreases the size of the comb-base domain. Whereas, dnHoxC8 has the inverse effect. (D) Quantification of the effects of ectopic HoxC8 expression and dnHoxC8 on the size of the apteric comb-base domain (mm^3); N = 5 for each condition. **P < 0.01. (E) An example of ectopic expression of a gene up-regulated by HoxC8. RCAS-MSX2 decreases the comb-base domain. (F) An example of suppression of a loop anchor point. Expressing a dominant negative form of ATF7 increases the size of the apteric comb-base domain. (G) RNAscope analysis shows ectopic HoxC8 expression (second row) leads to reduced comb-base domain size accompanied by decreased SRFP4, FN, and SHH expression compared to controls (first row). In the dominant negative ATF7-treated sample, the comb-base domain is increased while SPRF4, FN, and SHH expression is reduced (third row). White arrow indicates the reduced comb expression. (H) Summary diagram showing that the WL comb base size is decreased by ectopic HoxC8 or MSX2 expression but increased by HoxC8 or ATF7 inhibition. (I) In normal WL scalp, the after-feather is about half the length of the main feather. Ectopic HoxC8 expression increased the after-feather length to be two times the length of the main feather. (J) Feathers in the WL scalp are exclusively pennaceous. Ectopic expression of pooled HoxC genes caused the plumulaceous feather region to occupy ~50% of the total length. c, comb domain; f, feather domain; cb, comb-base domain, Pen, Pennaceous; Plu, Plumulaceous. To test HoxC8 gene function, replication-competent avian sarcoma (RCAS) virus was used to modulate HoxC8 expression. The virus was injected into WL E4 scalp skin, and the embryos were harvested at E13 to characterize skin phenotype changes ([163]Fig. 5B). The ectopic HoxC8 expression was traced using a mCherry fluorescent protein, SISH, and AMV-3C2 immunostaining to detect RCAS virus infection (fig. S8, A to C). Our data show that, compared to controls, ectopic HoxC8 expression in the scalp shrinks the comb base but expands the feather domains. In contrast, when HoxC8 activity is suppressed with a dominant negative construct (dnHoxC8), the comb-base domain expanded dramatically ([164]Fig. 5, C to D, and fig. S8C). RCAS-dnHoxC10 also caused a similar phenotype that increased the comb-base domain (fig. S8D). Hence, modulating HoxC8/C10 expression affects scalp skin regional specificity. To identify downstream mediators and possible signaling pathways associated with the phenotypic change, we performed RNA-seq on HoxC8 overexpressing WL samples (fig. S8, E and F). The control was a sample with a virus expressing only an mCherry fluorescent protein. Differential expression analysis was performed to find significant DEGs (FC >2 or <−2 and a P value <0.05). Up-regulated genes (1303) and 954 down-regulated genes were identified (fig. S8G). Several Wnt signaling pathway members (Wnt4, Wnt6, Wnt10A, Wnt16, AXIN2, LGR5, etc.) were up-regulated (fig. S8H). This parallels the finding of high Wnt expression in Koala mouse mutant ears that grow long hair ([165]31). Intriguingly, two genes, MSX2 and MAT1A, containing the HoxC8-binding motif (fig. S8I) in their upstream segments, were also up-regulated (fig. S8H) and are likely direct downstream targets of HoxC8. In addition, we identified canonical pathways for the DEGs using Ingenuity Pathway Analysis software. Wnt signaling is the most significantly enriched among the common morphogenic pathways (fig. S8J). To examine molecules acting downstream of HoxC8 on scalp fate determination, we overexpressed MSX2 in WL. We found that RCAS-MSX2 decreased the apteric comb-base domain size but increased the feather domain size ([166]Fig. 5, E versus C, control), mimicking one aspect of the PC crest phenotype. Furthermore, we tested the function of ATF7, a loop anchor point that is upstream of HoxC10 and C8 ([167]Fig. 3F). We suppressed ATF expression with a dominant negative form (dnATF7) ([168]58). After injecting RCAS-dnATF7 in the WL scalp, we observe shortened feather growth and an enlarged apteric comb-base domain ([169]Fig. 5, F versus C, control) (N = 3/3). Further characterization by RNAscope analysis revealed that the smaller comb-base domain in RCAS-HoxC8–infected samples is accompanied by the reduced expression of SFRP4, FN and SHH in the comb ([170]Fig. 5G, arrows). These molecules were chosen because they play important roles in tissue morphogenesis ([171]59, [172]60). Their expression was also reduced in the comb and comb-base in RCAS-dnATF7–transfected samples ([173]Fig. 5G, third row). Together, these results suggest that variation of HoxC8 expression can modulate the chicken scalp domain position by modulating Wnt signaling. These functional studies suggest that the levels of HoxC8 expression and their downstream genes (i.e., MSX2) regulate the WL scalp domain size (i.e., comb, comb-base, or feather). Furthermore, ATF7 acts upstream to regulate HoxC8 expression and changes the scalp domain size ([174]Fig. 5H). ATF7, MSX1, and MSX2 expression shows dramatic differences between the PC and WL crest region To further understand the roles of upstream regulators, we characterized the expression of ATF7, MSX1, and MSX2 in E9 WL and PC head by WISH. Faint expression of ATF7 and MSX1 and strong expression of MSX2 were detected in the PC crest region but not in the WL head (fig. S9, A to C). RNAscope analysis on E11 head revealed that these molecules are highly expressed in PC crest feathers but are faintly expressed in the WL head (fig. S9, D to K). HoxC modulates feather phenotypes in adult WL scalp To test Hox gene function in adult feather follicles, we injected RCAS-HoxC8 to adult WL scalp feather follicles near the crown. Normally, head feathers are pennaceous with a long feather vane and a short after-feather growing from the same follicle. RCAS-HoxC8 perturbation of regenerating feathers produced a shorter feather vane and a longer after-feather (4/9) ([175]Fig. 5I). To test the involvement of HoxC family members in regulating feather morphology, we injected HoxC virus combinations (including HoxC5/C8/C9/C10/C11/C12), which result in expanding the plumulaceous region at the expense of the pennaceous region (5/20) ([176]Fig. 5J). This result suggests that adult feather follicles still have the plasticity to respond to altered HoxC and downstream gene expression, which may provide positional information. Thus, HoxC gene perturbations can alter fate determinations along both the anterior-posterior and distal-proximal axis in adult feather follicles. PC crest feathers express high levels of HoxC PC non-crest feathers are similar to WL short cranial feathers ([177]Fig. 1, A to B). The prominent long PC crest feathers occupy most of the scalp area, which dramatically diminishes its comb domain. The scalp structures of both chickens were determined using coronal sections with hematoxylin and eosin (H&E) staining (fig. S10, A to B). Higher HoxC8 and HoxC10 expression levels are detected by RNAscope in the PC crest feathers (fig. S10C, red arrows) compared to their adjacent non-crest feathers and their WL counterparts. Compared with non-crest feathers, strong HoxC8 expression appears in the pulp, dermal papilla (DP), and papilla ectoderm (the epithelium surrounding the DP), while high HoxC10 expression was localized within the pulp of PC crest feathers, suggesting that HoxC genes most likely function in dermal components. However, the HoxC8 expression pattern suggests that it may regulate both dermal and epidermal progenitor cells ([178]14). The small comb and high HoxC8 expression in the PC scalp echo the results from our HoxC8 modulation study, which changed the comb-base domain size in the WL scalp. Next, we examined HoxC cluster regulation differences between PC long crest feathers and non-crest feathers at transcriptomic levels. RNA-seq comparisons of juvenile PC crest and non-crest feathers found 1061 significantly up-regulated and 1265 down-regulated genes (fig. S10D). Canonical pathway analysis using these DEGs indicates that Wnt signaling is the most significantly enriched morphogenic pathway (fig. S10E), mimicking our HoxC8 overexpression result (fig. S8H). In the HoxC cluster, only HoxC8, HoxC9, and HoxC10 are identified as DEGs and are up-regulated in the crest. Wnt signaling pathway members (Wnt4, Wnt6, Wnt10A, LGR5, etc.), MSX2 and MAT1A, up-regulated in the HoxC8 over-expressing sample are also up-regulated in the PC crest and were confirmed by SISH (fig. S10, F and G). These data suggest that different HoxC expression profiles and Wnt signaling activity could contribute to distinct chicken scalp skin appendages (fig. S10H). Modulation of HoxC8/C10 expression in the PC regulates feather length We next asked whether modulating HoxC8/C10 expression in PC crest could mimic the CRISPR–195 bp result. To do this, we injected RCAS-dnHoxC8 or dnHoxC10 into PC head skin at E4. Both dnHoxC8 (N = 4/5) and dnHoxC10 (N = 3/3) transfected samples showed reduced feather growth in the crest region (fig. S11, B and C), compared to the normal PC control (fig. S11A). The phenotype resembles the type II CRISPR–195 bp result ([179]Fig. 4D). HoxC expression patterns in JAQ and CAQ suggest that their effect on scalp morphology is conserved across avian species Hox gene expression patterns are largely conserved during vertebrate development ([180]61). We next asked whether the HoxC gene expression domains might show conserved characteristics in the scalp of birds without combs [PC, Japanese quail (JAQ, with short crest feathers) and California quail (CAQ, with six elongated crest feathers)] ([181]Fig. 6A). WISH reveals a similar Hox gene (C8, C10, and C12) expression pattern within the trunk regions of these three species. Specifically, HoxC8 exhibits heightened expression in the anterior-back skin, HoxC10 displays pronounced expression in the posterior-back skin, and HoxC12 demonstrates elevated expression in the tail region (fig. S2, B to D′). However, diverse scalp HoxC expression patterns were present in these birds ([182]Fig. 6B). PCs show high HoxC8 and HoxC10 in the whole crest region (top row, red dotted line). JAQ faintly express HoxC8 and HoxC10 in scalp feather buds (second row), whereas CAQ express high HoxC8 levels but not HoxC10 in the elongated crest feather buds (third row). These results suggest that, compared to the trunk skin, the scalp is a more adaptable region that accommodates a variety of HoxC expression patterns that contribute to diverse skin appendage phenotypes. Fig. 6. Birds exhibiting exotic scalp feathers express diverse HoxC expression and contain the evolutionarily conserved 195 bp sequence. [183]Fig. 6. [184]Open in a new tab (A) PC, JAQ, and CAQ show different scalp feather phenotypes. Inset, enlargement of the head region. (B) WISH data show different HoxC expression patterns on the scalp of E9 PC, E8.5 JAQ, and E10 CAQ. The PC crest expresses high levels of HoxC8 and HoxC10 (outlined with a dotted red line). Insets shows strong HoxC8 but no detectable HoxC10 expression in CAQ crest feathers. (C) Phylogeny of tetrapods and evolution of HoxC10 and the 195-bp sequence. The phylogeny was constructed from the NCBI common tree using MEGA software ([185]86). The conservation analysis is based on UCSC multiple alignments (fig. S12). The 195-bp sequence is present in all birds and some reptiles but is absent in mammals. CAQ, California quail; JAQ, Japanese quail. The presence of 195 bp across different vertebrate species Last, we queried whether the 195-bp sequence was evolutionarily conserved. The 195-bp sequence localizes to the first intron of HoxC10 (fig. S12). We examined the DNA sequence conservation spanning exons 1 and 2 among 77 representative vertebrate species (cons77way) using UCSC multiple alignments. We found that HoxC10 exon 2 is conserved throughout the vertebrates examined. HoxC10 exon 1 is conserved in actinopterygian (zebra fish and stickleback fish) and tetrapod lineages, while the 195-bp and the flanking regions are present in turtles, crocodiles, and birds (i.e., archelosaurs) ([186]Fig. 6C and fig. S12). The results suggest that the 195-bp region evolved during the diversification of archelosaurs. Despite the importance of the HoxC cluster in mice ([187]31), the 195-bp sequence is absent in fish, amphibians, lizards, snakes, and mammals (fig. S12, blue asterisks), suggesting that chromatin loops within the HoxC cluster might form using alternative mechanisms to regulate expression in these species. DISCUSSION The higher-order chromatin organization of HoxC-containing TADs is skin region specific Hox genes are expressed in a colinear fashion and have been proposed to encode regional specificity in the spinal cord and limb ([188]32–[189]34). Here, we show that colinearity of HoxC genes along the AP axis of the chicken skin parallels chromatin accessibility. While the anterior limit of Hox expression in the central nervous system is at the level of the medulla oblongata in mice ([190]62), there is more anterior expression of HoxC genes in craniofacial regions of birds, likely mediated by an evolutionarily adaptive mechanism. HoxC is not expressed in the early development of the WL scalp because it is negatively regulated by a closed TAD in HoxC. This region will form a comb. In contrast, in the PC, there is high HoxC expression in the scalp, due to an open TAD in HoxC, resembling those in the back and tail skin of wild-type chickens that form feathers. These contact sites are 5 to 10 kb in size as determined by Micro-C and contain noncoding regions near sequences encoding potentially interesting genes including ATF7, Twist3, SP1, and PRKAG1 (table S1). These long-distance contact sites partially overlap in WL back, tail skin, and PC scalp skin. Using in vivo CRISPR technology, we removed the duplicated 195-bp region from the developing PC embryo scalp. This converted the long PC crest feathers to form short, disoriented feathers or apteric regions. Removing the duplicated 195 bp also reduced distant contact sites in the PC scalp. TADs affect gene expression by modulating the probability of enhancer-promoter interactions that control the spatiotemporal expression of developmental genes ([191]63–[192]66). Epigenetic landscapes can be changed by forming chromatin loops between TADs. For example, the HoxD cluster, important for tetrapod limb development ([193]67), is regulated by two flanking enhancer element containing TADs ([194]68). This suggests that enhancers and genes normally separated by TAD boundaries may contact each other as the TAD borders are lost, causing ectopic gene expression profiles ([195]49, [196]69). This enhancer adoption (or enhancer hijacking) has been shown to result in changes in chromatin configurations. The connection sites identified here in the skin help to open putative promoter regions within HoxC gene clusters. For example, in one of the contact sites, six TFs (Dlx6, MSX1, Dlx3, Lmx1b, En1, and Srebf2) may regulate the HoxC cluster. Future work is required to characterize these sites. Misexpression of HoxC genes results in homeotic transformations on the scalp of embryonic WL skin and alters feather forms in adults Our results show that dramatic transformations of skin and feathers can result from misexpressing members of the HoxC gene cluster. Similar effects of Hox gene expression have been shown in other tissue models. For example, Hox genes are involved in mammalian vertebral specification ([197]70, [198]71). Another example shows that a mutation in the HoxC region led to increased Wnt expression and the growth of long hairs in ear skin of Koala mice ([199]31). In this study, TAD disruption in the PC scalp enables promoter/enhancer interactions that increase HoxC8/C10 expression and result in the crested phenotype. Overexpressing HoxC members led to increased expression of Wnt, MSX2, etc., and overexpressing both HoxC genes and MSX2 resulted in feather formation in the WL embryo scalp, partially mimicking the crested PC phenotype. Overexpression of multiple HoxC members in the adult WL feather follicles adjacent to the comb did not alter the comb phenotype but did change feather forms, showing an altered rachis position and feather vane/fluffy region ratio. These results suggest that altering the HoxC expression profile in developing skin may result in changes of skin appendage types at the macroscale (feather versus comb). The adult region-specific epigenetic profile may be more stable, so changing the adult skin HoxC profile only caused microscale regional specific changes in feather morphology rather than macroscale regional specific changes of comb-feather transition. The evolution of feather regional specificity and potential roles of HoxC gene cluster in regional skin morphology We aim to delve deeper into the evolution of feather forms, their regional specificity, and functional significance. The main feather types are downy, contour, and flight feathers, serving endothermic insulation, display, and flight, respectively. Today, most birds show these feathers distributed in the ventral trunk, dorsal trunk, and wings, with scales on the feet to optimize function. However, achieving this regional specificity was not a simple process; rather, it was shaped by evolutionary selection. For instance, in ornithischian dinosaurs, we see filamentous appendages and scales, with some species like Tianyulong exhibiting a heavily feathered body and a scaled tail, while others, such as Psittacosaurus, display a largely scaled body with a partially feathered tail ([200]3). In non-avian theropods, which emerged later, simple branched feathers covered most of the body ([201]8). From this basic form, vaned feathers evolved, leading to complex flight feathers appearing on limbs and tails ([202]10). Thus, the evolution of feather phenotypes should be examined from two perspectives. The first is the branching complexity of individual feathers, progressing from radially symmetric to bilaterally symmetric, and eventually to bilaterally asymmetric forms ([203]7), a process controlled by morphogens like FGF, Wnt, BMP, and retinoic acid ([204]5, [205]6). Since feathers function as populations within tracts, the second essential perspective is the evolution of regional specificity, where groups of feathers are arranged to function collectively, such as those seen in the flight feather tract on the wing. The establishment of regional specificity may involve region-specific epigenetic profiles that regulate the spatial and temporal expression of downstream morphogens, resulting in optimized, localized functional forms. Although regional specificity in the integument is crucial, its regulatory mechanisms in amniotes remain largely unknown. In this study, we provide evidence that there is an anterior-posterior increasing HoxC gradient along the body and that HoxC gene clusters play a role in scalp specification. However, our findings also suggest that HoxC genes must interact with additional molecular factors to achieve regional specificity. Thus, this study represents a step toward uncovering the mechanistic control behind regional specificity of integumentary appendages. Compared to the trunk, the scalp seems to be a more flexible region in which to grow different crest feather types as seen in different birds such as PC, CAQ, and hoopoe. This could be due to the absence of other redundant Hox members seen in the trunk skin. It also implies that HoxC downstream pathways can be coupled to different morphogenetic modules to form diverse feather types. PC crest feathers are HoxC8^+ and HoxC10^+, while CAQ are HoxC8^+ but HoxC10^− ([206]Fig. 6, A and B). Other genes may also be involved in the crest phenotype. A crested duck study suggested the involvement of HoxC8, EphA2, EphA3, and EphB2 ([207]72), and a crested pigeon study found EphB2 as a strong candidate gene for crest development ([208]73). Thus, multiple factors in the putative crest feather–forming pathways may be involved in producing different crest feather morphologies ([209]74). Previous studies suggest that it is important to maintain TADs and that there is an evolutionary preference for TAD shuffling over TAD destruction ([210]50). Given that the 195-bp intron is functionally important, this sequence would be expected to be stable. This raises the intriguing question of whether the 195 bp–mediated HoxC regulation is only present in chicken or a trait in other species as well. If it is functionally important in HoxC gene regulation, the sequence should be stable in evolutionary lineages. Our comparative genomic analysis reveals that the entire first intron of HoxC10 is conserved in Archelosauria (avian, alligator, and turtle) but not in other vertebrates (lamprey, amphibians, lizards, and mammals) ([211]Fig. 6C). Thus, this intronic region may have evolved from a common Archelosauria ancestor. We speculate that the presence of the 195-bp sequence may affect the regulation of HoxC cluster genes, conferring more flexibility, thus allowing the generation of a large spectrum of diverse feather phenotypes and regional specificity during the evolution of feathered dinosaurs and Mesozoic birds. Perturbation of this 195-bp sequence, as seen in the PC here, leads to changes of integumentary appendage phenotypes. Mammals diverged from reptiles/birds more than 300 million years ago before the appearance of this 195-bp sequence. Birds/reptiles and mammals may use different mechanisms to achieve higher-order HoxC regulation. Summary, limitation, and future directions These results show that in WL scalp, a TAD forms a barrier that blocks HoxC gene expression by restricting interactions between promoters and distal enhancers residing outside of the TAD. The HoxC10 195-bp intron is unique in archelosaurs (turtles, crocodiles, and birds). It provides a clue that an enhancer adoption mechanism may be involved in regulating high-order chromatin interactions of the HoxC gene cluster with distant regulators and increase diverse skin appendage phenotypes. Thus, the current study inspires us to propose a working model that integumentary regional specificity may have evolved through a novel mechanism to couple region-specific epigenetic profiles with distinct morphogenetic modules. The region-specific epigenetic profiles involve the HoxC cluster, but the HoxC is not sufficient and additional loop anchor points are required and will be further investigated. The interactions of the region-specific epigenetic profiles with distinct morphogen modules led to the formation of specific feather types in different body regions or different phenotypes in different species (e.g., unique appendages on the scalp of CAQ, cockatiel, and red-crowned crane). This study opens a frontier for future studies on the evolution of regional feather tracts, on how loop anchor sites regulate the HoxC gene cluster, and how the morphology of each feather is regulated at the epigenetic level. MATERIALS AND METHODS Ethics statement All the animals studied in this work were raised and handled according to an authorized protocol (21499) of the Institutional Animal Care and Use Committees of the University of Southern California (USC, Los Angeles, CA). Eggs and chickens Pathogen-free fertilized WL chicken eggs, JAQ eggs, and CAQ eggs were purchased from Charles River Laboratories (CT, USA), AA laboratories (CA, USA), and Murray McMurray Hatchery (IA, USA), respectively. The PC eggs were obtained from Meyer Hatchery (OH, USA) and Frizzled & Fabulous Farm (CA, USA). All eggs were incubated at 37.5°C. Adult chickens were hatched from the eggs. All chickens were raised at 22°C in the USC Department of Animal Resources. Embryonic skin and feather follicle dissection Scalp, neck, anterior back (between the wings), posterior back (between the thighs), and tail skins were dissected from E9 chicken embryos under a stereo microscope and kept in phosphate-buffered saline (PBS) on ice. To harvest the dermis part, the skins were treated with calcium-magnesium–free solution (2 × CMF) solution with 0.25% EDTA for 15 to 25 min on ice, and afterward, the epidermis was removed. For adult feather follicles, growth-phase feathers regenerated for 4 weeks after plucking with intact DP were collected and immersed in PBS on ice. The proximal 0.5-cm-long follicles were kept for further dissection. The follicles were cut open by a micro-scissor under a dissecting microscope and submerged in 2 × CMF containing 0.25% EDTA for 20 min on ice. Dermal components (DP and pulp) were then harvested by separating them from epidermis using forceps. Cell dissociation Dissected dermal samples were treated with 0.1% trypsin (Life Technologies, catalog. 27250-042) and 0.1% collagenase type I (Worthington, catalog [212]LS004196) in PBS at 37°C for 30 to 60 min with gentle rotation. Stiff feather follicle samples were further cut into small pieces by scissors, and the cells/tissues were gently pipetted up and down to facilitate disaggregation every 5 min. Dulbecco’s modified Eagle’s medium (3× volume) supplemented with 10% fetal bovine serum, 2% chicken serum, and 1% penicillin-streptomycin was added to stop the collagenase digestion. The cell suspension was mixed by gently pipetting up and down and filtered through a 70-μm and then a 40-μm cell strainer. After the dissociated cells were pelleted by centrifugation at 200g for 15 min at room temperature (RT), the samples were immediately placed on ice. Supernatants were removed, the cell pellets were resuspended and washed with cold PBS, and cell numbers were calculated. The harvested cells were used for RNA extraction, ATAC-seq, and Micro-C experiments. Total RNA extraction and RNA-seq Each group contained two biological replicates. Total RNAs were isolated using TRIzol reagent (Invitrogen, catalog 15596026). RNA quantity and purity were measured by a NanoDrop, and RNA integrity was analyzed using agarose gel electrophoresis (18S and 28S rRNA band ratio) as well as Agilent 2100 Bioanalyzer. For E9 dermal samples, 1 to 2 μg of total RNA from each sample was used to construct an RNA-seq library using TruSeq RNA sample preparation v2 kit (Illumina) and sequencing (50-cycle single read) was performed using Hi-seq2000 by the USC Molecular Genomics Core. Other RNA-seq library preparation and sequencing were conducted at Novogene Corporation (CA, USA). These libraries were generated using NEBNext Ultra RNA library prep kit for Illumina (NEB, USA) and then sequenced on the Illumina NovaSeq 6000 platform (150-cycle paired-end reads). RNA-seq analysis The RNA-seq analysis was performed on the Partek Flow platform (Partek Inc.) using default parameters. The Ensembl galGal6 and annotation were used. Mapping using STAR-2.7.3 was followed by post-alignment QA/QC (table S2). The aligned reads were quantified to annotation model (Ensembl release 103) using Partek E/M method. For plotting heatmaps and bar charts, the raw gene counts were normalized with TPM. DEG analysis To reduce noise, the genes with low aligned read counts were filtered out until around 80% of the genes passed the filter. Then, the filtered gene counts were normalized with median ratio. DESeq2 ([213]75) was used to identify DEGs that were filtered with |FC| >2 and either P value ≤0.05 or FDR ≤0.05. GO analysis to determine pathway enrichment was conducted using Ingenuity Pathway Analysis software (QIAGEN Inc.). Principal components analysis (PCA) was used to determine sample consistency. Volcano plots were generated in Partek Flow. The deepTools ([214]76) bamCoverage was applied to convert Bam files to bigWig files for visualization with the Integrative Genomics Viewer (IGV) ([215]77). Assay for transposase-accessible chromatin using sequencing Two biological replicates of each group were used. All procedures were performed on ice. ATAC resuspension buffer (ARB) was composed of 3 mM MgCl[2], 10 mM NaCl, and 10 mM tris-HCl (pH 7.4). Dissociated cells (50,000) of indicated tissues were lysed in 50 μl of cold ATAC lysis buffer (0.1% NP-40, 0.1% Tween-20, and 0.01% digitonin freshly added in ARB) for 5 min. One milliliter of ARB (with 0.1% Tween-20) was added and mixed well by inverting the tube. Nuclei were pelleted by 2500 rpm at 4°C for 10 min, and supernatants were discarded. Tn5 transposase (2.5 μl) was added to tagment chromatin at 37°C for 30 min in 50 μl of 1× Tagment DNA buffer (Illumina, catalog. FC-121-1030) with 0.1% Tween-20, 0.01% digitonin, and 0.33× PBS. The tagmented DNA was purified using a DNA Clean & Concentrator-5 kit (Zymo Research, catalog. D4014). To prepare the library, the tagmented DNA fragments were amplified in a thermocycler for 5 cycles using NEBNext High-Fidelity 2× PCR Mater Mix (NEB, catalog. M0541S) in a total volume of 50 μl. Five microliters of the PCR products was used in qPCR to determine the number of extra amplification cycles before saturation, corresponding to 25% of maximum fluorescent intensity. The remaining 45 μl of the PCR products was amplified using the previously calculated cycle numbers and purified using 1× AMPure XP beads (Beckman, A63880). For sequencing, the libraries were submitted to Novogene Corporation (CA, USA) and quantitated and analyzed using Qubit and Agilent 2100 Bioanalyzer, respectively. The libraries passing the quality control were sequenced in the Illumina NovaSeq 6000 system using 150-cycle paired-end reads. ATAC-seq analysis The ATAC-seq analysis was conducted on Partek Flow platform (Partek Inc.). Default parameters were applied except for those with additional notes. Prealignment QA/QC algorithms were run first. Mapping BWA-MEM algorithm v.0.7.17 (Genome: Ensembl galGal6) was adopted and then analyzed by post-alignment QA/QC (table S2). Filter alignment Aligned reads were further filtered to remove duplicate alignments (defined as the same start and same sequence). Peak-calling MACS2 - 2.1.1 (peak cutoff q value < 0.05; no shifting model; shift size, 37; extension size, 73) was used. The identified peaks were quantified (minimum read overlap with feature 80%), normalized with TPM, and annotated with TSS upstream limit 2000, TSS downstream limit 0, TTS upstream limit 0, and TTS downstream limit 300. DAR analysis Gene set analysis was used to detect DARs that were filtered with |FC| >2 and P value ≤0.05. GO term enrichment analyses were performed using PANTHER classification system ([216]78) and KEGG database, respectively. Motif prediction was determined using JASPAR database ([217]79). PCA plots were used to validate replicate consistency of each group. Bigwig files were used to visualize data in IGV. Whole-mount in situ hybridization RNA probes used in WISH were generated by PCR reactions from an E9 chicken skin cDNA template. The primers are listed in table S3. The PCR fragments were purified and inserted into the pDrive vector (Qiagen). Digoxigenin-labeled probes were transcribed using SP6 or T7 RNA polymerase (Roche). WISH was conducted as described ([218]21). Paraffin section, SISH, IM, and RNAscope Paraformaldehyde fixed chicken tissue samples were embedded in paraffin blocks. Seven-micrometer paraffin sections were prepared for staining with H&E, ISH, immunostaining (IM), and RNAscope according to the procedures described before ([219]14). HoxC8 (Proteintech, 15448-1-AP), AMV-3C2 (Developmental Studies Hybridoma Bank), and proliferating cell nuclear antigen (DAKO, clone PC10) antibodies were used for IM. The Multiplex Fluorescent v2 system (ACD, 323100) and the manufacturer’s instructions were adopted for RNAscope. We used the following probes: HoxC6 (1079941-C1), HoxC8 (1079951-C2), HoxC10 (1079961-C3), LGR4 (1097771-C2), LGR5 (480781-C1), LGR6 (1097781-C3), SFRP4 (1173691-C1), FN (489431-C2), SHH (551581-C3), MSX1 (480091-C2), MSX2 (1593561-C2), and ATF7 (1593571-C3). Construction of RCAS virus The target gene inserts were synthesized by PCR reactions using cDNA of E7 embryonic chicken skin as template. For HoxC8 cloning, PCR was conducted using primer set (forward, 5′-aaatatgcggccgcaccatgagttcctactttgtaaatc; backward, 5′-tgcactagtgtccttgctttcttcttttt) to generate the chicken HoxC8 CDS. The 5′-aaatatgcggccgcaccatgagttcttacgtcgcca (forward) and 5′-tgcactagtaagagcctctttgcttttcaa (backward) primers were used for the HoxC5 CDS, 5′-aaatatgcggccgcaccatgtcggcttctggccccat (forward) and 5′-tgcactagttggctgttccttatcggttttctc (backward) were for HoxC9 CDS, 5′-aaatatgcggccgcaccatgtcatgtcccaacaatg (forward) and 5′-tgcactagtggtgaaattaaaattggacgtca (backward) were for the HoxC10 CDS, 5′-aaatatgcggccgcaccatgtttaactcagtgaatctgg (forward) and 5′-tgcactagtcaataagggatttccagagaaata (backward) were for HoxC11 CDS, 5′-aaatatgcggccgcaccatgggcgagcacaacctc (forward) and 5′-tgcactagtaaagaaagacagggcttgctcc (backward) were for the HoxC12 CDS, and 5′-aaatatgcggccgcaccatggcttctccttccaaag (forward) and 5′-tgcactagtagataagtggtacatgctatatcc (backward) were for the MSX2 CDS. These gene inserts were ligated into the RCAS vector carrying a mCherry fluorescent tracer (RCAS-2A-mCherry). RCAS viruses were prepared according to the protocol described ([220]2) and concentrated at 20,000 rpm for 2 hours at 4°C. For dominant-negative construct, primers (forward, 5′-aaatatgcggccgcaccatgagttcctactttgtaaatc; backward, 5′-tgcactagtccagatcttcacttgtctct) was used to produce dominant-negative form of HoxC8 (RCAS-dnHoxC8), which is a truncated HoxC8 missing a C-terminal part of the homeodomain ([221]80). Another set of primers (forward, 5′-aaatatgcggccgcaccatgtcatgtcccaacaatg; backward, 5′-tgcactagtccagattttgacttgtctgt) was used to yield the dominant-negative form of HoxC10. To make a RCAS-dominant negative ATF7 (RCAS-dnATF7), we designed a chicken mutant form based on a dominant-negative human mutation in which a mitotically nonphosphorylatable version of ATF7 acts as a dominant-negative protein ([222]58). The cloning procedure includes two PCR steps and the following Fusion step. For the first PCR, the primer set is 5′-ACTCTGCTGGCGGCCGCGCCACCatgggcgacgacagaccctt-3′ (forward) and 5′-gcagttcttcaggaaccgagtcggcgccggggcctgatc-3′ (backward). For the second PCR, the primer set is 5′-actcggttcctgaagaactgcgaggaggt-3′ (forward) and 5′-CTCTGCCCTCACTAGTTctgcccgcagactgcgact-3′ (backward). The two PCR products and RCAS-2A-mCherry backbone are fused using the in-fusion cloning method (Takara, catalog 638947) to generate RCAS-dnAFT7. Gene misexpression For gene misexpression in embryos, 1 μl of concentrated RCAS virus was injected to the HH24 chicken embryo scalp skin and the embryos were harvested at HH39. RCAS-2A-mCherry without an insert was used as a control. Five replicates were harvested for each condition. For gene misexpression in a juvenile chicken, 1-month-old male scalp feathers were plucked and RCAS virus (100 μl) was injected into the empty feather follicles. Feathers were regenerated for 50 days before collection. At least nine feathers were harvested for each condition. Micro-C sample preparation A Dovetail Micro-C Kit (Dovetail, catalog 21006) was used for Micro-C experiments. The protocol was modified from the Dovetail Micro-C kit user guide (version 1.0). Cross-linking, digestion, and lysis Dissociated cells (5× 10^6) from tissues were pelleted down in a 15-ml tube at 200g for 10 min at 4°C. The supernatant was discarded, and the cell pellet was frozen at −80°C for at least one night. The cell pellet was thawed at RT and then resuspended in 1 ml of PBS. Ten microliters of fresh 0.3 M Disuccinimidyl Glutarate (DSG; Thermo Fisher Scientific, catalog A35392) in dimethyl sulfoxide was added, and the tube was rotated for 10 min at RT. Twenty-seven microliters of fresh 37% formaldehyde (Sigma-Aldrich, F8775-4X25ML) was added, and the tube was rotated for 10 min at RT. The cross-linked cells were pelleted at 3000g for 5 min, and then the supernatant was discarded. The cell pellet was washed with 200 μl of 1× wash buffer by gently pipetting to resuspend the pellet, the cells were spun down at 3000g for 5 min, and the supernatant was discarded. The wash step was repeated with 400 μl of 1× wash buffer, 2 × 10^6 cells were centrifuged in a 1.5-ml tube, and the cell pellet was resuspended in 50 μl of fresh 1× nuclease digest buffer. Micrococcal nuclease (0.02 to 0.5 μl) enzyme mix was added, and the tube was incubated in a thermal mixer at 22°C for 30 min at 1250 rpm. Then, 5 μl of 0.5 M EGTA was added to stop the reaction, and 3 μl of 20% SDS was added to lyse the cells, followed by incubation at 22°C for 5 min at 1,250 rpm. Lysate QC A total of 2.5 μl of the lysate was used for QC test, while the remaining 56 μl of the lysate was stored at −80°C for the next stage, proximity ligation. The 2.5-μl lysate was mixed with 51.5 μl of cross-link reversal master mix (1× cross-link reversal buffer with 1.5 μl of Proteinase K) and then incubated in a thermal cycler for 15 min at 55°C, 45 min at 68°C, and held at 25°C. The QC DNA was purified using a DNA Clean & Concentrator-5 kit (Zymo Research, catalog D4014) and eluted in 10 μl of DNA elution buffer. The quantity and fragment size distribution of the purified QC DNA were measured by a NanoDrop and Bioanalyzer (HS DNA kit), respectively. Proximity ligation One microgram of the lysate was mixed with 100 μl of resuspended Chromatin Capture Beads and then incubated at RT for 10 min. The tube was placed in a magnetic rack for 5 min, and the supernatant was discarded. The beads were washed with 150 μl 1× wash buffer and placed on the magnetic rack for 3 min. The supernatant was discarded. The wash step was repeated once. To optimize DNA termini, 53.5 μl of End Polishing Master Mix was added to the beads and incubated in the thermal mixer at 1250 rpm for 30 min at 22°C, followed by 30 min at 65°C. The tube was placed in the magnetic rack for 5 min at RT, and the supernatant was discarded. The wash step was repeated. To perform bridge ligation, 50 μl of Bridge Ligation Mix and 1 μl of Bridge Ligase were added, and the tube was incubated in the thermal mixer at 22°C for 30 min at 1250 rpm. The supernatant was discarded by pelleting beads in the magnetic rack. The wash step was repeated once. Fifty-two microliters of Intra-Aggregate Ligation Master Mix was added to the beads, and the tube was incubated in the thermal mixer at 22°C overnight at 1,250 rpm. The supernatant was discarded in the magnetic rack. A total of 51.5 μl of cross-link reversal master mix (1× cross-link reversal buffer with 1.5 μl of Proteinase K) was added and then incubated in a thermal mixer for 15 min at 55°C, 45 min at 68°C, and held at 25°C. The beads were pelleted in the magnetic rack, and 50 μl of the supernatant was transferred to a new 1.5-ml tube. The beads were discarded. To purify the DNA, SPRIselect beads were used. Ninety microliters of the resuspended SPRIselect beads were mixed with 50 μl of the sample by vortexing. The tube was quickly spun down and incubated at RT for 15 min. The supernatant was discarded in the magnetic rack. The beads were washed twice with 200 μl of fresh 80% ethanol without bead resuspension. The beads were air dried for 10 to 15 min in the magnetic rack. The beads were resuspended in 52 μl TE buffer (pH 8.0) and incubated at RT for 5 min. In the magnetic rack, 50 μl of the supernatant (purified DNA) was transferred to a new tube. The beads were discarded. The purified DNA was quantified by a Qubit fluorometer with Qubit double-stranded DNA (dsDNA) HS kit. Library preparation To repair the DNA termini, 50 μl of purified DNA (150 ng) was mixed with 10.5 μl of End Repair Master Mix in a PCR tube and then incubated in a thermal cycler for 30 min at 20°C, 30 min at 65°C, and held at 12°C. To ligate adapters to both ends of DNA, 31 μl of Ligation Enzyme Master Mix and 2.5 μl of adaptor for Illumina were added and incubated in the thermal cycler for 15 min at 20°C and held at 12°C. Three microliters of USER Enzyme Mix was added and incubated in the thermal cycler for 15 min at 37°C and held at 12°C. Eighty microliters of resuspended SPRIselect beads were mixed with the 97 μl of the sample to purify the DNA. The purification steps mentioned above were followed, the DNA was eluted with 100 μl of TE buffer (pH 8.0), and 95 μl of the supernatant (purified adaptor-ligated DNA) was transferred to a new tube. The beads were discarded. Ligation capture and amplification To prepare streptavidin beads, 25 μl of the resuspended streptavidin beads were transferred to a new 1.5-ml tube and pelleted in the magnetic rack, and the supernatant was discarded. The following buffer names are from the commercial Dovetail Micro-C kit version 1.2 (Cantata Bio). The beads were washed twice with 200 μl of TWB Solution and then resuspended in 100 μl of 2× NTB Solution. The 95-μl purified adaptor-ligated DNA was added to the beads, mixed well, and incubated in the thermal mixer at 25°C for 30 min at 1250 rpm. The magnetic rack was used to do a series of washes. The beads were washed once with 200 μl of LWB Solution, twice with 200 μl of NWB Solution, and twice with 200 μl of 1× wash buffer. Forty-five microliters of Index PCR Master Mix and 5 μl Index Primer were added and incubated in the thermal cycler for 3 min at 98°C, 12 cycles of amplification (20 s at 98°C, 30 s at 65°C, and 30 s at 72°C), 1 min at 72°C, and held at 4°C. To perform size selection, in the magnetic rack, 47 μl of the supernatant was transferred to a new 1.5-ml tube, and 53 μl TE buffer (pH 8.0) was added to bring the volume to 100 μl. Fifty microliters of resuspended SPRIselect beads was added and incubated at RT for 10 min. In the magnetic rack, the 145 μl of supernatant was transferred to a new tube. Thirty microliters of the resuspended SPRIselect beads was added and incubated at RT for 10 min. The supernatant was discarded in the magnetic rack. The beads were washed twice with 200 μl of fresh 80% ethanol, air dried for 5 min, resuspended in 30 μl TE buffer (pH 8.0), and incubated at RT for 5 min. In the magnetic rack, the 28 μl supernatant (size selected library) was transferred to a new tube. The libraries were quantified by a Qubit fluorometer with Qubit dsDNA HS kit. The fragment size distribution of the libraries was analyzed by Bioanalyzer (HS DNA kit). The ideal size range is between 350 and 1000 bp. Sequencing The libraries were shipped to Novogene Corporation (CA, USA) and sequenced on the Illumina NovaSeq 6000 platform with paired-end 150 cycles. Micro-C analysis Micro-C data were preprocessed according to Dovetail Micro-C Data Processing Guide. Briefly, BWA-MEM was used to map the sequencing reads to chicken galGal6 genome. Valid pairs from the mapped reads were identified by pairtools parse pipeline. Unmapped, low mapping quality (quality <40), and PCR duplicate read pairs were removed. Library QC and library complexity were shown in table S2. To visualize Micro-C data, juicertools.jar and Cooler tools were used to generate 1-kb resolution hic and cool contact matrices, respectively. Hic files can be visualized on Juicebox ([223]81). Cool files were visualized using the Galaxy platform ([224]82), version 22.1. The snapshots of contact matrices shown in this study were produced by Galaxy. TADs were identified by juicertools.jar arrowhead with a resolution of 10 kb. Loops were discovered by Mustache ([225]83) with 5-kb resolution and an FDR smaller than 0.2. Tol2-CRISPR assay To build a high efficiency Tol2-CRISPR-CTNNB1 plasmid, we modified the single CRISPR plasmid that targets chicken CTNNB1 ([226]54) by inserting a Tol2 autonomous transposon. To do this, we used BamHI and SpeI restriction enzymes to remove the Amp gene from CRISPR-CTNNB1. pT2AL200R175-CAGGS-LynEGFP ([227]84) was used as a PCR template to amplify the region between Tol2-F and Tol2-R regions. The primer set is Tol2-F (5′-tcgtgccagcggatcCagatctAATACTCAAGTACAATTTTAATGGAGTAC-3′ forward) and Tol2-R (5′-actattaataactagtcaataatcaatgtcCgggcccAAGTGATCTC-3′ backward). The two fragments are linked by in-fusion cloning method (Takara, catalog 638947) to generate the Tol2-CRISPR-CTNNB1 plasmid. We used a website [228]benchling.com to design CRISPR targeting the 195-bp duplicated region in the PC genome. The target sequence (CTGTGTGCGTTTGAAGACGT) was cloned to the single CRISPR plasmid according to the published method ([229]54). After the cloning step 5 [Amplification of HH-gRNAf+e-HDV precursor in ([230]54)], the product was cut with EcoRI/NotI and ligated with EcoRI/NotI-digested Tol2-CRISPR-CTNNB1 plasmid to generate the final product (Tol2-CRISPR–195 bp). A sequence targeting mCherry gene (CGCCCTCGATCTCGAACTCG) was used as CRISPR control (Tol2-CRISPR-Cherry). Tol2-CRISPR-CTNNB1 in vitro testing To test the efficiency of the Tol2-CRISPR-CTNNB1 construct, we added the plasmid (2 μg/μl) to DF1 cells [plasmid DNA (24 ug) and Lipofectamine 2000 (60 μl) in 10-cm dish] in vitro culture for 2 days. Testing efficiency was determined using the T7 endonuclease assay ([231]85). Densitometry analysis was conducted with ImageJ. Tol2-CRISPR-CTNNB1 in vivo testing We injected 1 μl of CRISPR plasmid (2.5 μg/μl) with the tol2 transpose (T2TP) plasmid (2.5 μg/μl) to the HH16 WL somite. Four somites on the left side were injected. The samples were collected at HH35. Tol2-CRISPR–195 bp in vitro testing E10 PC crest fibroblast cells were cultured until 90% confluency. The cells were electroporated by using Neon Transfection System (1200 V/20 ms/2 pulses) with the Tol2-CRISPR–195 bp plasmid (1.2 μg/μl) and cultured for 2 days. Green fluorescent protein–positive cells were sorted, and DNA was extracted. PCR with primers flanking the duplicated 195-bp region was performed [Forward (5′-TCCTTGGCGAATTGAAGTCT-3′) Reverse (5′-TTCGATCCGCTTAAATCTGC-3′)]. The PCR products were cloned to pDrive (Qiagen), and 10 colonies were collected for sequencing. The Tol2-CRISPR-Cherry plasmid (1.2 μg/μl) was used as a negative control for CRISPR activity. Tol2-CRISPR-195 bp in vivo testing We injected 0.5 μl of CRISPR plasmid (2.5 μg/μl) with T2TP plasmid (2.5 μg/μl) together to HH11 PC neural crest. The injection was followed by electroporation with three pulses of 15 V/50 ms. The samples were collected at HH39. Real-time qPCR The skin was dissected from the region between the beak and the top of the crest region from treated (N = 7) and control samples (N = 4). TRIzol reagent (Invitrogen, 15596018) was used to extract RNA, and All-In-One 5X RT MasterMix (abm, G592) was used to generate cDNA. Primers for qPCR are HoxC6-F (5′-CCAGTGACCAAAACAGGAACAC-3′), HoxC6-R (5′-GGTAACGGGAATAAATCTGGCG-3′), HoxC8-F (5′-GTCTCCCAGTCTCATGTTCCC-3′), HoxC8-R (5′-GGGCGTGAGAGACCTCAATC-3′), HoxC10-F (5′-GCACCCAAAATCAGCCCTTC-3′), HoxC10-R (5′-TTCTTCCGCTCTTTGCTGTCA-3′), glyceraldehyde-3-phosphate dehydrogenase F (GAPDH-F) (5′-CGATCTGAACTACATGGTTTACATG-TT-3′), and GAPDH-R (5′-CCCGTTCTCAGCCTTGACA-3′). qPCR was performed using Maxima SYBR Green/ROX qPCR Master Mix (Thermo Fisher Scientific, K0221) on the QuantStudio 3 Real-Time PCR System (Applied Biosystems, [232]A28567). The samples were incubated at 50°C for 2 min, 95°C for 10 min, followed by 40 cycles of 95°C for 15 s, 55°C for 30 s, and 72°C for 30 s. The delta-delta CT method was used to calculate the fold gene expression. Statistical analysis To plot heatmaps for [233]Fig. 1 (D and G) and fig. S8H, standardized gene expression values (z scores) were calculated from [TPM − (mean TPM)/SD]. In [234]Figs. 4G and [235]5D, data are means ± SD. Significance was determined by a paired two-tailed Student’s t test. Acknowledgments