New insight into the development of synpolydactyly caused by expansion of HOXD13 polyalanine based on weighted gene co-expression network analysis

Abstract

Background

Synpolydactyly (SPD) is mainly caused by mutations of polyalanine expansion (PAE) in the transcription factor gene HOXD13 and the involved cell types and signal pathway are still not clear possible pathways and single-cell expression characteristics of limb bud in HOXD13 PAE mice was analyzed in this study.

Method

We investigated a previous study of a mouse model with SPD and conducted weighted gene co-expression network analysis (WGCNA) using a single-cell RNA sequencing dataset from limb bud cells of SPD mouse model of HOXD13 + 7A heterozygote.

Results

Analysis of WGCNA revealed that synpolydactyly-associated Hoxd13 PAEs alter the immune response and osteoclast differentiation, and enhance DNA replication. Bmp4, Hand2, Hoxd12, Lnp, Prrx1, Gmnn, and Cdc6 were found to play potentially key roles in synpolydactyly.

Conclusions

These findings evaluated the main genes related to SPD with PAE mutations in HOXD13 and advance our understanding of human limb development.

Backgound

Synpolydactyly (SPD; MIM 186000), also known as syndactyly type II, is a distal limb deformity characterized by fusion of the third and fourth fingers, the fourth and fifth toes, and the presence of extra fingers in the syndactylous web [1]. Synpolydactyly is caused by mutations in transcription factor gene HOXD13. Such mutations include polyalanine expansion (PAE), point mutations and a 14 bp small deletions, and majority of SPD patients result from an incomplete trinucleotide repeat expansion within the polyalanine coding region in the first exon of the HOXD13 gene. However, the involved cell types and signaling pathways remain unclear. Correlation between the severity of SPD and the PAE length of HOXD13 has been discussed for many years [2]. The expansion of two alanines does not have pathological consequences [3], and, to our knowledge, PAE of less than seven alanines has not been reported in SPD patients. This means that synpolydactyly can be caused only when polyalanine expansion reaches a certain threshold, such as a PAE of seven alanines or more [4].

A recent study suggested that a novel missense mutation of Hoxd13 (NM_000523: exon2: c.G917T: p.R306L) promoted osteoclast differentiation by regulating the Smad5/p65/c-Fos/Rank axis, which might provide a new insight into SPD development [5]. Few prior studies focused on osteoclasts. A triple of HOXD10, HOXD11, and HOXD13 knockdown significantly decreased the number of TRAP⁺-stained osteoclasts within tumor tissue while their number in bones was not affected [6]. The effect of SPD related HOXD13 mutation on osteoclasts may be a parallel regulatory process of multi-gene changes.

Although we found that PAE of HOXD13 affected limb development, it was unclear which cells co-expressed common transcription factors or coordinated regulation by independent transcription factors with similar activation patterns during this process. Thus, in this study we performed bioinformatics analysis using a single-cell RNA sequencing (scRNA-seq) dataset of limb bud cells from a previous study in order to analyze PAE-associated cell type–specific gene expression changes in wild-type (WT) and + 7A repeat-expanded HOXD13 allele (spdh) mice [7].

Methods

Database retrieval

Based on a retrieval of polyalanine-expanded HOXD13 allele related databases, a single-cell RNA sequencing (scRNA-seq) dataset of limb bud cells from a previous study were obtained from GSE128818. Advanced analysis was performed in these data.

Single-cell analysis

Single-cell analysis was used to assess and validate the expression level of hub genes. Quality control, dimensional reduction, and clustering of the scRNA-seq data of the mouse limb (GSE128818) were performed using Seurat (v3.1.0). After obtaining a gene expression matrix generated by Cell Ranger software [8], based on the filtered cell–gene expression matrix, we used Seurat to standardize the data and perform cell typing through reduction and clustering and to then identify characteristic genes of the different cell populations. Genes with equal or above two-fold change in relative expression levels between the mutant group and wildtype group, along with an FDR-adjusted p value < 0.05 were considered differentially expressed. Further details are listed in Supporting information Data S2.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) is a systematic biological method that describes the pattern of gene associations between different samples. It can be used to identify highly synergistically altered gene sets as well as candidate biomarker genes or therapeutic targets based on the associations of gene sets, and associations between gene sets and phenotype. WGCNA was performed using a WGCNA package in R [9] to classify genes into modules. The analysis prevealed modules, or groups of related genes, which can reflect either shared regulation by common transcription factors or coordinated regulation by independent transcription factors with similar patterns of activation. Specifically, these groups were determined based on the scRNA-seq data of E12.5 limb buds from wild type (WT) and spdh mice. The selection of a β value for establishing an adjacent matrix makes our gene distribution conform to a scale-free network according to the degree of connectivity: the soft threshold (power) represents the weight, and the ordinate represents the correlation between the connectivity k and p (k). It is generally required that the correlation between k and p (k) reaches 0.85 as the value of β. Modules with a correlation coefficient greater than 0.8 (i.e., the coefficient of dissimilarity is less than 0.2) are combined [9]. Gene expression modules and different cell-type relationships are presented between spdh and WT small conditional (sc)RNA-sequence (seq) data. The eigenvector gene of each module is calculated as the first principal component gene E of a particular module and represents the overall level of gene expression within the module. We then calculated Student asymptotic P values for these correlations and corrected for the effects of multiple tests using the false discovery rate (FDR) [10]. Modules of highly coexpressed genes with |correlation coefficients|> 0.6 and an FDR < 0.05 were considered cell-type–related modules, and genes in these modules were considered cell-type–associated genes. Finally, the module eigengene was calculated and the Pearson correlation analysis was performed; hub modules and cell cluster were thus identified.

To evaluate the function of putative genes, gene ontology (GO) enrichment analysis was performed using the R package ClusterProfiler [11]. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed using g: Profiler [12]. To determine the biological processes and involved signaling pathways. FDR < 0.05 was set up to identify the statistically significant enrichments in the gene clusters of the modules. Significance was assessed using Fisher’s exact test.

Statistical analysis

The R software (version: 4.0.5) was used to perform all statistical analyses and graphical representations. For continuous variables, Student’s t-tests were used to compare two distinct groups. For categorical variables, chi-square tests were performed.

Results

Gene modules derived from WGCNA based on gene expression of a scRNA-seq dataset of limb bud cells

To define the cell types, we annotated each cell type based on the expression levels of canonical cell-type–specific markers according to a previous study [7]. In total, we identified 16 cell states (erythrocytes, perichondrium, proximal chondrocytes, proliferating cells, interdigital mesenchymal cells, myoblasts, hematopoietic cells) using Seurat analysis. The expression levels of specific marker genes for each cell type are shown in dot plots (Fig. 1A). The top four differentially expressed genes (DEGs; sorted by adjusted p value, with an adjusted p value of < 0.05) in each cell type are shown (Fig. 1B).

Fig. 1

The expression levels of specific marker genes were labeled for each cell type using Seurat analysis. A Dot plots showing the expression levels of specific marker genes for each cell type using Seurat analysis. B Heatmap showing the top four differentially expressed genes (DEGs) (sorted by adjusted p value, with an adjusted p value of < 0.05) (column) in each cell type (row). The color key from purple to yellow indicates gene expression levels from low to high

Based on the characteristics of different genes in the limb bud, we used WGCNA to construct a co-expression network to identify modules and genes that were significantly related to cell-type–specific gene expression and HOXD13 PAE. The scale-free topologies of the networks were assessed for various β shrinkage parameter values according to the WGCNA user manual, wherein a β = 5 value provided a satisfactory fit to scale-free topology. Both the adjacent matrix and topological matrix were obtained according to the β value and whether the storage network under the selected β value approximated the scale-free topology (Fig. 2A). A topologically overlapping heat map of 1227 highly variable gene pairings was selected from scRNA-seq data. The genes were clustered using dissimilarity between them, and then the tree was cut into different modules using dynamic shearing: seven modules for spdh and WT limbs (Fig. 2B and C). Both rows and columns represent a single gene, with dark, turquoise, yellow, and blue modules indicating a high degree of topological overlap (Fig. 2D). Myoblasts had the most significant correlation with the turquoise module (correlation coefficient = 0.83, Fig. 2E). Interestingly, we found that the proportion of differentially expressed genes in both significant (correlation coefficients > 0.6) and other modules was not significantly different (Fig. 2F), suggesting that the different single-cell types do not exhibit expression changes. The proportion of differentially expressed genes in the blue module, which included Hoxd13, was also not significantly different from that of other modules (Fig. 2G), suggesting that genes in the blue module may not be involved in SPD.

Fig. 2

Weighted correlation network analysis of the scRNA-seq data. A Plot of scale-free topology and mean connectivity with regards to soft-thresholding power for samples from scRNA-seq data. Red line indicates an R² cut-off of 0.87. B Visualization of the eigengene network representing the relationships among the modules and the topology of genes. C Clustering dendrogram showing correlations between modules from scRNA-seq data. D A gene co-expression network of single cells from the limb bud was constructed by weighted correlation network analysis. The hierarchical clustering and module assignment of genes are shown along the left side and top. E Gene expression modules showing correlations in relation to cell type. Positive correlations are red while negative correlations are green. The module and the sample traits are related to the heat map, the row represents the module, and the column represents the trait. The values in the box indicate correlation and FDR. F Relative numbers of DEGs identified in each module. Black bars represent |correlation coefficients|> 0.6 and FDR < 0.05 modules, and gray bars represent all other modules. G The box plot shows the relative number of DEGs in significant (sig.) modules and not significant (not sig.) modules compared with the blue module (including Hoxd13)

Synpolydactyly-associated polyalanine expansions alter the immune response and differentiation

The gene significance across modules for myoblasts from spdh and WT scRNA-seq data was most significant for the turquoise module (Fig. 3A). The correlation coefficient between the turquoise module and the gene significance for myoblasts from spdh and WT scRNA-seq data was 0.97 (p < 1e-200, Fig. 3B). The top four differentially expressed genes of myoblasts in each module in spdh and WT scRNA-seq data were Tmsb4x, Fcer1g, C1qb, and Apoe (Fig. 3C).

Fig. 3

Synpolydactyly-associated polyalanine expansions alter the immune response and osteoclast differentiation. A Gene significance across modules for myoblasts from spdh and WT scRNA-seq data. B Module membership in the turquoise module and gene significance for myoblasts from spdh and WT scRNA-seq data. C Expression of four differential genes of myoblasts in each module from spdh and WT scRNA-seq data. D Treemap of the top 10 condensed gene ontology (GO) terms for the turquoise module from spdh and WT scRNA-seq data, with box size corresponding to the number of significant terms associated with the GO category. E Representative KEGG pathways of the turquoise module from spdh and WT scRNA-seq data. ‘‘Count’’ stands for gene number. The color keys from yellow to blue indicate the range of the p value

Examination of common transcription factors or coordinated regulation by independent transcription factors with similar patterns of activation within the turquoise module revealed immune response changes in spdh and WT limbs, such as the immune effector process, inflammatory response, and innate immune response (Fig. 3D, Supplementary Data S1).

Moreover, to further assess turquoise module–related signaling pathways, we performed KEGG enrichment analysis. As shown in Fig. 3E, the turquoise module was mainly enriched in cytokine–cytokine receptor interaction, osteoclast differentiation, leishmaniasis, Staphylococcus aureus infection, tuberculosis, natural killer cell–mediated cytotoxicity, chemokine signaling pathway, phagosome, Fc gamma R–mediated phagocytosis, and lysosome (Supplementary Data S1).

Synpolydactyly-associated polyalanine expansions alter DNA replication

The top four differential genes of proximal chondrocytes in each module in spdh and WT scRNA-seq data were Hoxd13, Hoxd12, Aldh1a2, and Creb5 (Fig. 4A). The gene significance across modules for proximal chondrocytes from spdh and WT scRNA-seq data showed that the blue module, which included Hoxd13, was the most significant module (Fig. 4B). The correlation coefficient between the blue module and the gene significance for proximal chondrocytes from spdh and WT scRNA-seq data was − 0.049 (p = 0.43, Fig. 4C).

Fig. 4

Synpolydactyly-associated polyalanine expansions alter the DNA replication. Comparing the scRNA-Seq data from the limb buds of spdh and WT indicated: A Expression of four differential genes of proximal chondrocytes in each cell type. B Module membership in the blue module and gene significance for proximal chondrocytes. C Gene significance across modules for proximal chondrocytes related to Hoxd13. D Protein–protein interaction network of the blue module and the schematics how might Hoxd13 and its partners might be involved in SPD phenotype. E Treemap of the top 10 condensed gene ontology (GO) terms, with box size corresponding to the number of significant terms associated with the GO category

We used STRING [13] to obtain protein–protein interactions. Bmp4, Hand2, Hoxd12, Lnp, Prrx1, Gmnn, and Cdc6 were identified as hub genes (Fig. 4D). Examination of common transcription factors or coordinated regulation by independent transcription factors with similar patterns of activation within the blue module revealed DNA replication changes between the limb buds of spdh and WT, such as regulation of DNA metabolic process, DNA repair, DNA conformation change, and DNA-dependent DNA replication (Fig. 4E).

Discussion

Synpolydactyly (SPD) is an autosomal dominant disease caused by mutations in HOXD13 gene. Although the genetic etiology of SPD is well known, but the detail functional mechanism of HOXD13 mutations remains unknown. Point mutations, frameshift mutations, or PAE mutations in HOXD13 can lead to the typical SPD phenotype [14]. Some findings highlight the value of long-read whole genome sequencing in elucidating the molecular etiology of congenital limb malformation disorders [15]. Other mutations, unlike polyalanine extension, which tends to form α-helix and causes protein aggregation in the cytoplasm as shown by molecular simulation and immunofluorescence, the c. 925A > T mutation of HOXD13 impairs downstream transcription of EPHA7 [16]. There is a report speculate that the 11,451 bp microdeletion at chr2:176,933,872–176,945,322 (GRCh37), which is located upstream of HOXD13 gene (associating with SPD1), which was identified as the smallest deletion upstream of the HOXD13 gene and not altering the sequence of the HOXD13 gene [17]. We identified a new manifestation of preaxial polydactyly in both hands in a pediatric patient with an expansion of seven alanines, a phenotype not previously noted in SPD patients [18]. Other results showed a 24bp duplication (c.183_206dupAGCGGCGGCTGCGGCGGCGGCGGC) in the exon 1 of HOXD13 in heterozygous form which was predicted to result in eight extra alanine (A) residues in N-terminal domain of HOXD13 protein [19].

HOXD13 is the 5’ terminal member of the HOXD cluster and has two coding exons. Exon 1 contains an incomplete trinucleotide repeats of GCN (which encodes polyalanine, with N representing A, C, G, or T) and encodes a highly conserved homeobox domain. HOXD13 protein is a transcription factor that must enter the nucleus to regulate downstream target genes [20]. Therefore, the HOXD13 protein may have a nuclear localization signal [21]. Additional amplification of more than seven alanines disrupts this functional region, causing the protein aggregation to form inclusion body in the cytoplasm, thereby preventing the HOXD13 protein from entering the nucleus [22]. These data facilitate the interpretation of synpolydactyly radiographs by future automated tools [23]. In addition, a deep understanding of the key roles of HOXD13 in the most cell types in the bone marrow microenvironment need to explore.

Here, we provided a deep understanding of HOXD13 PAE based on scRNA-sequencing of limb bud samples from synpolydactyly and WT mice. Based on this dataset, we identified osteoclast differentiation related to different cell-type–specific co-expression genes. In line with our results, previous results demonstrated that posterior HOXD gene expression similarly stimulated RUNX2 expression and enhanced the expression of osteoblast-specific osteocalcin (BGLAP), pre-/osteoclast-specific PDGFB, and PTHLH [6]. HOXD13 PAE may alter osteoclast-specific gene expression. Our study is consistent with previous study about the association between osteoclast and SPD in the limb bud of mice [5].

Gene ontology terms of the turquoise module revealed an immune response difference between spdh and WT limbs, which may be supported by a previous study in which knockdown of HOXD11 or HOX13 significantly suppressed lung metastasis in immune-deficient mice [6], which is consistent with our results. In addition, our results indicate that Hoxd13 PAE altered the expression of Bmp4, Hand2, Hoxd12, Lnp, Prrx1, Gmnn, and Cdc6. This effect is consistent with recent reports that Gmnn is required for Hox gene to regulate limb pattern [24]. Many HOX proteins interact with the DNA replication licensing regulator geminin and bind a characterized origin of replication [25]. DNA replication and immune systems are the primary defenses against endogenous and exogenous threats, and defects in these systems are at the core of many health problems, including infections, autoimmune/inflammatory diseases, cancer, and other age-associated diseases. Replication stress can elicit an innate immune response by releasing DNA particles into the cytoplasm. HOXD13 interacts with the CDC6 loading factor, promotes pre-replication complex protein assembly at origins, and stimulates DNA synthesis in an in vivo replication assay [26]. Lnp (Lunapark) was strongly expressed in both limbs and genital buds, and co-expressed with Hoxd genes, especially with Hoxd13 [27]. Prrx1 and Hoxd13 were identified as the partner genes, which play important roles in the development of limb [28]. HOXD13 represses the transcription of SMAD1 and blocks BMP4-SMAD1-induced epithelial mesenchymal transition [29]. In the mouse model of synpolydactyly, polyalanine expansions in Hoxd13 alter the transcriptional program. Our results combine to other's finding, prompt that HOXD13 involved in the DNA synthesis and regulation of hub genes. PAE is believed to modify DNA replication, which subsequently results in DNA degradation and boost the immune system.

There are two limitations in the present study. First, the generalizability of findings to human synpolydactyly and the specific conditions under which the study was conducted. Second, we did not provide experimental verifying, and further validation of the related mechanism should be done through bench experiments on cell and animal. In this study, we aimed to explore the molecular mechanism using weighted gene co-expression network analysis (WGCNA) of a single-cell RNA sequencing dataset from synpolydactyly limb bud cells. Thus, we mainly focus on the results of bioinformatics analysis.

Conclusions

In summary, WGCNA analysis revealed that synpolydactyly-associated Hoxd13 PAEs alter the immune response and enhance DNA replication. Bmp4, Hand2, Hoxd12, Lnp, Prrx1, Gmnn, and Cdc6 were found to play key roles in the development of SPD. Such findings identify possible signaling pathways involved in the pathogenesis of spdh induced by Hoxd13 PAE mutations and advance our understanding of human limb development.