This study totally collected over 50 families with congenital limb malformations from outpatient clinics between 2023 and 2025. Clinical testing was performed by physical and X-ray examinations before surgical operation. Peripheral blood samples (3–5 mL) were collected from the probands, and all available family members and genomic DNA was extracted using a conventional proteinase K‐phenol‐chloroform method from peripheral blood (Li and Durbin, 2009). Among these, ten families were confirmed to have the variants of GLI3 gene and were included in this study.

Next-generation sequencing (NGS) used in this study included whole-exome sequencing (WES) and whole-genome sequencing (WGS). WES was performed in all probands, and the schedule for WES involved exome capture performed using an Agilent Sure Select Human All Exon kit (Agilent Technologies, Wilmington, DE, United States) and DNA being fragmented into approximately 150–300 base-pair (bp) lengths by ultrasonic treatment. Subsequently, a capture DNA library for sequencing was established using the HiSeq 2000 platform (Illumina Inc., San Diego, CA, United States) to a mean coverage of ×100. Raw data (FASTQ files) were generated using the Illumina base-calling pipeline. For Proband 10 in whom conventional pathogenic variants have been ruled out by WES, WGS was conducted to identify genomic structural variations and deep intronic variations. WGS was conducted thus: genomic DNA was randomly sheared to approximately 350 bp using a Covaris shearing instrument. The sheared DNA underwent end repair and end blunting to ensure 5′-end phosphorylation. The repaired products were purified using Agencourt SPRIselect beads and size-selected with the Agencourt SPRIselect nucleic acid fragment selection kit to enrich fragments of the target length. The purified, end-blunted DNA fragments underwent A-tailing to add an A base at the 3′end, enabling ligation with the T-overhang of the sequencing adapter. The completed library, after adapter ligation, underwent PCR amplification for enrichment followed by quality control. Based on effective library concentration and data volume requirements, final sequencing was performed on the Illumina NovaSeq6000 platform.

The raw data obtained from next-generation sequencing were aligned to the human reference genome (hg19) using Burrows–Wheeler Aligner backtrack software (Li and Durbin, 2009). File sorting and duplicate marking were performed with the Sequence Alignment/Map (SAM) tools (Li et al., 2009) and Picard (http://broadinstitute.github.io/picard), respectively. Variant detection was conducted using the SAM tools “mpileup” command and BCF tools (http://samtools.github.io/bcftools/bcftools.html), which identified single-nucleotide polymorphisms (SNPs) and insertion–deletion variants (indels). Annotation of variant positions, types, conservation predictions, and other relevant information was performed through ANNOVAR software (Wang et al., 2010) by integrating multiple databases, including the SNP database (dbSNP), 1000 Genomes Project, the Exome Aggregation Consortium (ExAC; http://exac.broadinstitute.org/), and the Human Gene Variant Database (HGMD; http://www.hgmd.cf.ac.uk/ac/index.php). Based on the filtering strategy, synonymous single-nucleotide variants with a minor allele frequency (MAF) > 1% in the 1000 Genomes database were excluded from the candidate pathogenic variants.

PCR-Sanger sequencing was used for validation of candidate variants revealed by NGS and family co-segregation of the phenotype and genotype. Genomic DNA/DNA (cDNA) and protein reference sequences of GLI3 (hg19: NM_000168.6/NC_000007.14 and NP_000159.3) were obtained from the University of California, Santa Cruz (UCSC) Genome Browser database (http://genome. ucsc. edu/) and the National Center for Biotechnology Information (NCBI) Reference Sequence Project. Primers were designed via the Primer3 online tool (http://primer3.ut.ee/) and were then inspected using the UCSC Genome Browser BLAT and in silico PCR online tools (Table 1). TaKaRa LA Taq^® with GC Buffer (TaKaRa, Shiga, Japan) was used in all PCR volumes, with initial denaturation at 95 C for 3 min, followed by 35 cycles of 94 C for 30 s, 58 C for 40 s, and 72 C for 40 s, with a final extension at 72 C for 8 min. Sequencing results from Applied Biosystems 3730 × l DNA Analyzer (Thermo Fisher Scientific, Waltham, MA, United States) were aligned to reference sequence using Codon Code Aligner (v.6.0.2.6; Codon Code, Centerville, MA, United States).

The large fragment deletion (Family 10) found by WGS was verified by real-time quantitative PCR (qPCR) and Gap-PCR. qPCR was performed in a saturated fluorescence dye method using TB green (TaKaRa, Dalian, China) in a 96-well Thermal iCycler (Bio-Rad, Hercules, CA, United States). The qPCR condition was as follows: initial denaturation at 95 C for 30 s, followed by 39 cycles of 95 C for 5 s (denaturation), 58 C for 30 s (annealing), and 72 C for 5 s (extension, optimized for Taq enzyme activity). Subsequently, a melting curve analysis was conducted to assess amplification specificity: 60 C for 1 min, then 60–95 C on 0.5 C per cycle (continuous fluorescence collection), and final hold at 40 C for 30 s. Based on gross deletion detected by qPCR, the breakpoint was identified by Gap-PCR combined by a 2% agarose gel electrophoresis and DNA sanger sequencing. The primers used for qPCR and Gap-PCR are listed in Table 2.

The VarCards2 (https://genemed.tech/varcards2) platform was used to predict the potential impact of all candidate variants using more than 30 algorithms, including the tools of allele frequency, pathogenicity prediction scores, and databases of disease variants. Variant pathogenicity was assessed based on ACMG guidelines. For the missense variants, evolutionary conservation across multiple species was analyzed with MEGA-11 software (Molecular Evolutionary Genetics Analysis software); the spatial structure of the GLI3 protein was predicted with Swiss model by comparing wild-type and mutant molecules, and the three-dimensional structure of the proteins was displayed in Swiss-Pdb Viewer software. For frameshift variants, DNAMAN (version 5.2.2; Lynnon Bio soft) was used to predict the possible protein truncation.

Ten families had limb malformations caused by GLI3 variants. The pedigrees and clinical manifestations are presented in Figure 1 and Table 3. Based on whether there are disease manifestations beyond limb malformations, these patients were classified into two categories: syndromic and non-syndromic. Proband 1 was an abandoned infant diagnosed with Pallister–Hall syndrome (PHS), presenting with bilateral polydactyly, cryptorchidism, and an incomplete cleft lip/palate. The remaining nine probands displayed non-syndromic limb malformations. Proband 2 was a fetus of induced termination, with a unilateral broad thumb and bilateral central polydactyly of the feet and displayed clinical manifestations akin to those of his mother (blood sample not available). Proband 3 presented with bilateral preaxial polydactyly of the hands (surgically excised) and syndactyly of the second and third toes. Probands 4 and 5 had postaxial polydactyly (PAP-A/B). Proband 6 exhibited complex polydactyly–syndactyly solely in both feet, but his mother exhibited brachydactyly of both thumbs and syndactyly between the third and fourth fingers on her left hand (preoperative photographs not available). Probands 7 and 8 both presented with Wassel type V preaxial polydactyly (PPD-IV), but Proband 8 also presented with bilateral postaxial polydactyly. Isolated PPD-IV is often categorized within the GCPS spectrum as a subtype characterized by mild craniofacial features that can be difficult to distinguish from unaffected individuals. Probands 9 and 10 displayed foot characteristics similar to those of Proband 6, characterized by complex polysyndactyly, but Proband 9 exhibited syndactyly in both hands without polydactyly, whereas Proband 10 showed both preaxial polysyndactyly and postaxial polydactyly in both hands.

Ten heterozygous GLI3 variants were screened by NGS combined with Sanger DNA sequencing. These variants included two missense (c.1063G>A; c.1489C>A), four nonsense (c.2374C>T; c.2008C>T; c.1096C>T; c.2029C>T), three frameshift (c.600delC; c.1880_1881del; c.811_812delCT) (Figure 2A), and one gross deletion (NC_000007.14: g.42061081_42069739del), and each patient was a heterozygote of one pathogenic GLI3 variant (Table 3). The familial samples of Probands 1, 2, and 3 were not available, so the co-segregation of pathogenic variant genotypes and phenotypes cannot be established in these families. In the other seven pedigrees, co-segregation of genotype and phenotype was observed. The pathogenic alleles of Probands 4, 5, and 6 were inherited from their mothers, who exhibited hand and foot malformations. The pathogenic allele in Proband 8 originated from his father, who also showed hand and foot malformations. Probands 7, 9, and 10 were sporadic cases and caused by de novo variants (c.1489C>A, c.600delC, g.42061081_42069739del, respectively). qPCR was employed to compare the copy numbers of exons 6-10 of GLI3 between Proband 10 and his normal parents. The results showed that the dosage of exons 7–9 in Proband 10 was merely half of that in normal parents. Gap-PCR was utilized to pinpoint the breakpoint of the gross deletion, uncovering a deletion of 8,658 bp in this patient. Gap-PCR can specifically amplify a 239 bp DNA fragment in cases of gross deletion in GLI3, whereas no such fragment was observed in his parents (Figure 2B). In the ten variants identified above, three (c.1880_1881del, c.2374C>T, and c.1096C>T) had been recorded by HGMD (Professional 2024.4), while the other seven (c.1063G>A, c.600delC, c.2029C>T, c.1489C>A, c.2008C>T, c.811_812delCT, NC_000007.14: g.42061081_42069739del) were novel and reported for the first time.

The evidence codes applied and the final pathogenicity classification for each variant are summarized in Table 3. According to database gnomAD v2.1.1, GLI3 is intolerant of protein-truncating variations (pLI = 1.00, oe = 0.09), indicating that loss of function is a pathogenic mechanism of the truncating variants of GLI3. Briefly, null variants (nonsense, frameshift, large deletion) were classified as pathogenic or likely-pathogenic, supported by very strong (PVS1) and supporting (PS4, PM2, and PP1) evidence. The two missense variants were assessed using computational prediction (PP3), allele frequency (PM2), and segregation (PP1) data, resulting in a classification of likely pathogenic.

GLI3 is a transcriptional activator, its structure is shown in Figure 3. Variant p. Tyr200* located within the repressor domain is predicted to result in a truncated protein lacking functional repressor and subsequent domains, including the zinc finger motif, proteolytic cleavage site, and activation domain, thereby impairing transcriptional repression. The p. Leu271Serfs*5 variant occurs within the SUFU-binding sequence (SUFU) and is likely to disrupt GLI3–SUFU interaction and subsequent processing. Variants p. Val355Ile and p. Arg366*, situated near the junction between the SUFU-binding sequence and the zinc finger domain, may affect repressor function. Variants p. Leu497Ile and p. His627Argfs*48 within the zinc finger domain are expected to impair DNA binding and transcriptional regulation. Variants p. Gln670* and p. Gln677 located close to the proteolytic cleavage site may interfere with normal protein processing. The p. Arg792 variant lies between the cleavage site and the CREB-binding domain and may reduce CBP-binding affinity. The large deletion (NC_000007.14: g.42061081_42069739del) is predicted to cause haploinsufficiency through protein truncation.

Protein conservation analysis and a three-dimensional structure prediction were also conducted for the two missense variants (Figure 4). The conservation analysis of variant sites using multiple sequence alignments indicated that GLI3 variants c.1063G>A (p.Val355Ile) and c.1489C>A (p.Leu497Ile) were highly conserved across different species (monkey, mouse, elephant, chicken, rabbit, and cow). The two missense variants were analyzed by three protein deleteriousness prediction software (Table 4). The CADD scores for the two missense variants are 22.4 and 25.2, indicating relatively strong predictive scores for harmfulness. Variants c.1063G>A (p.Val355Ile) and c.1489C>A (p. Leu497Ile) were predicted to be potentially deleterious and disease-causing (Polyphen-2 and Variant Taster, respectively). Variant p. Val355Ile may disrupt the packing of the local hydrophobic core by introducing a larger side chain, indirectly affecting the hydrogen bonding network and structural stability. Variant p. Leu497Ile induces alterations in the α-helices, side chains, and hydrogen bonds compared with normal condition, directly interfering with the spatial conformation of the GLI3 protein and decreasing the affinity of its ZNF-domain-binding DNA.