The establishment of Hoxd13Q50R mutant mice by TALEN technique

In order to achieve efficient genome editing of the murine ortholog of the human HOXD13 gene, we designed a pair of TALEs targeting the second exon of the murine Hoxd13 gene. Our objective was to induce an A-to-G mutation at position 1 769 (50 position in Hoxd13 homeodomain). Additionally, to facilitate genotyping of the mice in subsequent experiments, we incorporated a single NdeI site into the donor DNA template by introducing a T-to-A mutation at the eighth nucleotide upstream of the target site (as depicted in Fig. 1a). It is worth noting that this introduced variant corresponds to a synonymous change, thereby encoding the same wild-type amino acid.

Fig. 1

The establishment of Hoxd13Q50R knock-in mice. a Schematic illustration of mutation knock-in using the TALEN system. b No.10, 18, 19, and 21 of first-generation mice were confirmed to carry a heterozygous mutation of Hoxd13Q50R by NdeI digestion. (M = DNA marker, N = naïve mice). c Sanger-sequencing result of wild-type (WT/WT), Hoxd13Q50R heterozygous mice (Q50R/WT), and Hoxd13Q50R homozygous mice (Q50R/Q50R)

In total, we obtained 21 F0 generation mice and subsequently identified four individuals (designated as F0-10, F0-18, F0-19, and F0-21) harboring the Hoxd13Q50R mutation via single restriction enzyme NdeI digestion (Fig. 1b). Subsequently, PCR-Sanger sequencing was used to identify potential off-target effects in exon2 of genes(Hoxa2, Hoxa10, Hoxa11, Hoxb9, Hoxb13, Hoxc8, Hoxc11, Hoxd9) exhibiting high homology with Hoxd13. Our results confirmed that the Hoxd13Q50R mutation was the sole mutation induced in the four identified F0 generation mutant mice. By mating Hoxd13Q50R heterozygous female and male mice, we obtained Hoxd13Q50R homozygous mice. Figure 1c shows the genotypes of the wild-type, heterozygote, and homozygote.

Hoxd13Q50R mutant mice display limb malformations

The limb of the mouse consists of five digits that follow an antero-posterior sequence (denoted as 1, 2, 3, 4, and 5). Of these, digits 2-5 exhibit a high degree of morphological similarity due to each possessing three phalanges.8 Notably, we have observed a mild manifestation of brachydactyly and syndactyly, either cutaneous or bony, between digits 2 and 4 in adult heterozygous mice with Hoxd13Q50R mutation (Fig. 2a1, a2). As shown in Table 1, we performed a statistical analysis of the proportion of different degrees of syndactyly in 32 heterozygous Hoxd13Q50R mice. All homozygous mice of Hoxd13Q50R display a more severe phenotype characterized by syndactyly with webbing and lack of distinct phalanges in the forelimbs, and central syndactyly and brachydactyly in the hindlimbs (Fig. 2a3, a4, and b1). Micro computed X-ray tomography(micro CT) scanning revealed the presence of six digits in the hindlimbs of the homozygous mice (Fig. 2b2).

Fig. 2

Morphological analyses of the mice with different genotypes. a Macroscopic results. 2a1 and 2a2 show cutaneous or bony syndactyly between digits 2–4 in Q50R/WT mice forelimb (F) and hindlimb (H), respectively. 2a3 and 2a4 show truncated forelimb and central syndactyly hindlimb in Q50R/Q50R mice (red boxes marked syndactyly in Q50R/WT). b Micro CT scanning results. 2b1 show syndactyly with webbing and no distinct phalanges in Q50R/Q50R mice forelimb. 2b2 show central syndactyly and brachydactyly in the hindlimb of Q50R/Q50R mice. c Skeletal preparations of WT/WT and Q50R/Q50R at postnatal day 4(P4). d Skeletal preparations of WT/WT, Q50R/WT, and Q50R/Q50R at postnatal day 0(P0) (red arrows marked different ossification levels of phalanges between WT/WT and Q50R/WT). 2d1, 2d2, 2d4, and 2d5 show the delayed ossification of phalanges in the mice of Q50R/WT compared with the wild-type. 2d3 and 2d6 show the dysplasia and disorganization of bone and cartilage in Q50R/Q50R’s autopods

Table 1 Distribution of syndactyly manifestation in 32 heterozygous mice with Hoxd13Q50R

The animals at postnatal day 4 (P4) were subjected to skeletal preparations and subsequently analyzed. No significant differences were found in the development and structure of long bones, and cranial bones between wild-type mice and homozygous mice. However, statistical analysis revealed that the vertebral column length of homozygous mice was 1–4 mm shorter than that of wild-type mice, with a significance level of P < 0.001. (Fig. 2c, Table S1 and Fig. S1). Furthermore, the comparison between Hoxd13Q50R heterozygotes and their wild-type counterparts showed a delayed ossification of phalanges in the former (Fig. 2d1, d2, d4, and d5). In contrast, the autopod of Hoxd13Q50R homozygotes exhibited dysplasia and disorganization of bone and cartilage (Fig. 2d3, d6).

Decreased BMP signaling regulates interdigit PCD, leading to syndactyly in Hoxd13Q50R

Firstly, we aimed to elucidate the pathogenic correlation between SDTY5 and the Hoxd13Q50R mutation during embryonic development. In a developing limb, critical regions of programmed cell death (PCD) are situated in the ectoderm of the AER and the undifferentiated mesoderm, concomitant with the establishment of prechondrogenic condensation of the skeleton.22 The bone morphogenetic protein (BMP) signaling pathway plays a positive role in regulating interdigital PCD.23

To explore the potential involvement of BMP signaling in syndactyly development in Hoxd13Q50R mutant mice, we conducted an analysis of Bmp2 mRNA expression. Specifically, we utilized whole-mount in situ hybridization (WISH) to visualize and quantify Bmp2 expression patterns across several developmental stages of the mutant mice. The results indicated a clear down regulation of Bmp2 expression in conjunction with ectopic expression in affected tissues (Fig. 3). At E11.0 (TS = 18), the LBs of wild-type exhibited posterior and distal expression of Bmp2, whereas the mutant mice showed only minor posterior expression. At E11.5 (TS = 19), Bmp2 expression was observed between digits in the forelimbs and hindlimbs of both wild-type and Q50R/WT mice, but was not observed in Q50R/Q50R mice. The expression level of Bmp2 in the distal regions of Q50R/WT and Q50R/Q50R mice was found to be gradiently decreased compared to the wild-type. At E12.0(TS = 20), Bmp2 was expressed in posterior and interdigital regions, with downregulation of Bmp2 expression in Hoxd13Q50R mutant mice, particularly in interdigital regions. In Q50R/Q50R mice, Bmp2 expression in interdigital regions was barely detectable. Bmp4 expression pattern was also examined at E11.5, but no transcriptional changes were observed (data not provided).

Fig. 3

Gene expression analysis in BMP signaling. a The expression levels of Bmp2 during various developmental stages of limb buds. The experimental groups include the mice of wild-type (WT), heterozygotes (Q50R/WT), and homozygotes (Q50R/Q50R). The forelimbs and hindlimbs are represented by “F” and “H,” respectively. “TS” denotes the time with somite numbers. The signals of Bmp2 expression were detected by WISH indicated by red arrows. b Western blot analysis using SMAD1/5 and Phospho-SMAD1/5 in forelimbs (F) and hindlimbs (H) of WT/WT, Q50R/WT and Q50R/Q50R mice at E11.5. c TUNEL staining in interdigital regions of forelimbs from WT/WT and Q50R/Q50R mice at E13.5(TS = 22). Apoptotic cells were red. DAPI was used as a counterstain

The critical event in BMP signal transduction involves the phosphorylation of transcription factors Smad1, 5, and 8 (referred to collectively as Smad1/5/8) through BMP receptor mediation, ultimately resulting in the transcription activation of BMP-induced gene in the nucleus.24 Total SMAD1/5 and phospho-SMAD1/5 were examined in hindlimbs of WT/WT, Q50R/WT and Q50R/Q50R mice at E11.5 by western blot and the ratio of pSMAD1/5:total SMAD1/5 in the limb buds of Q50R/Q50R was significant decreased than that in the wild-type (Fig. 3b). These data demonstrated that the decreased expression of Bmp2 in the limbs of the mutant mice regulates the subsequent BMP signaling pathway by influencing the phosphorylation levels of SMAD1/5.

Subsequently, we employed the TUNEL staining method to detect apoptotic cells in the interdigital regions of E13.5 wild-type and homozygous mice. The signal intensity of apoptotic cells in the interdigital regions of wild-type mice was significantly stronger than that in the homozygous mice (Fig. 3c). This further elucidates that the downregulation of BMP signaling indeed affects the PCD in the interdigital regions, thereby resulting in syndactyly.

Based on these findings, it can be inferred that the reduced expression of Bmp2 in the distal regions of Hoxd13Q50R mutant mice during the onset of limb patterning results in the untimely and irregular PCD of the interdigital tissue. This, in turn, interferes with the usual process of chondrogenesis and digit differentiation, ultimately resulting in anomalies such as syndactyly and limb malforamtions.

The expression of the SHH/GREM1/AER-FGF epithelial-mesenchymal feedback signaling in Hoxd13 mutant mice

As previously mentioned, the A-P axis is formed through the activity of ZPA signaling centers that are under the control of the morphogenetic hormone SHH. The regulation of distal LB growth occurs through the SHH/GREM1/AER-FGF epithelial-mesenchymal (e-m) feedback signaling system17 (Fig. 4a). Within this feedback signaling system, the 5′-HOXD transcription factor initially regulates SHH expression, which subsequently promotes the expression of GREM1. GREM1 then represses the expression of BMPs in the interphalangeal region, acting as a BMP repressor.22 In Q50R/Q50R mice, limb patterning along both the A-P and P-D axes fails completely. Therefore, experiments were conducted to identify the expression patterns of the genes involved in the SHH/GREM1/AER-FGF e-m feedback signaling system.

Fig. 4

Gene expression analyses of SHH/GREM1/AER-FGF e-m signaling in limb buds by WISH. a The illustration of apical ectodermal ridge (AER) and zone of polarizing activity (ZPA) with A-P axis, P-D axis and D-V axis and the SHH/GREM1/AER-FGF e-m signaling feedback loop. b Grem1 expression pattern in WT/WT, Q50R/WT, and Q50R/Q50R at E10.5 (green boxes marked Grem1 signal). c Shh expression pattern in WT/WT, Q50R/WT, and Q50R/Q50R at E10.5 (red arrows marked Shh signal). d Shh-Fgf8 expression pattern in WT/WT, Q50R/WT, and Q50R/Q50R at E11.5 (red arrows marked Shh signal; blue arrows indicated Fgf8 signal). e Fgf8 expression pattern in WT/WT, Q50R/WT, and Q50R/Q50R at E11.5 (black arrows marked Fgf8 signal; blue double arrows were used as scale bars to indicate that the hindlimbs of WT/WT, Q50R/WT, and Q50R/Q50R were of the same size)

Based on the previous analysis of Bmp2 expression, we performed WISH to detect the expression pattern of BMP antagonist Grem1 at the RNA level, and we did find a significant difference between Q50R/Q50R and wild-type mice at E10.5. The expression level of Grem1 in Q50R/Q50R mice forelimbs was higher than that in wild-type (Fig. 4b). However, no obvious expression changes was detected at later development phases.

To evaluate the expression pattern of the Grem1 upstream morphogen Shh and Bmp2 downstream Fgf8, we performed a combined WISH analysis. Our findings revealed that Shh was strictly expressed in the zone of polarizing activity (ZPA) in wild-type mice. However, in Hoxd13Q50R mutant mice, we observed a slight diffusion of Shh expression sites around the ZPA at E10.5 (Fig. 4c). At E11.5, Q50R/Q50R mice had a shorter but broader Shh expression than WT/WT and Q50R/WT mice (Fig. 4d). Furthermore, our results indicated that Fgf8 had a broader expression in the AER of Hoxd13Q50R mutant mice, which suggested that Fgf8 was upregulated in the mutant mice (Fig. 4e) when compared to wild-type mice.

The ZPA is essential for developing the asymmetry in digits across the A-P axis, and Shh acts as a mediator of A-P patterning.14,25Fgf8 is one of the AER-Fgfs that is necessary for normal outgrowth and patterning of the limb.13 Early activated HoxA;D genes promote the maintenance of AER-FGFs by providing proximal identity to cells where they are expressed, and early activated HoxA;D genes ensure proper AER-FGF-dependent expression of late/distal HoxA;D genes.26 Our analysis indicates that the abnormal expression of Shh and Fgf8, combined with the downregulation of Bmp2 in Hoxd13 mutant mice, contributes to the failure of A-P and P-D axis patterning and a reduction in interdigital PCD. Ultimately, this leads to the development of the syndactyly and brachydactyly phenotype in Q50R/WT mice and severe limb malformations in Q50R/Q50R mice.

RNA-sequencing analysis of LBs of mice at E11.5

Transcriptomic profiling was conducted on murine buds of forelimbs and hindlimbs from the mice of wild-type and Hoxd13Q50R homozygotes at E11.5. A total of three individual biological replicates were selected for both groups to ensure reproducibility and statistical significance.

The RNA-sequencing analysis revealed significant differences in the expression of differentially expressed genes (DEGs) between the wild-type and homozygotes of Hoxd13Q50R (Fig. 5a). Based on the volcano plot analysis, there were 115 DEGs in the forelimbs between mice of Q50R/Q50R and wild-type, including 42 upregulated and 73 downregulated genes, and there were 274 DEGs in the hindlimbs between the mice of Q50R/Q50R and wild-type, with 213 upregulated genes, and 61 downregulated genes. 31 common DEGs were found, by combining the DEGs of the forelimbs and hindlimbs. According to P-values (P < 0.05), the 31 overlapped DEGs between the forelimbs and hindlimbs comprised 1190002N15Rik, Pcsk6, Plxna4, Tspan18, Lhx9, Cdc42ep3, Olfm1, Hoxd12, Evx2, Draxin, Smarcd2, Hey1, Tfap2b, Polg, Ppfibp2, Crmp1, Pam, Utrn, Bambi, Tbx4, Slc1a3, Cited2, Lhx2, Plod2, Klhl29, Slit3, Fzd8, Lix1, Pdgfrb, Sall1, and Scmh1 (Fig. 5b).

Fig. 5

The results of RNA-sequencing analysis. a Clustering of DEGs in different groups. Red indicated highly expressed genes; blue indicated poorly expressed genes. b Volcano plot and venn diagram show the number of upregulated and downregulated genes in Q50R/Q50R_F vs WT/WT_F and Q50R/Q50R_H vs WT/WT_H, and their overlapping genes. DEGs were considered at Padj value < 0.05. c GO analysis of DEGs in the forelimbs of Q50R/Q50R and WT/WT mice. d GO analysis of DEGs in the hindlimbs of Q50R/Q50R and WT/WT mice

The GO functional enrichment analysis revealed that the DEGs in the forelimbs of Q50R/Q50R and WT/WT mice were primarily associated with various developmental processes such as epithelial cell morphogenesis, embryonic organ development, stem cell differentiation, regulation of organ morphogenesis, arterial morphogenesis, leg-tail morphogenesis, limb morphogenesis, mesenchymal development, arterial development, growth and development involved in morphogenesis, limb development, and vascular morphogenesis (Fig. 5c). The DEGs between Q50R/Q50R and WT/WT mice hindlimbs were mainly involved in processes such as limb morphogenesis, Wnt signaling pathway, regulation of cell morphogenesis, epithelial cell morphogenesis, mesenchymal cell development, transmembrane receptor protein tyrosine kinase signaling pathway, and embryonic limb morphogenesis (Fig. 5d).

Subsequently, a subset of common DEGs of forelimbs and hindlimbs between Q50R/Q50R and wild-type were validated using WISH and RT-qPCR. Specifically, 1190002N15Rik encodes Golgi protein GoPro49, which is expressed in embryonic mesenchymal tissues. GoPro49 is known as a marker of the dental follicle and may have a function in the secretory pathway.27 In this study, we utilized WISH to investigate the expression pattern of 1190002N15Rik and discovered that it exhibited initial expression in the distal 2–4 digit regions at E11.5, followed by strong expression between the metatarsals at E12.5. RT-qPCR analysis indicated a decreased expression level of 1190002N15Rik in mutant mice compared to WT/WT, which suggests that mutant mice may display an abnormal rate of limb patterning (Fig. 6c). Furthermore, our observations revealed the absence of 1190002N15Rik expression at E11.5, but we noted widespread ectopic expression in the distal regions of the limb at E12.5 in Q50R/Q50R mice (Fig. 6a). Microfibril-associated protein 3-like (MFAP3L), belonging to the MAGP family, is involved in regulating cell proliferation, migration, and invasion in cancer, as reported in previous studies.28 Our findings indicated a strong signal of Mfap3l in interdigit regions of the limbs in WT/WT, moderate expression in Q50R/WT, and a slight distal expression in Q50R/Q50R at E11.5 (Fig. 6b). Meanwhile, RT-qPCR results indicated that the expression levels of Mfap3l in the LBs of wild-type at E11.5 were significantly higher compared to Hoxd13 mutant mice (Fig. 6c).

Fig. 6

The expression of 1190002N15Rik and Mfap3l in the LBs of mice with different genotypes. a The WISH results of 1190002N15Rik expression pattern in WT/WT, Q50R/WT, and Q50R/Q50R hindlimbs at E11.5 and E12.5. b The WISH results of Mfap3l expression pattern in WT/WT, Q50R/WT, and Q50R/Q50R hindlimbs at E11.5. c RT-qPCR detection of the expression of 1190002N15Rik and Mfap3l in forelimbs and hindlimbs of WT/WT, Q50R/WT, and Q50R/Q50R mice at E11.5. At least three replicates were performed. d RT-qPCR detection of the expression of Mfap3l and 1190002N15Rik after overexpressing Bmp2 in primary cultured chondrocytes from Q50R/Q50R mice (*P < 0.05)

To evaluate the correlation between 1190002N15Rik, Mfap3l, and Bmp2, we introduced a Bmp2 overexpression vector into primary cultured chondrocytes derived from homozygous mutant mice. RT-qPCR analysis revealed a significant upregulation of Mfap3l expression and a concurrent downregulation of 1190002N15Rik upon Bmp2 overexpression (Fig. 6d). Conversely, overexpressing MFAP3L and 1190002N15Rik in the human chondrocyte cell line C28/I2 did not induce significant alterations in BMP2 expression compared to the control group (Fig. S2). These findings suggest that Bmp2 may function as an upstream regulator of Mfap3l and 1190002N15Rik during limb development. However, the precise regulatory associations between Mfap3l and 1190002N15Rik and Bmp2, as well as Hoxd13, warrant further investigation and validation.