Modeling uniquely human gene regulatory function via targeted humanization of the mouse genome – Nature.com

Generating a HACNS1 knock-in mouse model

We designed a targeting construct for homologous recombination including a 1.2kb human sequence encompassing HACNS1 that was previously shown to encode human-specific enhancer activity in transgenic mouse embryos10. We replaced the orthologous mouse locus using homology-directed repair in C57BL6/J-AwJ/J (B6 agouti) embryonic stem (ES) cells (Fig.1B, Supplementary Figs.1A, B, Supplementary Data1; Methods). To provide a control that would enable us to distinguish bona fide human-specific functions of HACNS1 from possible primate-rodent differences, we also generated a mouse model for the orthologous chimpanzee sequence using the same approach (Supplementary Figs.1AC). The 1.2kb chimpanzee sequence shows no evidence of evolutionary acceleration and includes 22 single nucleotide differences relative to the human sequence3; 15 of these are human-specific based on comparisons to other primate genomes (Methods, Supplementary Data2). Previous studies indicate that multiple human-specific substitutions contribute to the gain of function in HACNS110. Twelve of the 15 substitutions introduce one or more predicted transcription factor binding sites to the human sequence (Supplementary Data2, Supplementary Fig.1A). An extensive comparison of sequence conservation and divergence among the human, chimpanzee, and mouse sequences is provided in theSupplementary Note (Supplementary Materials).

In order to verify the integrity of the edited loci, we sequenced a 40kb region encompassing the human or chimpanzee sequence replacement, the homology arms used for targeting, and flanking genomic regions in mice homozygous for either HACNS1 or the chimpanzee ortholog (Supplementary Fig.1C; Methods). We found no evidence of aberrant editing, sequence rearrangements, or other off-target mutations at either edited locus. We also verified that each homozygous line carried two copies of the human or chimpanzee sequence using quantitative real-time PCR (RT-qPCR) (Supplementary Fig.1D; associated Source Data).

We used chromatin immunoprecipitation (ChIP) to determine if HACNS1 exhibits epigenetic signatures of increased enhancer activity in mice. We first performed epigenetic profiling in the developing mouse limb bud based on prior evidence that HACNS1 drives increased reporter gene activity in transgenic enhancer assays and exhibits increased H3K27ac marking in the human embryonic limb10,11. We profiled both H3K27ac and H3K4 dimethylation (H3K4me2), which is also associated with enhancer activity, in embryonic day (E) 11.5 limb buds from embryos homozygous for HACNS1, embryos homozygous for the chimpanzee ortholog, and wild type embryos. We found a strong signature of H3K27ac marking at HACNS1 in the limb buds of HACNS1 homozygous embryos (Fig.2, Supplementary Fig.2). The chimpanzee and mouse sequences both showed significant but weaker H3K27ac enrichment relative to the human sequence, supporting the conclusion that HACNS1 maintains its human-specific enhancer activity in the mouse genomic and developmental context.

Epigenetic profiling in the HACNS1 homozygous E11.5 limb bud compared to the chimpanzee ortholog line and wild type. The normalized H3K27ac signals are shown for the HACNS1 line (in dark green), the chimpanzee ortholog line (in olive), and wild type (in teal) (see Methods). The location of the edited HACNS1 locus in the human ortholog line relative to nearby genes is shown above the track. The double slanted lines indicate intervening H3K27ac signal data between the edited and wild type loci and Gbx2 that were removed for clarity; see Supplementary Fig.2 for complete views for each line as well as input signals. H3K27ac peak calls showing significant increases in signal between HACNS1 homozygous and wild type, and the corresponding peak regions compared between the chimpanzee control line and wild type, are shown below the signal track. Litter-matched embryos were used for each comparison (see Methods). N.S. not significant. All peak calls for each line are shown in Supplementary Fig.2. Adjusted P values were obtained using DESeq2 (implemented in HOMER) with a Wald test followed by Benjamini-Hochberg correction23,24.

We used DESeq2 implemented in HOMER (Methods) to identify genome-wide significant differences in H3K27ac and H3K4me2 levels in E11.5 limb buds from mice homozygous for HACNS1 or the chimpanzee ortholog versus wild type23,24. We found that H3K27ac and H3K4me2 levels were significantly increased at the edited HACNS1 locus compared to the endogenous mouse locus (Fig.2, Supplementary Fig.2A, B,Supplementary Data3). In contrast, the level of H3K27ac at the edited chimpanzee locus was not significantly different than that at the endogenous locus (Fig. 2,Supplementary Fig.2A, Supplementary Data3). The levels of H3K4me2 were significantly increased at both the humanized and orthologous chimpanzee loci in each respective line compared to mouse (Fig.2, Supplementary Fig.2, Supplementary Data3). As high levels of H3K4me2 coupled with low levels of H3K27ac are associated with weak enhancer activity25, it is likely that the chimpanzee sequence is not acting as a strong enhancer in the limb bud overall, a finding further supported by the gene expression analyses described in Fig.3.

A Spatial and temporal expression of Gbx2 in HACNS1 homozygous, chimpanzee ortholog line, and wild type E11-E12 embryos visualized by whole-mount in situ hybridization (ISH). Representative images are shown for each genotype at three fine-scale time points; see text and associated Source Data for details on staging. Magnified views of Gbx2 expression in limb buds are shown to the right of each embryo. Annotations of anatomical structures and developmental axes are indicated at the top right: FL forelimb, HL hindlimb, DI diencephalon, NT neural tube, PA pharyngeal arch, A anterior, P posterior. The arrows at the top far right indicate the anterior-posterior (A-P) and proximal-distal (Pr-D) axes for the magnified limb buds. Bottom right: Crown-rump lengths for all embryos assayed for Gbx2 mRNA expression by ISH. Each point indicates a single embryo. Colors denote each fine-scale time point (T1-T6). B Left: representative images of anterior, posterior, proximal, distal (top), and strong versus weak Gbx2 staining patterns (bottom). Anterior (A), posterior (P), and body wall (BW) domains are denoted on top left limb bud. Right: Gbx2 ISH staining pattern data across 6 developmental timepoints from each of three independent, blinded scorers (marked at top as counting replicates 13; see text and Fig. 3A for timepoint scheme and associated Source Data for annotations). The darkest shade for HACNS1 homozygous (dark green), chimpanzee ortholog line (olive), and wild type (teal) represents percentage of forelimbs or hindlimbs showing strong anterior and posterior limb bud staining. Medium-dark shade, as shown inthe legend on the left, denotes strong anterior staining only, while the lightest shade denotes weak staining in any domain. Independent biological specimens were analyzed for n=139 (wild type), n=103 (chimpanzee ortholog line), and n=106 (human ortholog line) forelimbs and n=137 (wild type), n=103 (chimpanzee ortholog line), and n=102 (human ortholog line) hindlimbs. For body wall and pharyngeal arch scoring data see Supplementary Fig.3A, B and associated Source Data.

Previous transgenic mouse enhancer assays indicated that HACNS1 drives increased reporter gene activity in the pharyngeal arch compared to the chimpanzee ortholog10. We therefore profiled H3K27ac and H3K4me2 in pharyngeal arch tissue from E11.5 embryos homozygous for either HACNS1 or the chimpanzee ortholog. We detected reproducible, significant enrichment of H3K27ac in the pharyngeal arch at the humanized and orthologous chimpanzee loci compared to input controls, but the H3K27ac signal at neither the human nor the chimpanzee ortholog locus was significantly different compared to the mouse endogenous locus (Supplementary Fig.2). We did, however, identify a significant gain of H3K4me2 signal in the pharyngeal arch at the humanized and orthologous chimpanzee loci compared to the mouse locus (Supplementary Fig.2, Supplementary Data3).

In order to identify downstream epigenetic changes resulting from HACNS1 activation in the limb bud, we searched for other genome-wide gains of H3K27ac and H3K4me2 at enhancers and promoters. We identified a significant gain of H3K27ac in HACNS1 versus wild type limb buds at the promoter of the nearby gene Gbx2 (Fig.2, Supplementary Fig.2, Supplementary Data3). While significant H3K27ac enrichment was found in all three lines at the Gbx2 promoter compared to input controls, H3K27ac levels were not significantly increased at Gbx2 in limb buds with the chimpanzee ortholog compared to wild type, indicating the gain of activity is specific to HACNS1. H3K4me2 was also enriched at the promoter of Gbx2 in all three lines compared to input controls (Supplementary Fig.2B). After multiple testing correction, we did not identify any significant differentially marked regions outside of the HACNS1-Gbx2 locus between either the HACNS1 homozygous or the chimpanzee ortholog line compared to wild type for each chromatin mark in either tissue (Supplementary Data3). Moreover, we did not detect a significant increase in either H3K27ac or H3K4me2 levels at the Gbx2 promoter in the pharyngeal arch in HACNS1 homozygous embryos (Supplementary Figs.2C, D). This may be due to a lack of statistical power to detect small differences in histone modification levels given the number of replicates in the analysis, or our use of whole tissues to map histone modification profiles, which could obscure spatially restricted changes.

Gbx2 encodes a transcription factor with a highly complex expression pattern and multiple functions in the developing embryo. GBX2 has been implicated in midbrain and hindbrain development26,27, guidance of thalamocortical projections28,29, ear development30, and pharyngeal arch patterning31. Gbx2 is expressed in developing mouse limb at E10.5; however, its role in limb development remains undetermined as no limb phenotype has been reported in Gbx2 knockout mice27. HACNS1 and GBX2 are located in the same topologically associated domain (TAD), and TADs have been shown to restrict enhancer interactions to genes within their boundaries32,33. In E11.5 limbs from HACNS1 homozygous mice, the only significant increases in H3K27ac we detected in this TAD were at the HACNS1 knock-in locus and at the endogenous Gbx2 promoter (Supplementary Data3). Together, these results suggest that Gbx2 is a regulatory target of HACNS1, evoking the hypothesis that the gain of function in HACNS1 might alter Gbx2 expression in the limb.

To visualize potential expression changes resulting from HACNS1-driven upregulation of the Gbx2 promoter in HACNS1 homozygous mouse embryos, we used in situ hybridization (ISH) (Fig.3). We analyzed Gbx2 expression in >50 E11.5 embryos per genotype (Fig.3A, B; Methods and associated Source Data). In wild type embryos, we observed single foci of Gbx2 expression in forelimb and hindlimb (Fig.3A, left). In contrast, embryos homozygous for HACNS1 showed substantially increased Gbx2 expression in both forelimb and hindlimb (Fig.3A, right). Embryos homozygous for the chimpanzee ortholog showed a weak increase in Gbx2 expression compared to wild type (Fig.3A, B). Gbx2 expression in HACNS1 homozygous embryos was increased in two distinct anterior and posterior regions in the forelimb and hindlimb bud, as well as an anterior proximal region in the latter. Overall, Gbx2 expression in the limb bud was temporally dynamic in embryos of all genotypes. Embryos from the same litter vary in developmental age such that individual embryos collected at E11.5 range from E11 to E12. Therefore, we established a fine staging scheme to characterize changes in Gbx2 expression within this short developmental interval. We assigned embryos to 6 temporally ordered groups (designated T1-T6, and ranging from ~40 to 48 somites, although we did not use somite counts to stage embryos; see Methods for more details) according to crown-rump length and used a blinded approach to qualitatively assess staining patterns (Fig.3A, B and associated Source Data)10,34,35,36.

We identified differences in the distribution of Gbx2 expression in the forelimb and hindlimb buds of HACNS1 compared to both chimpanzee ortholog and wild type embryos across all 6 developmental time points (Fig.3B). At the earliest time point (T1), we found that Gbx2 was strongly expressed in distinct anterior-distal and posterior domains in HACNS1 forelimb and hindlimb buds (Fig.3A, B). Robust expression of Gbx2 in HACNS1 limb buds persisted through the remaining time points (up to T6), though the size of the anterior and posterior domains decreased over time. Strong expression of Gbx2 in HACNS1 homozygous embryos persisted for a longer period of time in hindlimb than in forelimb, consistent with the delayed developmental maturation of the former37. In addition, HACNS1 homozygous embryos showed a hindlimb-specific anterior-proximal expression domain adjacent to the body wall across all 6 time points (Fig.3A, B, Supplementary Fig. 3A).

In contrast to the robust expression observed in limb buds of HACNS1 homozygous mice, Gbx2 expression in limb buds from both the chimpanzee ortholog and wild type lines was weak and mostly evident at early time points (Fig.3). Chimpanzee ortholog line embryos and wild type embryos both showed weak distal Gbx2 expression foci in early forelimb and hindlimb that were generally restricted to the anterior limb bud (Fig.3A, B). Weak distal expression was primarily restricted to approximately T1-T2 in wild type forelimb but persisted until approximately T4 in a subset of embryos with the chimpanzee ortholog (Fig.3A, B). Weak distal expression persisted in hindlimb through T5-T6 in both the chimpanzee line and wild type (Fig.3B, bottom). These findings suggest that the chimpanzee ortholog line exhibits a modest increase in Gbx2 expression compared to wild type, potentially due to primate-rodent sequence differences affecting enhancer activity that our experimental design was intended to control for (see Supplemental Note). However, the HACNS1 knock-in line exhibits profound changes in Gbx2 limb bud expression compared to both. Together, these findings suggest that HACNS1 drives spatial and quantitative changes in Gbx2 expression in the limb, as well as a temporal extension of expression compared to wild type.

In addition to the forelimb and hindlimb bud, Gbx2 was also expressed in the neural tube, diencephalon, and pharyngeal arches of embryos homozygous for HACNS1 or the chimpanzee ortholog, and in wild type embryos (Fig.3A)27,29,31. Whereas Gbx2 expression was primarily restricted to the first pharyngeal arch in embryos with the chimpanzee ortholog and in wild type embryos, it expanded dorsally into the second pharyngeal arch in HACNS1 homozygous embryos during T1-T5 (Supplementary Fig.3B). However, we chose to focus on limb, due to the absence of a significant gain in H3K27ac marking at the HACNS1 knock-in locus or the Gbx2 promoter in HACNS1 homozygous versus wild type pharyngeal arch (Supplementary Fig.2).

In order to quantify the gain of Gbx2 expression in HACNS1 knock-in mice, we used real-time quantitative reverse transcription PCR (RT-qPCR) in pooled forelimb and hindlimb buds from embryos homozygous for HACNS1, the chimpanzee ortholog, and the endogenous mouse locus at time points T1-T6. We found that Gbx2 expression was increased in forelimb and hindlimb of HACNS1 embryos versus both chimpanzee ortholog and wild type at all 6 time points (Supplementary Fig.3C and associated Source Data). Although we detected an increase in Gbx2 expression in forelimb and hindlimb of embryos with the chimpanzee ortholog versus wild type at early time points, this change was substantially weaker than that between HACNS1 and wild type or HACNS1 and the chimpanzee line. Consistent with our ISH results, Gbx2 expression in HACNS1 knock-in line forelimb and hindlimb was strongest at the earliest time points and persisted longer in hindlimb than in forelimb (Fig. S3C). While Gbx2 expression declined over time in all three genotypes, it persisted longer in HACNS1 knock-in line forelimb and hindlimb.

In order to identify the specific cell types expressing Gbx2 as well as genes that are co- expressed with Gbx2 in the developing limb of all three genotypes, we performed single-cell RNA-sequencing (scRNA-seq) in E11.5 hindlimbs, which showed the most pronounced upregulation of Gbx2 in spatial and quantitative expression analyses (Fig.3, Supplementary Fig.3). Using the 10x Genomics scRNA-seq platform for cell barcoding, library preparation, and sequencing, we obtained transcriptomes from ~10,000 cells per genotype. We used the Seurat toolkit for data preprocessing and library size normalization (Methods)38. During pre-processing, we removed endothelial and blood cells (Cd34-positive; Pf4-positive; Hbb-positive), as our analysis indicated that these developmentally and transcriptionally distinct cell types do not express Gbx2 in any of our datasets39,40. After normalization of the filtered data using the SCTransform method in Seurat and integration of data from all samples into a single dataset using the Seurat v3 integration workflow (Methods), we performed clustering analysis on the integrated dataset to identify cell type categories present in all three genotypes38. To visualize similarities between cells, we used Uniform Manifold Approximation and Projection (UMAP), a dimensionality reduction method for data visualization, followed by the Louvain method for community detection to identify cell subtypes41,42.

This analysis revealed three distinct groups: (a) mesenchymal cell subtypes based on expression of the known markers Sox9 (clusters 1, 3, 4), Bmp4 (cluster 2), Shox2 (cluster 1), and Hoxd13 (clusters 14) (Fig.4A); (b) non-mesenchymal cell types, including myogenic cells (cluster 5, Myod); and (c) ectodermal cells (cluster 6, Fgf8) (Fig.4A)43,44,45,46,47,48. Furthermore, our analysis revealed finer separation of mesenchymal cells according to known limb patterning markers. We first examined the expression of known proximal-distal limb bud markers Meis1, Hoxa11, and Hoxd13 (proximal, medial, and distal, respectively)49. Cells expressing each of these markers showed a distinct localization in the UMAP embedding (Meis1+cells in the top left, Hoxa11+cell in the center, and Hoxd13+cells in the lower right; Fig.4B, upper left), suggesting our analysis recovered transcriptional and cell-type transitions along the proximal-distal patterning axis (Fig.4B).

A Left: UMAP embedding of HACNS1 homozygous, chimpanzee ortholog line, and wild type cells. The colors indicate cell clusters identified by Louvain clustering. Right: Expression of known limb bud cell-type marker genes in each cluster. Black dots denote cluster mean expression. B UMAP embedding of hindlimb bud cells from HACNS1 homozygous, chimpanzee ortholog line, and wild type, showing expression of proximal-distal, anterior-posterior, chondrogenesis-apoptosis, and non-mesenchymal markers. See text and Supplementary Fig.4 for details. C Expression of Gbx2 in each Louvain cluster, separated by genotype. Dots denote cluster mean expression. D UMAP embeddings illustrating cells expressing Gbx2 (indicated in red) in HACNS1 homozygous, chimpanzee ortholog line, and wild type cells. All gene expression data shown in plots and UMAP embeddings (AD) were imputed using ALRA and centered and scaled using z-scores (see Methods)81.

We also found that the first axis of the UMAP embedding clearly recapitulated known gene expression gradients along the anterior-posterior limb bud axis based on expression of the anterior-proximal marker Irx3, the anterior marker Zic3, and the posterior-proximal marker Shh (Fig.4B, upper right)43,50,51,52. Using markers of chondrogenic (Sox9, Shox2) versus non-chondrogenic (Bmp4) mesenchyme, we found that the second UMAP axis followed the chondrogenic versus interdigital apoptotic fate gradient (Fig.4B)43,45,48. We also found that the expression patterns of these markers were broadly conserved between genotypes, with each genotype showing comparable subsets of proximal, distal, anterior, posterior, chondrogenic, and non-chondrogenic cell types (Supplementary Figs.4A, B). Collectively, our scRNA-seq analyses identified specific conserved cell types and spatial transcriptional gradients in the developing hindlimb bud across all three genotypes.

We then sought to define genotype-specific differences in Gbx2 expression. In order to identify the cell types expressing Gbx2 in HACNS1 homozygous hindlimb buds, we examined the distribution of Gbx2-positive cells across cell clusters in all three genotypes. We found that Gbx2 was upregulated in HACNS1 homozygous hindlimbs versus chimpanzee ortholog and wild type hindlimbs, primarily in the mesenchymal cell clusters (clusters 14), consistent with the ISH and RT-qPCR expression analyses (Fig.3, Supplementary Fig. 3Cand Fig.4C). In HACNS1 homozygous hindlimbs, 24% of cells expressed Gbx2, 96% of which were mesenchymal cells, whereas less than 1% of cells in chimpanzee ortholog and wild type hindlimbs expressed Gbx2 (Methods; Supplementary Data4). Greater than 98% of Gbx2-positive cells in the chimpanzee ortholog line and wild type hindlimb were mesenchymal, and the majority of Gbx2-positive cells in each of these lines were assigned to Louvain cluster 2 (70% and 68%, respectively; Supplementary Data4). Only one non-mesenchymal cell from chimpanzee ortholog hindlimb and one non-mesenchymal cell from wild type hindlimb was Gbx2-positive (Fig.4C; Supplementary Data4). UMAP embedding of cells revealed that Gbx2-positive cells in HACNS1 homozygous hindlimb buds largely clustered within a distinct subset of mesenchymal cells belonging primarily to Louvain clusters 1, 2 and 4 (Fig.4C, D; Supplementary Data4). The genotype-specific differences in Gbx2 expression were consistent between the imputed and unimputed data as well as across individual replicates (Supplementary Fig.4C, D).

To identify genes whose expression is associated with Gbx2, we used k-Nearest-Neighbors Conditional-Density Resampled Estimate of Mutual Information (kNN-DREMI), which computes scores quantifying the strength of the relationship between two genes53,54. Using kNN-DREMI scores, we ranked each gene expressed in HACNS1 homozygous limbs by the strength of its association with Gbx2. To determine if genes associated with Gbx2 were enriched in particular functions, we then performed Gene Set Enrichment Analysis (GSEA) on this set of ranked genes. We found that Gbx2 expression was associated with genes in several limb development-related ontologies, including Cartilage Morphogenesis (KolmogorovSmirnov (KS) P=1.61103) and Regulation of Chondrocyte Differentiation (KS P=5.84103); the latter overlapped considerably with Collagen Fibril Organization (KS P=1.30104) (Fig.5A, Supplementary Data5). These results indicate that in the HACNS1 homozygous hindlimb, Gbx2 is co-regulated with genes expressed in condensing mesenchymal cells destined to become chondrocytes (Fig.5A).

A Ontology enrichments of genes with expression associated with Gbx2 expression (top) and the relative likelihood of the HACNS1 knock-in condition (HACNS1 RL, bottom) in HACNS1 homozygous mesenchymal cells. The log-transformed Gene Set Enrichment Analysis KolmogorovSmirnov P value for each category is plotted on the x-axis. Ontologies shown are those overlapping in the Gbx2 expression and HACNS1 RL ontology enrichment lists. See also Supplementary Data5 and 6. B HACNS1 RL and Gbx2 kNN-DREMI scores are plotted for all genes. Genes ranked in the top 20% of kNN-DREMI scores in the Chondrocyte Differentiation ontology (GO:0002062) for the union of the HACNS1 RL and Gbx2 kNN-DREMI analysis gene lists are colored in red and labeled. Dotted lines indicate the top 20% of values for each dataset. C Heatmap showing expression of genes belonging to the ontology Chondrocyte Differentiation (GO:0002062) in all HACNS1 homozygous mesenchymal cells (Louvain clusters 14). Hierarchical clustering was used to determine the order of cells (in columns) and genes (in rows). The bar at the top of the heatmap shows Gbx2-positive and Gbx2-negative cells in red and gray, respectively. D Expression of selected genes in Gbx2-positive (red) versus Gbx2-negative (gray) mesenchymal cells belonging to Louvain clusters 1 and 2, partitioned by proximal-distal axis markers as follows: Proximal cells (Prox) are Meis1+, Hoxd13, Hoxa11; distal cells (Dist) are Hoxd13+, Hoxa11 and Meis1; and intermediate cells (Mid) are all remaining Hoxa11+cells. Cells were randomly down-sampled to enable comparison of equal numbers of Gbx2-positive and Gbx2-negative cells. Larger red and gray dots respectively denote mean expression of each indicated gene in each group in Gbx2-positive and Gbx2-negative cells. All gene expression values shown in C and Dwere imputed using ALRA and centered and scaled using z-scores (see Methods)81.

We also used Manifold Enhancement of Latent Dimensions (MELD) to quantify the differences in the transcriptional profiles of HACNS1 homozygous limb buds compared to the chimpanzee ortholog line and wild type. MELD is an unsupervised learning algorithm and is therefore an orthogonal approach that is nave to our identification of Gbx2 as the target of HACNS155. MELD uses graph signal processing to quantify the relative likelihood of observing each cell in each of multiple experimental conditions based on its transcriptional profile. In this case, MELD is used to quantify the relative likelihood (RL) of observing a cell in the HACNS1 hindlimb versus the chimpanzee ortholog or wild type hindlimb. Rather than explicitly classifying genes as differentially expressed in one condition versus another, the RL value can be used to identify trends in gene expression across cells that are associated with the HACNS1 knock-in condition55.

To identify overall gene expression patterns characteristic of HACNS1 homozygous hindlimb bud cells, we used kNN-DREMI to associate gene expression with the HACNS1 knock-in condition (HACNS1 RL) calculated by MELD. We then used the resulting gene rankings to identify enriched biological functions via GSEA, as described above for Gbx2 expression. We found that genes associated with both HACNS1 RL and Gbx2 expression converged on related biological processes. Performing GSEA using genes ranked by mutual information with HACNS1 RL revealed significant enrichment of the Chondrocyte Differentiation ontology (KS P=7.00105), along with four other categories also significantly enriched in the Gbx2 expression analysis: Hindlimb Morphogenesis (KS P=1.42103), Cartilage Development Involved in Endochondral Bone Morphogenesis (KS P=1.11103), Collagen Fibril Organization (KS P=1.40105), and Positive Regulation of Phosphatidylinositol 3-kinase Signaling (KS P=4.40105), of which the last two are also implicated in chondrocyte differentiation (Fig.5A, Supplementary Data6)56,57. Collagen Fibril Organization is the most significantly enriched GO term for genes associated with HACNS1 RL and is the second most enriched for genes associated with Gbx2 expression. The top GO term for genes associated with Gbx2 expression, Roof of Mouth Development, (KS P=1.10105), shares >20% of its genes with Embryonic Hindlimb Morphogenesis (KS P=1.18103)58. This illustrates that many genes involved in limb development are also implicated in craniofacial development, and likely accounts for why craniofacial development-related GO terms were also enriched in our analysis.

These findings led us to examine the expression patterns of chondrocyte differentiation-related genes in HACNS1 homozygous mesenchymal cells belonging to Louvain clusters 14 (Fig.5C). We clustered HACNS1 homozygous mesenchymal cells by humanized RL and Gbx2 expression and examined the expression of the Chondrocyte Differentiation ontology genes within Gbx2-positive cells (Fig.5C). This clustering analysis revealed higher expression of positive regulators of chondrocyte differentiation (e.g., Sox9, Col2a1, Bmp2, and Runx2) specifically in Gbx2-positive versus Gbx2-negative humanized mesenchymal cells, supporting that HACNS1-driven upregulation of Gbx2 occurs in chondrogenic cells (Fig.5C)57,59,60,61. We also identified a subset of Gbx2-positive humanized cells that were also positive for Bmp4, which is expressed in the apoptotic interdigital domains62. These findings suggest that upregulation of Gbx2 is associated with the interdependent pathways of digit condensation and interdigital cell fate specification required for digit morphogenesis.

The gain of a Gbx2 anterior proximal hindlimb expression domain in HACNS1 knock-in mice suggests that Gbx2 upregulation may impact development of multiple segments in the limb (Fig.3A, B). To explore this further, we characterized the distribution of Gbx2 expression along the proximal-distal (PD) axis in the HACNS1 homozygous hindlimb single-cell dataset. We defined three groups of mesenchymal cells along the (PD) axis based on expression of the patterning markers shown in Fig.4: distal (Hoxd13+, Hoxa11, Meis1); proximal (Hoxd13, Hoxa11, Meis1+); and intermediate (all remaining Hoxa11+) cells. Gbx2 was expressed in all three subtypes (Fig.5D, top). We also examined expression of several Chondrocyte Differentiation ontology genes (selected from Fig.5C) in Gbx2+and Gbx2- cells in each subtype. Cells expressing Gbx2 also expressed chondrocyte differentiation markers along the PD axis, with a subset of markers showing modestly increased expression in Gbx2+cells (Fig.5D, middle and bottom). Together, these results link Gbx2 upregulation to chondrocyte differentiation in multiple developing regions of the hindlimb.

To determine if Gbx2 upregulation and downstream transcriptional changes in HACNS1 limb buds affect digit formation or overall limb morphology, we performed morphometric analysis of skeletal preparations for embryos homozygous for HACNS1, the chimpanzee ortholog, and wild type embryos (Methods). We performed morphometric analysis at E18.5 in order to capture any major phenotypic effects of targeted humanization that occurred by the end of embryonic skeletogenesis22. We did not detect gross morphological differences among genotypes; the three major limb segments (autopod, zeugopod, and stylopod) were present in both HACNS1 skeletons and chimpanzee ortholog skeletons (Supplementary Fig.6). We also examined digit length (normalized to body size based on the length of the ossified humerus), and intradigital (phalange to metacarpal or metatarsal length) and interdigital ratios. Again, we found no significant differences in digit length or autopod proportions between genotypes (Supplementary Fig.6, Supplementary Data79; Source Data). Although these analyses failed to detect morphological differences, it is possible that subtle phenotypes do indeed exist in the HACNS homozygous limb, which may be revealed in future studies.

See the original post here:

Modeling uniquely human gene regulatory function via targeted humanization of the mouse genome - Nature.com