Chao Fang
Ling-Ling Liu
Wu-Jun Liu
Frederic Farnir- 1College of Animal Science, Xinjiang Agricultural University, Urumqi, China
- 2Faculté de Médecine Vétérinaire, Université de Liège Quartier Vallée, Liège, Belgium
Mandibular prognathism is a craniofacial trait that negatively affects feed intake in sheep and thus influences body weight and production performance. A genome-wide association study (GWAS) was conducted to identify genetic variants associated with mandibular prognathism (MP) in Duolang sheep. Phenotypes were classified into a binary trait: normal and MP. The analysis was based on whole-genome resequencing data from 221 individuals. We identified 48 potentially associated SNPs and 77 candidate genes related to mandibular prognathism. Based on gene functional analysis, we found five genes (PAX7, PLXNA4, CHD2, DZIP1, FBXW7) that may be involved in mandibular prognathism.Notably, seven SNPs on chromosome 4 were located within the PLXNA4 gene, and seven SNPs on chromosome 17 were located upstream of the FBXW7 gene. These findings provide novel insights into the genetic architecture of mandibular prognathism in sheep and offer potential molecular targets for future breeding and selection programs.
1 Introduction
Mandibular prognathism (MP), commonly referred to as underbite, is a craniofacial developmental anomaly characterized by anterior displacement of the mandible relative to the cranial base, resulting in malocclusion (1, 2). Clinically, according to Angle’s classification, it is categorized as class III malocclusion, in which the mesiobuccal cusp of the upper first molar occludes distal to the mesiobuccal groove of the lower first molar (3). MP may occur with or without maxillary hypoplasia (4). This disproportionate mandibular overgrowth relative to the craniofacial skeleton, historically termed the “Habsburg jaw,” which accounts for approximately 63–73% of all class III malocclusions, can cause functional impairment and reduced feeding efficiency in both humans and animals (5, 6).
The prevalence of MP shows marked variation across populations and species. In humans, it is most frequent among Asian populations (up to 15%), moderate in African populations (10–16.8%), and relatively rare among Caucasians (1%) (7). Population-level variation suggests a substantial genetic contribution. Environmental influences likely modulate dental/occlusal features rather than basal mandibular length. In livestock, jaw deformities have also been documented: for instance, the incidence of brachygnathia inferior in Brown Swiss cattle ranges from 0.1–5% (8), while an incidence of 29.5% of mandibular prognathism has been reported in Swiss Black Brown Mountain sheep (9). In wild ruminants, the frequency is typically below 1% (10). These differences highlight the complex genetic and environmental contributions underlying MP.
Research on craniofacial malformations in sheep has been sparse and mainly limited to descriptive and morphometric studies. Previous reports have largely focused on mandibular shortening in East Friesian sheep and their crossbreds, where affected individuals displayed distinct brachygnathia inferior phenotypes (10–12). However, no systematic investigations of mandibular prognathism have been conducted in native sheep breeds from China.
The Duolang sheep breed is one such population and generally produces one lamb per ewe. It is mainly distributed in Maigaiti, Bachu, Yuepuhu, and Shache counties of Xinjiang Province, where it has evolved specific adaptations to arid desert conditions (13). Despite its economic importance, genetic studies in Duolang sheep remain limited. In recent years, genome-wide association studies (GWAS) in sheep have predominantly focused on economically important traits, such as growth and body size, reproductive performance, wool and cashmere quality, and health-related resilience. For instance, Li et al. identified candidate genes associated with growth traits in Hu sheep (14), while Liu et al. reported loci related to body-size traits in Tibetan sheep (15). Beyond production-related traits, GWAS has also been extended to longevity and reproductive traits (16). Nevertheless, systematic GWAS addressing craniofacial abnormalities, particularly mandibular prognathism, remain extremely scarce in sheep. Here, we conducted a GWAS comparing normal and MP phenotypes in sheep from Kashi, Xinjiang, to identify candidate genes associated with this trait. Our findings provide new insights into the genetic basis of MP, offering potential targets for future therapeutic interventions and genetic improvement strategies.
2 Materials and methods
2.1 Sample collection and phenotyping
A total of 242 blood samples were collected from sheep in Kashi, Xinjiang, China. Each animal was screened for mandibular prognathism (MP, commonly referred to as underbite). The trait was coded as binary (1 = affected; 0 = unaffected). Phenotypic features are shown in Figure 1. MP was assessed in the field using a simple, standardized procedure: one handler restrained the head, while a trained examiner checked the occlusion in natural intercuspation with a neutral head position. A case was recorded when the lower incisors lay anterior to the upper incisors on direct inspection, consistent with the study criteria, while upper and lower incisors were aligned in controls. Husbandry practices were routine and comparable between groups. No targeted nutritional interventions were reported or used. Sampling took place in the morning and was performed by the same team following the standard protocol. Due to oversight by the sampling team, the phenotypes of 21 individuals could not be determined; thus, 221 individuals were included in the GWAS analysis.
Figure 1. Phenotypic characteristics of Duolang sheep. (A) MP phenotype. (B) Normal phenotype. Reproduced from (61), licensed under CC-BY-4.0.
2.2 DNA extraction and sequencing
Genomic DNA was extracted using a standard phenol-chloroform protocol. DNA quality was evaluated via 1% agarose gel electrophoresis, and concentrations were quantified using a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA). Purified DNA samples were shipped on dry ice to Novogene Bioinformatics Institute (Beijing, China). Following fragmentation, adapter ligation, and PCR amplification, 350-bp libraries were constructed for each individual. Sequencing was performed on an Illumina NovaSeq 6,000 platform (Illumina Inc., San Diego, CA, USA) in 2 × 150-bp paired-end mode. Libraries were sequenced to a target depth of ~10 × per individual.
2.3 Read mapping and variant detection
Sequenced reads were processed using Trimmomatic v0.39 (17) with the parameters “LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:90” to remove low-quality bases. The filtered reads were aligned to the reference genome assembly (ARS-UI_Ramb_v3.0) using the BWA-MEM algorithm v0.7.17-r1188 (18) with default settings. The resulting alignment files were coordinate-sorted by chromosome using Samtools v1.17 (19), duplicates were marked using the MarkDuplicates module in GATK v4.4.0.0 (20), and mate-pair information was synchronized using GATK’s FixMateInformation module.
Both raw and clean sequencing reads were evaluated using Seqkit v2.6.127 (21) and FastQC v0.11.928 to assess base counts and quality metrics. Variant calling was performed for each individual using GATK’s HaplotypeCaller module, followed by merging variant files with the CombineGVCFs module. Joint genotyping was conducted with the GenotypeGVCFs module. Hard filtering criteria were applied using the VariantFiltration module, with thresholds as follows: for SNPs, “QD<3.0 || FS>30.0 || SOR > 4.0 || MQ<30.0 || MQRankSum<-10.0 || QUAL<50.0 || ReadPosRankSum<-5.0,” for InDels, “QD<3.0 || FS>100.0 || ReadPosRankSum <-10.0.” Subsequently, population-level filtering was applied using VCFtools v0.1.17 (22) to ensure data quality for downstream analyses. This step involved the removal of variants with excessive missing data (missing genotype rate > 20%, parameter: --max-missing 0.8) and those with a minor allele frequency (MAF) below 0.03 (−-maf 0.03). This filtering pipeline identified a total of 35,663,272 SNPs and 4,161,247 InDels.
The preliminarily quality-controlled VCF files were converted to binary format using PLINK v1.9 (23). Further quality control for autosomal SNPs was performed with PLINK v1.9. The following criteria were applied: individual call rate ≥ 80%; SNP call rate ≥ 80%; minor allele frequency ≥ 0.03; and Hardy–Weinberg equilibrium p-value ≥ 1 × 10−5.
2.4 PCA analysis
PCA was conducted using PLINK v1.9. To reduce linkage disequilibrium, SNPs were pruned with the command: “--indep-pairwise 50 5 0.2.” Pruned markers were extracted and used to compute the top 20 principal components. The resulting eigenvectors were visualized using R v4.5.0 to assess population structure.
2.5 Genome-wide association study (GWAS)
GWAS for MP was performed using PLINK v2.0. A logistic regression framework with an additive genetic model was applied to test the association between each SNP and the phenotype. Because the case–control ratio in our dataset was unbalanced, PLINK’s firth-fallback option was enabled to provide penalized maximum-likelihood estimates when complete or quasi-complete separation occurred. To control for population structure, the top 10 principal components (PC1–PC10) derived from genome-wide genotype data were included as fixed-effect covariates. The final model can be written as:
Where is the case–control status of individual i, is the SNP genotype coded additively as 0,1, or 2 copies of the minor allele, represents the kth principal component of individual i, and coefficients represent the estimated fixed effects. Firth penalized logistic regression was automatically applied to SNPs for which the standard logistic regression failed to converge. Only autosomal chromosomes were retained for GWAS after recoding to a sheep-specific 1–26 chromosome set. Additionally, a suggestive significance threshold was set at p < 1.0 × 10−5 (15). LD structure was visualized using LDBlockShow based on SNPs extracted from the GWAS-associated regions.
2.6 Gene annotation
Genomic annotation of suggestive SNPs was performed using the Ensembl Ovis aries reference genome (Oar_v3.1; release 114). The corresponding GTF annotation file (Ovis_aries. Oar_v3.1.114.gtf) was used to obtain genomic coordinates for all annotated genes. SNP chromosomal positions were first standardized to match the reference genome. For each SNP, its genomic position was compared with annotated gene intervals to determine whether it was located within a gene. SNPs falling inside gene boundaries were assigned directly to the corresponding gene. For intergenic SNPs, all genes located within a ± 200 kb window centered on the SNP position were retrieved and recorded. Due to the limited number of annotated candidate genes near the significant SNPs, conventional GO and KEGG analyses had low statistical power. In addition, 51.9% of the genes in this study were uncharacterized. Therefore, we did not rely on enrichment statistics. Instead, we performed manual biological interpretation of individual genes based on known craniofacial developmental mechanisms.
3 Results
3.1 Results of the phenotypic analysis
A total of 221 Duolang sheep were included in the genome-wide association study, consisting of 182 individuals with the prognathism phenotype and 39 phenotypically normal controls (Figure 1).
3.2 Resequencing and sequencing depth results
The whole-genome sequencing yielded 6,837.92 Gb of raw sequence data. After stringent quality filtering, 6,780.34 Gb of high-quality clean reads were obtained, with mean Q20 and Q30 values of 96.985 and 92%, respectively. The average GC content was 43.58%, and the mean sequencing depth was approximately 11×. After quality control, one sample was removed due to missing genotype data (−-mind), 220 individuals and 27,172,983 SNPs were retained for subsequent analyses. The SNP numbers for each chromosome are shown in Figure 2.
Figure 2. Distribution of SNP counts across chromosomes.
3.3 Principal component analysis (PCA)
PCA was performed to assess population structure. The analysis did not reveal distinct clustering of individuals, suggesting minimal stratification within the analyzed samples. The first two principal components explained 8.43 and 7.02% of the genetic variance, respectively (Figure 3). As mentioned above, we nevertheless included the first 10 principal components as covariates in the GWAS model to account for potential subtle population structuration.
Figure 3. PCA plot.
3.4 Genome-wide association study results
We identified 48 SNPs potentially associated with the trait (Supplementary Table S1). Gene annotation within ±200 kb of these SNPs yielded 77 genes, of which 40 were uncharacterized (Supplementary Table S2). The SNPs were unevenly distributed, with chromosome 4 carrying the most (11 SNPs), followed by chromosome 17 (7 SNPs) (Figure 4). The QQ plot showed that the observed p-values closely followed the expected distribution, with a modest deviation in the upper tail. The genomic inflation factor was λ = 1.15, indicating acceptable control of population structure and a small number of true association signals (Figure 5).
Figure 4. Manhattan plots indicate the suggestive threshold for the GWAS of MP trait.
Figure 5. QQplot.
3.5 Regional Manhattan plots and gene annotation
To further characterize the significant signals identified in the GWAS, we performed regional association analyses for the top loci on chromosomes 4 and 17 (Figure 6). On chromosome 4 (Figure 6A), we identified a cluster of 11 significant SNPs. The association signals showed a clear peak structure, which mapped directly to the PLXNA4 (Plexin A4) gene region. On chromosome 17 (Figure 6B), a total of 7 significant SNPs were detected. Notably, the genomic interval harboring these variants contains the FBXW7 gene, a key regulator of the Notch signaling pathway.
Figure 6. Regional Manhattan plot for the potentially associated region: (A) Chromosome 4 and (B) Chromosome 17.
3.6 LD analysis of two region
We then computed pairwise LD among the suggestive SNPs within each region and generated LD heatmaps (Figure 7), which revealed the LD structure and confirmed whether the associated SNPs represented a single haplotype block. Pairwise LD analysis among the 11 SNPs on chr4 and among the 7 SNPs on chr17 showed consistently high LD (mean r2 = 0.81, range 0.60–1.0 and mean r2 = 0.79, range 0.54–1.0, respectively). For the locus on chromosome 4 (Figure 7A), the heatmap revealed a distinct haplotype block structure. The leading SNPs (e.g., from 97,918,417 to 97,928,395) exhibited strong linkage disequilibrium (indicated by the red blocks), suggesting that these variants are inherited together and could therefore potentially tag the same functional mutation within the PLXNA4 gene. Similarly, on chromosome 17 (Figure 7B), a tight LD block was observed among the top associated SNPs (e.g., 5,354,228 to 5,360,940). These high degrees of linkage (r2 approaching 1) might contribute to narrow down the genomic interval bearing putative causal variants.
Figure 7. LD heatmap. (A) Chr 4 region. (B) Chr 17 region. The color scale represents the strength of pairwise LD (r2 value), ranging from blue (low LD) to red (high LD). The axes display the positions of the SNPs.
3.7 Gene annotation analysis
We annotated the 48 SNPs and identified a total of 77 genes. Twenty SNPs were located within annotated genes, and the remaining SNPs were positioned in upstream or downstream regions. On chromosome 4, seven SNPs were located within the PLXNA4 gene, and four SNPs were located in the intergenic region between PLXNA4 and CHCHD3. On chromosome 17, seven SNPs were located in the intergenic region between TMEM154 and FBXW7. Based on published literature, we classified the annotated genes according to their known functions. These candidate genes can be broadly grouped into four functional categories. Craniofacial, neural crest, and developmental regulators, including PAX7, PLXNA4, CHD2, CHCHD3, FBXW7, DZIP1, CEP135, and KHDRBS2, which are involved in neural crest cell fate, craniofacial patterning, major signaling pathways (e.g., Hedgehog, ERK, Notch), and cytoskeletal organization. Nervous system–related genes, such as GPR158, OPTN, UBR4, MPPED2, EXOC1, and DNAJC3, which participate in neural development, axon guidance, synaptic processes, and cellular stress responses. Epithelial, structural, and metabolic genes, including CLDN10, LAMC3, PLPP7, POMT1, PRRC2B, and AIF1L, associated with cell adhesion, epithelial integrity, extracellular matrix organization, glycosylation, energy metabolism, and cell proliferation. Poorly annotated and tRNA-related genes, such as TRNAC-GCA, and TRNAR-ACG, which remain uncharacterized.
3.8 Functional annotation of candidate genes
We searched the literature for the known functions of these genes. Some are involved in craniofacial nerve development. Others are related to the nervous system, or to epithelial cell structure and metabolism. Based on these functions, PAX7 (Paired Box Gene 7), PLXNA4 (plexin-A4), CHD2 (Chromodomain Helicase DNA Binding Protein 2), DZIP1 (DAZ interacting zinc finger protein 1) and FBXW7 (F-box and WD repeat domain-containing 7) are the most likely candidates for mandibular prognathism in Duolang sheep. PAX7 is a transcription factor belonging to the PAX gene familyand is characterized by a paired DNA-binding domain. It plays important roles in muscle formation, neural development, and craniofacial morphogenesis. Functional studies have demonstrated that PAX7 contributes to neural crest specification and craniofacial patterning. In chicken embryos, Basch et al. reported that PAX7 is required for neural crest specification during gastrulation (24). In mammals, Mansouri et al. reported that PAX7 knockout in mice leads to aberrant craniofacial neural crest cell development (25). Further, Cates et al. applied particle-based shape modeling to compare cranial base morphology between Pax7-deficient and wild-type neonatal mice, and identified systematic alterations such as increased cranial base width and anterior–inferior bending of the posterior cranial margin (26). Together, these studies indicate that PAX7 is involved in the regulation of cranial neural cell fate and craniofacial development. Mechanistically, its lineage descendants are widely incorporated into cranial neural crest–derived structures, including cranial cartilage, suggesting that disruption of PAX7 may cause abnormal cranial and mandibular development (27), In humans, PAX7 has been reported to be associated with mandibular prognathism (28).
CHD2 is a member of ATP-dependent CHD chromatin remodeling family (29, 30). The canonical function of CHD2, like other chromatin remodelers, is to regulate chromatin accessibility (31). Early evidence from a chromosomal deletion study proposed that CHD2 loss may contribute to congenital anomalies, with skeletal abnormalities noted among affected individuals (32). Subsequent clinical reviews have expanded the phenotypic spectrum, documenting patients with CHD2 variants who exhibit growth disturbances, skeletal or physical abnormalities, and, in some cases, subtle facial dysmorphism (33). Case-based studies similarly describe variable morphological findings, including craniofacial or growth-related anomalies, although such features are inconsistent and typically secondary to neurological presentations (34). Broader analyses of CHD2-related microdeletions further highlight substantial phenotypic variability, encompassing skeletal, muscular, and craniofacial manifestations across different patients (35). Previous studies have reported that variants in CHD2 are associated with developmental disorders.
DZIP1 encodes a zinc-finger protein originally identified in zebrafish as the iguana locus and acts as a core regulator of primary ciliogenesis and Hedgehog (Hh) signal transduction (36, 37). Subsequent studies established DZIP1 as an essential component of the ciliogenic pathway, required for axonemal biogenesis and basal-body function (38). To date, no studies have directly linked DZIP1 or its paralog DZIP1L to mandibular prognathism in humans or animal models. DZIP1 is a regulator of primary ciliogenesis and Hedgehog signaling (39), both of which play critical roles in cranial neural crest differentiation and craniofacial morphogenesis. Loss of dzip1 in zebrafish results in ocular coloboma and abnormal cranial neural crest–derived tissue patterning (40). These findings support a role for DZIP1 in craniofacial developmental processes.
4 Discussion
4.1 Genetic basis of mandibular prognathism in Duolang sheep
Similar domestication-driven selection on standing variation has been documented in horses, based on ancient genome time-series analyses of regulatory loci linked to behavior and body conformation (41), paralleling our findings that genetic variants affecting craniofacial morphology may have been targets of selection in domesticated sheep. In this study, we performed a GWAS using the presence or absence of MP in Duolang sheep as the phenotype. A total of 220 individuals were examined, potentially providing sufficient statistical power to detect effects of reasonable size. Defining the trait as binary simplified the complex craniofacial features. It also allowed clear distinction between individuals with different morphologies. With rigorous control of population structure and stringent statistical thresholds, we obtained reliable association signals. Notably, significant regions were detected on chromosomes 4 and 17, suggesting that MP has a clear molecular basis. A study in Scottish Blackface sheep showed that feeding level and pasture type can affect some dental traits, the authors noted that feeding conditions may change eruption timing or tooth size, but they do not meaningfully alter occlusal relationships, especially the anteroposterior jaw position (42). Taken together, these observations suggest that while environmental conditions may modulate dental characteristics, mandibular prognathism itself is more likely to reflect underlying developmental and genetic mechanisms rather than postnatal environmental effects.
The identification of a limited set of candidate genes in our GWAS is consistent with this characteristic. Several candidate genes detected in this study overlap with known craniofacial developmental genes in humans and mice, which might add support to our findings. MP is not only an anatomical feature but may also influence feeding efficiency, mastication, and overall health in animals. The key genes identified in this study could serve as potential molecular markers for future breeding programs. In addition, our results provide genetic evidence for the evolution of jaw morphology in livestock. Since some candidate genes identified here are also implicated in craniofacial development in humans and mice, this suggests that mandibular prognathism in sheep shares conserved genetic mechanisms with jaw development in other mammals.
4.2 Functional interpretation of PLXNA4 and FBXW7 identified by GWAS
Among the candidate genes identified in this study, PLXNA4 and FBXW7 stand out as biologically coherent and may be highly relevant to mandibular prognathism, as both genes participate in conserved developmental pathways that link neural patterning, tissue organization and skeletal growth across vertebrate species. PLXNA4 encodes the semaphorin receptor plexin-A4, a central component of the SEMA3A–NRP1–PLXNA4 signaling axis that governs axonal guidance and cranial nerve organization (43, 44). Functional evidence from multiple model organisms highlights the evolutionary conservation of this pathway. In mice, loss of PLXNA4 leads to defasciculation of cranial branchiomotor nerves, including facial nerves such as the greater superficial petrosal nerve (GSPN), resulting in aberrant midline crossing due to disruption of SEMA3A–NRP1–PLXNA4 signaling. Additional studies have shown that plexin-A4 is expressed in the ventral posteromedial nucleus (VPM) of the thalamus, where its deficiency disrupts the organization of thalamocortical projections, highlighting its essential role in cortical somatosensory patterning (45). In zebrafish embryos, loss of plexin-A4 reduces the branching of Rohon–Beard sensory neurons and trigeminal ganglion axons, suggesting that plexin-A4 also promotes axonal branching and elaboration (46). Collectively, these studies indicate that PLXNA4 is involved in semaphorin-mediated axonal guidance and neuronal patterning. The identification of a strong association signal spanning PLXNA4 in Duolang sheep therefore suggests that perturbation of conserved axon guidance pathways may contribute to variation in mandibular architecture.
FBXW7 (or hCdc4), a substrate-recognition component of the SCF ubiquitin ligase complex, represents a second compelling candidate linking developmental regulation to mandibular structure (47, 48). Many studies on this gene are related to cancer (20, 49, 50). In addition to its role in tumorigenesis, FBXW7 has been implicated in multiple developmental processes. In mesenchymal progenitor cells, FBXW7 has been shown to control osteogenic and chondrogenic differentiation by targeting the endoplasmic reticulum–anchored transcription factors OASIS and BBF2H7 for proteasome-mediated degradation, and conditional deletion of FBXW7 in mouse mesenchymal cells promotes both osteogenesis and chondrogenesis in vitro (51). In articular cartilage, FBXW7 has also been implicated in cartilage homeostasis and osteoarthritis: mechanical overloading or ageing leads to reduced FBXW7 expression in chondrocytes, which is associated with enhanced chondrocyte senescence, increased catabolic marker expression and accelerated cartilage degeneration in mouse models (52, 53). Additional work in osteoarthritis models indicates that FBXW7 can modulate chondrocyte degeneration at least in part through regulation of the HIF-1α/VEGF axis (54). Germline loss-of-function variants in FBXW7 cause a multisystem neurodevelopmental syndrome characterized by intellectual disability, hypotonia and structural brain anomalies (55). Taken together, these findings suggest that FBXW7 may influence mandibular prognathism through modulation of chondrogenesis, cartilage maintenance and growth dynamics of the mandibular skeleton.
Collectively, the identification of PLXNA4 and FBXW7 in this GWAS points to a developmental framework in which neuromuscular patterning and osteochondral regulation act in concert to shape mandibular morphology.
4.3 Integration of candidate pathways underlying mandibular prognathism in Duolang sheep
The mandible derives from cranial neural crest–derived mesenchyme of the first pharyngeal (mandibular) arch, and develops through Meckel’s cartilage and subsequent ossification (56), these processes are orchestrated by a network of signalling pathways, including Hedgehog (Hh), BMP/FGF, Notch, semaphorin–plexin guidance cues, which together control CNC survival, patterning and osteochondrogenic differentiation during craniofacial morphogenesis.
The Notch signalling pathway is associated with mandibular prognathism and related craniofacial morphology (57, 58). NOTCH3 is one of the substrates of FBXW7 (59). Loss or mutation of FBXW7 can lead to sustained activation of Notch signalling (60). Therefore, we hypothesize that FBXW7 may regulate mandibular growth through the Notch signalling pathway.
Our candidates may be conceptualized as lying along a single developmental axis in which PAX7 and CHD2 are thought to influence the specification and epigenetic priming of cranial neural crest cells, DZIP1-dependent ciliary HH signaling may help sustain their survival and early mandibular patterning, PLXNA4-mediated Semaphorin–Plexin cues are likely to contribute to the neuromuscular loading environment experienced by the jaw, and FBXW7-dependent control of Notch and cell-cycle dynamics in osteochondral cells may modulate matrix production and remodeling. Within such a multi-layered CNCC–HH–neuromuscular–osteogenic framework, small allelic effects at each node could act in concert to bias mandibular growth vectors, potentially offering a biologically plausible explanation for the polygenic associations detected in our GWAS.
4.4 Limitations
Due to the lack of mandibular bone samples, GWAS-identified candidate genes could not be experimentally validated, and gene-level mechanistic interpretations therefore remain provisional. To address this, future work will prioritize sheep-focused validation, including use a quantitative measure rather than a binary one, improve the statistical model, use haplotypes and spatial expression mapping in embryonic mandibular tissues (e.g., in situ hybridization or immunolabeling), perturbation studies in ovine cells or embryos (e.g., targeted editing within ethical and welfare constraints), and larger population genetic analyses with fine-mapping and replication.
5 Conclusion
This study provides new insights into the genetic basis of MP in Duolang sheep. We identified five genes that may be associated with MP trait. Future studies should focus on functional validation of these candidate genes, investigation of gene–environment interactions, and cross-validation in different populations to further elucidate the complex genetic mechanisms underlying MP.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://identifiers.org/ncbi/insdc.sra:SRP544918 and https://identifiers.org/ena.embl:PRJEB83806.
Ethics statement
The animal study was approved by the Animal Welfare and Ethics Committee of Xinjiang Agricultural University, Urumqi, China (Ethics Approval ID: 2021103). The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
CF: Formal analysis, Investigation, Visualization, Writing – original draft. L-LL: Data curation, Formal analysis, Methodology, Writing – original draft, Writing – review & editing. W-JL: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing. FF: Methodology, Visualization, Writing – review & editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This research was funded by the National Natural Science Foundation of China Project, grant number 32160775, Key R&D Program of Xinjiang Uygur Autonomous Region, grant number 2021YFD1600702, and the Tianshan Talent Program Project, grant number 2023TSYCLJ0017.
Acknowledgments
The authors would like to thank Guo-Zhi Sun at the Animal Diseases Research Institute for technical assistance, Qingyong Guo and Jianlong Li for expertise and advice.
Conflict of interest
The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that Generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2025.1719178/full#supplementary-material
References
1. Kantaputra, PN, Pruksametanan, A, Phondee, N, Hutsadaloi, A, Intachai, W, Kawasaki, K, et al. ADAMTSL1 and mandibular prognathism. Clin Genet. (2019) 95:507–15. doi: 10.1111/cge.13519,
2. Sun, R, Wang, Y, Jin, M, Chen, L, Cao, Y, and Chen, F. Identification and functional studies of MYO1H for mandibular prognathism. J Dent Res. (2018) 97:1501–9. doi: 10.1177/0022034518784936,
3. Ahmed, MK, Ye, X, and Taub, P. Review of the genetic basis of jaw malformations. J Pediatr Genet. (2016) 5:209–19. doi: 10.1055/s-0036-1593505
4. Doraczynska-Kowalik, A, Nelke, KH, Pawlak, W, Sasiadek, MM, and Gerber, H. Genetic factors involved in mandibular prognathism. J Craniofacial Surg. (2017) 28:e422–31. doi: 10.1097/SCS.0000000000003627,
5. Chang, H-P, Tseng, Y-C, and Chang, H-F. Treatment of mandibular prognathism. J Formos Med Assoc. (2006) 105:781–90. doi: 10.1016/S0929-6646(09)60264-3,
6. Fang, H, Li, P, Zhu, S, and Bi, R. Genetic factors underlying mandibular prognathism: insights from recent human and animal studies. Mamm Genome. (2025) 36:293–305. doi: 10.1007/s00335-024-10084-x,
7. Tassopoulou-Fishell, M, Deeley, K, Harvey, EM, Sciote, J, and Vieira, AR. Genetic variation in myosin 1H contributes to mandibular prognathism. Am J Orthod Dentofacial Orthop. (2012) 141:51–9. doi: 10.1016/j.ajodo.2011.06.033,
8. Widmer, S, Seefried, FR, Häfliger, IM, Signer-Hasler, H, Flury, C, and Drögemüller, C. WNT10B: a locus increasing risk of brachygnathia inferior in Brown Swiss cattle. J Dairy Sci. (2023) 106:8969–78. doi: 10.3168/jds.2023-23315,
9. Greber, D, Doherr, M, Drögemüller, C, and Steiner, A. Occurrence of congenital disorders in Swiss sheep. Acta Vet Scand. (2013) 55:27. doi: 10.1186/1751-0147-55-27
10. Eriksen, T, Ganter, M, Distl, O, and Staszyk, C. Cranial morphology in the brachygnathic sheep. BMC Vet Res. (2016) 12:8. doi: 10.1186/s12917-016-0634-7
11. Kerkmann, A, Kuiper, H, Ganter, M, and Distl, O. Review of literature and results from test matings of east Friesian milk sheep affected with brachygnathia inferior. Berl Munch Tierarztl Wochenschr. (2008) 121:292–305. doi: 10.2376/0005-9366-121-292,
12. Eriksen, T, Kuiper, H, Pielmeier, R, Ganter, M, Distl, O, and Staszyk, C. Ovine craniofacial malformation: a morphometrical study. Res Vet Sci. (2012) 93:1122–7. doi: 10.1016/j.rvsc.2012.03.010,
13. Fang, C, Druet, T, Cao, H, Liu, W, Chen, Q, and Farnir, F. Whole genome sequences of 297 Duolang sheep for litter size. Sci Data. (2025) 12:1086. doi: 10.1038/s41597-025-05448-0
14. Li, T, Xing, F, Zhang, N, Chen, J, Zhang, Y, Yang, H, et al. Genome-wide association analysis of growth traits in Hu sheep. Genes (Basel). (2024) 15:1637. doi: 10.3390/genes15121637
15. Liu, D, Li, X, Wang, L, Pei, Q, Zhao, J, Sun, D, et al. Genome-wide association studies of body size traits in Tibetan sheep. BMC Genomics. (2024) 25:739. doi: 10.1186/s12864-024-10633-3
16. Smitchger, JA, Taylor, JB, Mousel, MR, Schaub, D, Thorne, JW, Becker, GM, et al. Genome-wide associations with longevity and reproductive traits in U.S. rangeland ewes. Front Genet. (2024) 15:1398123. doi: 10.3389/fgene.2024.1398123
17. Bolger, AM, Lohse, M, and Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170,
18. Li, H, and Durbin, R. Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics. (2009) 25:1754–60. doi: 10.1093/bioinformatics/btp324,
19. Li, H, Handsaker, B, Wysoker, A, Fennell, T, Ruan, J, Homer, N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. (2009) 25:2078–9. doi: 10.1093/bioinformatics/btp352,
20. McKenna, A, Hanna, M, Banks, E, Sivachenko, A, Cibulskis, K, Kernytsky, A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. (2010) 20:1297–303. doi: 10.1101/gr.107524.110,
21. Shen, W, Le, S, Li, Y, and Hu, F. SeqKit: a cross-platform and ultrafast toolkit for FASTA/Q file manipulation. PLoS One. (2016) 11:e0163962. doi: 10.1371/journal.pone.0163962,
22. Danecek, P, Auton, A, Abecasis, G, Albers, CA, Banks, E, DePristo, MA, et al. The variant call format and VCFtools. Bioinformatics. (2011) 27:2156–8. doi: 10.1093/bioinformatics/btr330,
23. Purcell, S, Neale, B, Todd-Brown, K, Thomas, L, Ferreira, MAR, Bender, D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. (2007) 81:559–75. doi: 10.1086/519795,
24. Basch, ML, Bronner-Fraser, M, and García-Castro, MI. Specification of the neural crest occurs during gastrulation and requires Pax7. Nature. (2006) 441:218–22. doi: 10.1038/nature04684,
25. Mansouri, A, Stoykova, A, Torres, M, and Gruss, P. Dysgenesis of cephalic neural crest derivatives in Pax7 − / − mutant mice. Development. (1996) 122:831–8. doi: 10.1242/dev.122.3.831,
26. Cates, J, Nevell, L, Prajapati, SI, Nelon, LD, Chang, JY, Randolph, ME, et al. Shape analysis of the basioccipital bone in Pax7-deficient mice. Sci Rep. (2017) 7:17955. doi: 10.1038/s41598-017-18199-9
27. Murdoch, B, DelConte, C, and García-Castro, MI. Pax7 lineage contributions to the mammalian neural crest. PLoS One. (2012) 7:e41089. doi: 10.1371/journal.pone.0041089,
28. da Fontoura, CSG, Miller, SF, Wehby, GL, Amendt, BA, Holton, NE, Southard, TE, et al. Candidate gene analyses of skeletal variation in malocclusion. J Dent Res. (2015) 94:913–20. doi: 10.1177/0022034515581643,
29. Alendar, A, and Berns, A. Sentinels of chromatin: chromodomain helicase DNA-binding proteins in development and disease. Genes Dev. (2021) 35:1403–30. doi: 10.1101/gad.348897.121,
30. Li, W, and Mills, AA. Architects of the genome: CHD dysfunction in cancer, developmental disorders and neurological syndromes. Epigenomics. (2014) 6:381–95. doi: 10.2217/epi.14.31,
31. Clapier, CR, Iwasa, J, Cairns, BR, and Peterson, CL. Mechanisms of action and regulation of ATP-dependent chromatin-remodelling complexes. Nat Rev Mol Cell Biol. (2017) 18:407–22. doi: 10.1038/nrm.2017.26,
32. Kulkarni, S, Nagarajan, P, Wall, J, Donovan, DJ, Donell, RL, Ligon, AH, et al. Disruption of chromodomain helicase DNA binding protein 2 (CHD2) causes scoliosis. Am J Med Genet A. (2008) 146A:1117–27. doi: 10.1002/ajmg.a.32178,
33. Wilson, M-M, Henshall, DC, Byrne, SM, and Brennan, GP. CHD2-related CNS pathologies. Int J Mol Sci. (2021) 22:588. doi: 10.3390/ijms22020588,
34. Zhu, L, Peng, F, Deng, Z, Feng, Z, and Ma, X. A novel variant of the CHD2 gene associated with developmental delay and myoclonic epilepsy. Front Genet. (2022) 13:761178. doi: 10.3389/fgene.2022.761178,
35. De Maria, B, Balestrini, S, Mei, D, Melani, F, Pellacani, S, Pisano, T, et al. Expanding the genetic and phenotypic spectrum of CHD2-related disease: from early neurodevelopmental disorders to adult-onset epilepsy. Am J Med Genet A. (2022) 188:522–33. doi: 10.1002/ajmg.a.62548,
36. Sekimizu, K, Nishioka, N, Sasaki, H, Takeda, H, Karlstrom, RO, and Kawakami, A. The zebrafish iguana locus encodes Dzip1, a novel zinc-finger protein required for proper regulation of hedgehog signaling. Development. (2004) 131:2521–32. doi: 10.1242/dev.01059,
37. Wolff, C, Roy, S, Lewis, KE, Schauerte, H, Joerg-Rauch, G, Kirn, A, et al. iguana encodes a novel zinc-finger protein with coiled-coil domains essential for hedgehog signal transduction in the zebrafish embryo. Genes Dev. (2004) 18:1565–76. doi: 10.1101/gad.296004,
38. Tay, SY, Yu, X, Wong, KN, Panse, P, Ng, CP, and Roy, S. The iguana/DZIP1 protein is a novel component of the ciliogenic pathway essential for axonemal biogenesis. Dev Dyn. (2010) 239:527–34. doi: 10.1002/dvdy.22199,
39. Kim, HR, Richardson, J, van Eeden, F, and Ingham, PW. Gli2a protein localization reveals a role for Iguana/DZIP1 in primary ciliogenesis and a dependence of hedgehog signal transduction on primary cilia in the zebrafish. BMC Biol. (2010) 8:65. doi: 10.1186/1741-7007-8-65
40. Nandamuri, SP, Lusk, S, and Kwan, KM. Loss of zebrafish dzip1 results in inappropriate recruitment of periocular mesenchyme to the optic fissure and ocular coloboma. PLoS One. (2022) 17:e0265327. doi: 10.1371/journal.pone.0265327,
41. Liu, X, Jia, Y, Pan, J, Zhang, Y, Gong, Y, Wang, X, et al. Selection at the GSDMC locus in horses and its implications for human mobility. Science. (2025) 389:925–30. doi: 10.1126/science.adp4581,
42. Purser, AF, Wiener, G, and West, DM. Causes of variation in dental characters of Scottish blackface sheep in a hill flock, and relations to ewe performance. J Agric Sci. (1982) 99:287–94. doi: 10.1017/S0021859600030045
43. Suto, F, Murakami, Y, Nakamura, F, Goshima, Y, and Fujisawa, H. Identification and characterization of a novel mouse plexin, plexin-A4. Mech Dev. (2003) 120:385–96. doi: 10.1016/S0925-4773(02)00421-5,
44. Suto, F, Tsuboi, M, Kamiya, H, Mizuno, H, Kiyama, Y, Komai, S, et al. Interactions between plexin-A2, plexin-A4, and Semaphorin 6A control Lamina-restricted projection of hippocampal mossy Fibers. Neuron. (2007) 53:535–47. doi: 10.1016/j.neuron.2007.01.028,
45. Suto, F, Ito, K, Uemura, M, Shimizu, M, Shinkawa, Y, Sanbo, M, et al. Plexin-A4 mediates axon-repulsive activities of both secreted and transmembrane Semaphorins and plays roles in nerve Fiber guidance. J Neurosci. (2005) 25:3628–37. doi: 10.1523/JNEUROSCI.4480-04.2005,
46. Miyashita, T, Yeo, S-Y, Hirate, Y, Segawa, H, Wada, H, Little, MH, et al. PlexinA4 is necessary as a downstream target of Islet2 to mediate slit signaling for promotion of sensory axon branching. Development. (2004) 131:3705–15. doi: 10.1242/dev.01228,
47. Skaar, JR, D’Angiolella, V, Pagan, JK, and Pagano, M. SnapShot: F box proteins II. Cell. (2009) 137:1358–1358.e1. doi: 10.1016/j.cell.2009.05.040
48. Cardozo, T, and Pagano, M. The SCF ubiquitin ligase: insights into a molecular machine. Nat Rev Mol Cell Biol. (2004) 5:739–51. doi: 10.1038/nrm1471,
49. Yang, H, Lu, X, Liu, Z, Chen, L, Xu, Y, Wang, Y, et al. FBXW7 suppresses epithelial-mesenchymal transition, stemness and metastatic potential of cholangiocarcinoma cells. Oncotarget. (2015) 6:6310–25. doi: 10.18632/oncotarget.3355,
50. Xiao, G, Lu, W, Yuan, J, Liu, Z, Wang, P, and Fan, H. Fbxw7 suppresses carcinogenesis and stemness in triple-negative breast cancer through CHD4 degradation and Wnt/β-catenin pathway inhibition. J Transl Med. (2024) 22:99. doi: 10.1186/s12967-024-04897-2
51. Yumimoto, K, Matsumoto, M, Onoyama, I, Imaizumi, K, and Nakayama, KI. F-box and WD repeat domain-containing-7 (Fbxw7) protein targets endoplasmic reticulum-anchored osteogenic and chondrogenic transcriptional factors for degradation. J Biol Chem. (2013) 288:28488–502. doi: 10.1074/jbc.M113.465179,
52. Zhang, H, Shao, Y, Yao, Z, Liu, L, Zhang, H, Yin, J, et al. Mechanical overloading promotes chondrocyte senescence and osteoarthritis development through downregulating FBXW7. Ann Rheum Dis. (2022) 81:676–86. doi: 10.1136/annrheumdis-2021-221513
53. Shen, J, Fu, B, Li, Y, Wu, Y, Sang, H, Zhang, H, et al. E3 ubiquitin ligase-mediated regulation of osteoblast differentiation and bone formation. Front Cell Dev Biol. (2021) 9:706395. doi: 10.3389/fcell.2021.706395,
54. Zhu, W-J, Chang, B-Y, Wang, X-F, Zang, Y-F, Zheng, Z-X, Zhao, H-J, et al. FBW7 regulates HIF-1α/VEGF pathway in the IL-1β induced chondrocytes degeneration. Eur Rev Med Pharmacol Sci. (2020) 24:5914–24. doi: 10.26355/eurrev_202006_21484,
55. Stephenson, SEM, Costain, G, Blok, LER, Silk, MA, Nguyen, TB, Dong, X, et al. Germline variants in tumor suppressor FBXW7 lead to impaired ubiquitination and a neurodevelopmental syndrome. Am J Hum Genet. (2022) 109:601–17. doi: 10.1016/j.ajhg.2022.03.002,
56. Fabik, J, Psutkova, V, and Machon, O. The mandibular and hyoid arches-from molecular patterning to shaping bone and cartilage. Int J Mol Sci. (2021) 22:7529. doi: 10.3390/ijms22147529,
57. Han, X, Xiong, X, Shi, X, Chen, F, and Li, Y. Targeted sequencing of NOTCH signaling pathway genes and association analysis of variants correlated with mandibular prognathism. Head Face Med. (2021) 17:17. doi: 10.1186/s13005-021-00268-0
58. Gershater, E, Li, C, Ha, P, Chung, C-H, Tanna, N, Zou, M, et al. Genes and pathways associated with skeletal sagittal malocclusions: a systematic review. Int J Mol Sci. (2021) 22:13037. doi: 10.3390/ijms222313037,
59. Matsumoto, A, Onoyama, I, Sunabori, T, Kageyama, R, Okano, H, and Nakayama, KI. Fbxw7-dependent degradation of Notch is required for control of “stemness” and neuronal-glial differentiation in neural stem cells. J Biol Chem. (2011) 286:13754–64. doi: 10.1074/jbc.M110.194936,
60. Wang, Z, Inuzuka, H, Zhong, J, Wan, L, Fukushima, H, Sarkar, FH, et al. Tumor suppressor functions of FBW7 in cancer development and progression. FEBS Lett. (2012) 586:1409–18. doi: 10.1016/j.febslet.2012.03.017,
Keywords: candidate genes, Duolang sheep, enrichment analysis, GWAS, mandibular prognathism
Citation: Fang C, Liu L-L, Liu W-J and Farnir F (2026) Identification of key genes for mandibular prognathism in Duolang sheep via genome-wide association analysis. Front. Vet. Sci. 12:1719178. doi: 10.3389/fvets.2025.1719178
Edited by:
Xuexue Liu, UMR5288 Anthropologie Moleculaire et Imagerie de Synthese (AMIS), FranceCopyright © 2026 Fang, Liu, Liu and Farnir. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ling-Ling Liu, bGluZ2xpbmdsaXUxOTg4QHhqYXUuZWR1LmNu