Tissue-specific expression and functional role of keratin 1 in sheep horn development

Introduction: Sheep horn development has significant implications for animal welfare and farm management, yet its molecular mechanisms remain incompletely understood. Keratin 1 (KRT1), a key structural protein in epidermal keratinization, has been implicated in horn formation. This study systematically investigates the expression patterns, genetic variations, and functional role of KRT1 in sheep horn development.

Methods: We integrated RNA sequencing (RNA-seq) and whole-genome sequencing (WGS) data from Tibetan sheep and public databases. Multi-tissue expression profiling of KRT1 was performed across sheep, cattle, pigs, and humans. Phylogenetic and protein structural analyses identified conserved amino acid sites. Allele-specific expression (ASE) loci and functional SNPs were screened using population genetics approaches. Association analysis linked genotypes with horn length in Small Tail Han sheep.

Results: KRT1 expression was significantly higher in scurred (small-horned) sheep compared to SHE (large-horned) sheep (p = 0.024), and exhibited tissue-specific enrichment in horn, skin, and periosteum. Cross-species analysis confirmed high KRT1 expression in horn and skin tissues. We identified 11 horned-artiodactyl-specific amino acid sites, including K312, which forms a hydrogen bond with E262 of KRT10; mutations at this site disrupted the interaction. Four ASE loci showed strong bias toward reference alleles in horned phenotypes. Thirty-two functional SNPs were prioritized, and nine haplotype blocks contained 13 highly differentiated SNPs (Fst > 0.05). Four SNPs were significantly associated with horn length, with wild-type homozygotes exhibiting longer horns (p < 0.05).

Discussion: Our findings demonstrate that KRT1 plays a critical role in sheep horn development through its expression regulation, protein interaction stability, and genetic variation. The conserved K312 residue and associated SNPs may serve as potential molecular markers for horn phenotype selection. These results provide new insights into the keratin-based mechanisms underlying horn morphogenesis and offer a foundation for molecular breeding strategies aimed at horn size and type modulation in sheep.

1 Introduction

Sheep (Ovis aries) are globally important livestock, providing humans with high-quality mutton, nutritious milk, and wool, a key raw material for the textile industry (1). In intensive farming systems, maximizing production efficiency and ensuring animal welfare are core objectives. However, sheep horns—an anatomical structure composed of a bony core and an outer keratin sheath (2)—pose significant challenges in densely housed environments. In natural settings, horns play important roles in defense against predators, intraspecific competition, and sexual selection (3, 4). Under confined conditions, however, fighting among horned individuals frequently causes severe injuries, such as body surface bruising, udder damage (reducing milk yield), and deterioration of carcass quality (affecting meat quality) (5). These traumas not only directly impact individual growth, development, and reproductive performance but also increase the risk of injury to handlers during routine management operations (6). To mitigate these issues, early physical disbudding is commonly practiced in farming. However, this procedure is typically performed on lambs within 1–2 weeks after birth and often lacks effective analgesia due to technical limitations, safety concerns, and regulatory requirements, causing significant pain to the animals and constituting a serious animal welfare problem (7). Furthermore, routine disbudding has been banned in some regions, only permitted under specific conditions (8). Therefore, deciphering the molecular mechanisms of horn formation and development is crucial for enhancing animal welfare and ensuring sustainable production benefits in the livestock industry.

The KRT1 gene primarily encodes Keratin 1, a key structural protein in epidermal keratinization. Existing research indicates that KRT1 is closely associated with the formation of the keratin sheath and plays an important role in its development (9). Sheep horns, as a specialized organ, consist of an internal bony core and an external keratin sheath (2). Studies on sheep horn type differentiation have identified several key genes, including RXFP2, FOXL2, ACAN, and TNN, all of which have been confirmed to significantly influence horn phenotypes (10). Among these, the relaxin family peptide receptor 2 (RXFP2) gene, located on chromosome 10 and enriched in the neuroactive ligand-receptor interaction pathway, is considered a key gene controlling horn formation and is significantly associated with horn size and shape (11, 12). Notably, the KRT1 gene has also been identified as a candidate gene for horn development, potentially influencing sheep horn size (10). Homozygous nonsense mutations in KRT1 [e.g., c.457C>T (p. Gln153) and c.33C>G (p. Tyr11)] can lead to nonsense-mediated mRNA decay and KRT1 protein deficiency, resulting in epidermolytic palmoplantar keratoderma (EPPK) with knuckle pads (13); additionally, KRT1 mutations have been confirmed to cause epidermolytic hyperkeratosis, manifesting as ichthyosis-like hyperkeratosis (14). These findings consistently demonstrate the central role of KRT1 in epidermal keratinization. Given that horns are skin derivatives and the development of their keratin sheath is primarily regulated by skin-related pathways (15), we hypothesize that the KRT1 gene may regulate the generation, development, and final morphology of sheep horns by influencing skin keratinization pathways.

Based on this background, this study aims to systematically dissect the core regulatory role and molecular mechanism of the KRT1 gene in sheep horn morphological differentiation. To achieve this, we integrated RNA-seq transcriptome sequencing and whole-genome sequencing (WGS) data to conduct multi-dimensional analyses: (1) Investigate the differential expression patterns of the KRT1 gene in horn tissue and other key tissues (e.g., periosteum, skin, internal organs); (2) Systematically identify functional genomic variations affecting horn development [including allele-specific expression (ASE) loci, functional SNP loci, and linkage disequilibrium haplotype blocks]; (3) Conduct in-depth analysis of key amino acid loci and their mediated protein interaction networks (e.g., KRT1-KRT10 interaction). Through these integrated analyses, this study aims to elucidate the molecular pathway by which the KRT1 gene influences horn morphological development through regulating the keratin interaction network. This will provide crucial new evidence for a deeper understanding of the molecular basis of sheep horn development and lay the theoretical foundation for molecular marker-assisted breeding based on functional loci within the KRT1 gene.

2 Materials and methods

2.1 Ethical statement

The animal study protocol was approved by the Animal Ethics Committee of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS) under approval number IASCAAS-AE-03 (December 12, 2016).

2.2 Animals and sample collection

The RNA-seq dataset used in this experiment was previously collected by our research group (PRJNA1003277) (16). Soft horn tissue samples were collected from eight female Tibetan sheep (with an average age of 4.5 years) in the Tibet region of China. Samples were divided into two groups: the scurred group (small horns, n = 4, horn length 0–12 cm) and the SHE group (large horns, n = 4, horn length >12 cm; Supplementary Table S1). The collected tissues were placed in cryotubes and flash-frozen in liquid nitrogen.

2.3 RNA-seq data analysis

To explore the expression of the KRT1 gene across tissues, 2,915 publicly available high-quality RNA-seq datasets from sheep were collected, primarily sourced from the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov/, accessed on July 23, 2023) and the European Bioinformatics Institute (EBI) database (https://www.ebi.ac.uk/, accessed on July 23, 2023). These public datasets, along with the RNA-seq data from the 8 Tibetan sheep collected previously, were processed as follows. Adapters and low-quality data were removed using Trim Galore (v.0.6.10). Clean reads were mapped to the sheep reference genome ARS-UI_Ramb_v2.0 using STAR (v.2.7.10b) with parameters –chimSegmentMin 10 and –outFilterMismatchNmax 3. High-quality RNA-seq clean data meeting the criteria of unique mapping rate >85% and number of clean reads >20,000,000 were obtained for subsequent analysis (17) (Supplementary Table S2). Transcripts Per Million (TPM) and Fragments Per Kilobase of exon per Million mapped fragments (FPKM) values were calculated using StringTie (v.2.1.5) (18). Raw gene counts were subsequently extracted using featureCounts (v.2.0.1).

2.4 KRT1 gene expression in SHE and scurred groups

To explore differences in KRT1 gene expression between the SHE and scurred groups, violin plots were generated using the ggplot2 (v.3.4.2) package in R (v.4.3.1), and significance was assessed using a T-test. To further investigate differences in exon expression of the KRT1 gene between the SHE and scurred groups, we used Subread_to_DEXSeq (https://github.com/vivekbhr/Subread_to_DEXSeq) of dexseq_prepare_annotation2.py script to format the annotation (GTF) file of the genome, then used load_SubreadOutput.R in Rstudio reads the formatted GTF file and the counts' matrix output by featureCounts, constructs the DEXSeqDataSetFromFeatureCounts (dds) object, and performs exon difference analysis (19). To further visualize the differential expression of the KRT1 gene across the SHE group, scurred group, and different tissues, the KRT1 gene track was constructed and visualized on the UCSC Genome Browser (https://genome.ucsc.edu/s/yoser4/15horn_vs_145other).

2.5 KRT1 gene expression profiling

To explore KRT1 gene expression across species, publicly available RNA-seq data from four species (sheep, pig, cattle, and human) were downloaded. RNA-seq data from 2,651 pigs were obtained from the Pig GTEx project, with TPM gene expression values sourced from http://piggtex.farmgtex.org/ (accessed 10 September 2023) (20). RNA-seq data from 4,359 cattle were obtained from the Cattle GTEx project, with TPM gene expression values sourced from https://cgtex.roslin.ed.ac.uk/ (accessed 1 October 2023) (21). TPM gene expression values for RNA-seq data from 9,810 human samples were sourced from https://gtexportal.org/home/datasets (accessed 5 August 2023) (22, 23). Sheep, pig, cattle, and human RNA-seq data were pooled to align tissue samples, which were subsequently categorized into 17 distinct tissue types. The average TPM value of the KRT1 gene for each species was calculated. Bubble plots illustrating the tissue-specific expression of the KRT1 gene across the four species were generated using the pheatmap (v.1.0.12) package in R (v.4.3.1). Additionally, to investigate sex-specific differences in tissue expression of the sheep KRT1 gene, significant differential expression between rams and ewes across 12 tissues was calculated using a T-test. To study differences in KRT1 gene expression between scurred and SHE sheep breeds, data for the seven most representative sheep breeds from the aforementioned dataset were extracted and divided into horned and hornless groups for expression analysis. The groups were defined as follows: horned Group: Minxian Black Fur Sheep, Tibetan Sheep, Tan Sheep, Chinese Merino Sheep; hornless Group: Bashibai Sheep, Texel Sheep, Hu Sheep, using multiple comparisons test.

2.6 KRT1 protein sequence analysis and 3D structure prediction

To investigate differences in the KRT1 protein among different species, protein sequences of KRT1 homologs from 18 distinct species were retrieved from the NCBI Homologene database (Supplementary Table S3). Using MEGA11 software (24), multiple sequence alignment was performed using the ClustalW method (default parameters) (25). Based on the alignment results, a circular phylogenetic tree of the KRT1 protein was constructed using the Neighbor-Joining (NJ) method (26). Differences in the KRT1 protein sequences across species were analyzed, with a particular focus on differences between horned artiodactyl species and other species, to identify specific amino acid loci. The three-dimensional structure of the KRT1 protein was predicted using AlphaFold3 (default parameters) via NCBI, and the spatial locations of the specific amino acid loci were annotated (27). To identify spatial binding loci of the KRT1 protein, a Protein-Protein Interaction (PPI) network was constructed using the STRING database (28) to predict the interaction network of the KRT1 protein in sheep. The three-dimensional structure of the interaction between the KRT1 and KRT10 proteins was predicted using AlphaFold3 (default parameters). PyMOL was used for visualization analysis (29). Combining the specific amino acid loci identified above, the specific amino acid loci mediating the interaction between KRT1 and KRT10 proteins were identified. The effect of mutations at these specific binding loci amino acids on the hydrogen bonding interaction between the two proteins was observed. The resulting mutant protein complex structures were predicted using AlphaFold3 software and visualized using Pymol.

2.7 Whole-genome sequencing (WGS) analysis

Whole-genome sequencing (WGS) data comprised publicly available data and data collected in the laboratory. Publicly available data for 200 sheep were downloaded from NCBI, including projects PRJNA304478, PRJNA325682, PRJNA47952, PRJNA624020, PRJNA675420, PRJNA822017, PRJNA30931, PRJNA480684, PRJNA509694, PRJNA779188, and PRJNA783661 (30–41). Samples were classified as horned or polled. The FST values—a measure of genetic differentiation between populations—were calculated for samples from the public data using vcftools (v.0.1.16) with default parameters. An FST value >0.05 was set as the differentiation threshold. SNPs were annotated using SnpEff and categorized based on the annotation results (e.g., exon, intron, upstream, downstream) to filter for more specific SNP loci in these sheep. Linkage disequilibrium (LD) analysis of SNP loci within the KRT1 gene region was performed, and LD heatmaps were generated using LDBlockShow (v.1.40) software (42).

To validate the WGS findings, data from 35 Small Tail Han sheep were collected (horn length trait was defined as the average of left and right horn lengths). Raw data were trimmed and filtered for low-quality sequences and adapters using Trimmomatic (v.0.39), and quality control results were assessed using FastQC (v.0.12.1) (43). Clean reads were aligned to the reference genome. Qualified reads were aligned, sorted, and duplicates marked using BWA (v.0.7.17) and Picard (v.3.1.1) against the sheep reference genome (44, 45). SNP calling was performed using the Genome Analysis Toolkit GATK (v.4.2.5.0) with default parameters, and variant effects were annotated using SnpEff (v.5.2) (46). VCF files for the KRT1 gene region were extracted using bcftools (v.1.19) and filtered using vcftools (v.0.1.16) with parameters –maf 0.05 and –max-missing 0.8 (47). Horn length-associated SNP loci within the KRT1 gene region were extracted. Violin plots were generated to analyze the association between genotype mutations at these significant loci and the horn length trait in sheep.

3 Results

3.1 Expression of KRT1 gene in sheep with divergent horn phenotypes

Differential expression analysis revealed significantly higher KRT1 gene transcript levels in the scurred group (small-horned phenotype) compared to the SHE group (large-horned phenotype; p = 0.024; Figure 1A). Exon-level resolution demonstrated elevated expression of all nine exons in scurred horn tissue, with the highest expression observed in exon 9 (E009; Figure 1B). Transcriptomic profiling via the UCSC Genome Browser (Figure 1C) delineated the expression patterns, GC content distribution, and repetitive element landscape of KRT1 gene in horn tissues. Consistent with exon quantification data, all exons exhibited stronger expression signals in scurred vs. SHE samples, with maximal intensity at E009. Tissue-specific analysis across seven ovine tissues revealed pronounced KRT1 gene expression in horn tissue, moderate levels in periosteum, lower expression in skin, and negligible detection in non-keratinized visceral organs (heart, liver, lung, kidney) and muscle. This expression hierarchy—horn tissue > periosteum > skin > non-keratinized viscera/muscle—highlights the structural congruence between horn keratinization and dermal-derived tissues. The concordant expression patterns in horn and periosteum support the hypothesis that horn morphogenesis involves skin-like keratinization processes governing both the keratinous sheath and bony core development.

Figure 1

Figure 1. (A) Violin plot depicting differential KRT1 gene expression between scurred (small-horned, n = 4) and SHE (large-horned, n = 4) horn tissues (p = 0.024, T-test). (B) Exon-specific expression of KRT1 gene. Y-axis: Fitted expression estimates from GLM regression; X-axis: Exons E001–E009. Red and blue lines denote scurred and SHE groups, respectively. (C) UCSC Genome Browser view showing KRT1 gene expression profiles, GC content, and repetitive elements in horn tissues (scurred/SHE) and seven somatic tissues. Tissue labels are indicated on the left.

3.2 The tissue-specific function of the KRT1 gene in sheep

To gain deeper insights into the cross-species expression profile of the KRT1 gene, we examined its expression in sheep, pigs, humans, and cattle. Results (Figure 2A) revealed that KRT1 gene is highly expressed in the skin tissue of all four species, with the highest expression level observed in human skin compared to bovine and ovine skin. Notably, significant KRT1 gene expression was also detected in both sheep horn and cattle horn tissues, with higher expression levels in cattle horn than in sheep horn. Utilizing publicly available RNA-seq data, we further analyzed the sex-specific expression differences of KRT1 gene across various sheep tissues (Figure 2B). The analysis demonstrated that KRT1 gene expression is markedly tissue-specific, being detectable predominantly in skin tissue. Crucially, within skin tissue, KRT1 gene expression levels were significantly higher in male individuals compared to females (p = 0.00691). Furthermore, to assess the consistency of KRT1 gene expression across diverse sheep breeds, we compared its expression levels between horned and hornless sheep breed groups (Figure 2C). Among the hornless breeds, Hu Sheep exhibited the highest KRT1 gene expression. Conversely, within the horned breeds, Minxian Black Fur Sheep showed the lowest expression. Overall, the expression level of KRT1 gene was relatively higher in the hornless breed group compared to the horned breed group, with the difference being statistically significant, aligning with its potential biological role in horn development.

Figure 2

Figure 2. (A) The difference of RNA-seq expression of KRT1 gene in different species and tissues. (B) The difference of RNA-seq expression of KRT1 gene in tissues of Female and Male (p = 0.00691). (C) The difference of RNA-seq expression of KRT1 gene in skin tissues of horned and hornless sheep.

3.3 KRT1 protein evolutionary analysis

A phylogenetic tree of KRT1 protein sequences from sheep and 17 other species was constructed using the NG method in MEGA11 software (Figure 3A). The results indicate that the KRT1 protein sequences are most closely related between sheep and other horned artiodactyls, suggesting structural similarity of the KRT1 protein among horned animals. Concurrently, the KRT1 amino acid sequences from these 18 species were compared, and 11 divergent amino acid loci specific to horned artiodactyls were identified (Figure 3C). Amino acid loci 122, 221, 282, 305, 312, 328, 421, 433, 519, 560, and 592 are identical in horned artiodactyls but differ in other mammals. These results demonstrate the structural similarity of the KRT1 protein among horned animals. These conserved loci may play crucial roles in horn growth and development. Furthermore, it was found that amino acid loci 282, 305, 312, 328, 421, and 433 are located on the α-helix of the central rod domain of the KRT1 protein, while the other loci reside on linkers, with loci 221, 519, and 560 being aromatic rings (Figure 3B).

Figure 3

Figure 3. (A) Circular phylogenetic tree of KRT1 proteins from 18 species, with inner ring colored by group and outer ring colored by family. (B) 3D structure of the KRT1 protein highlighting the specific amino acid loci. (C) Comparative diagram of homologous KRT1 amino acid sequences from 18 species, with loci specific to horned artiodactyls marked in red. (D) Protein-protein interaction network of KRT1 in sheep. (E) Schematic diagram of the interaction between loci E262 of the sheep KRT1 protein and loci K312 of the KRT10 protein. (F) Schematic diagram illustrating the disruption of the hydrogen bond between KRT1 and KRT10 due to amino acid changes at loci 312 of the KRT1 protein.

The STRING database was used to predict the protein-protein interaction (PPI) network of KRT1 with other proteins in sheep (Figure 3D). Interactions were observed between KRT1 and proteins KRT10, KRT78, UNC79, GRHL3, DSC1, DSG1, PNPLA1, TGM1, PKP1, and DSP. Among these, KRT10 showed the strongest connection to KRT1. The study also identified a specific divergent amino acid sequence region involved in the interaction between KRT1 and KRT10 proteins. Specifically, the K312 amino acid loci of the KRT1 protein interacts with the E262 loci of the KRT10 protein via a hydrogen bond with a strength of 2.9 Å (Figure 3E). Based on the divergent amino acid loci specific to horned artiodactyls (referenced from Figure 4B), position 312 in hornless animals is occupied by T or S. It was observed that substitution of K312 with either T or S disrupts the hydrogen bond between K312 and E262 (Figure 3F). This further underscores the importance of the KRT1-KRT10 interaction in horned artiodactyls. This specific amino acid loci may play a vital role in horn growth and development.

Figure 4

Figure 4. (A) Pie chart of differential ref/alt allele counts at ASE1 (chr3:133703627) between Scurred and SHE groups; blue = ref allele counts, purple =alt allele counts. (B) Pie chart of differential ref/alt allele counts at ASE1 (chr3:133705599) between Scurred and SHE groups; blue = ref allele counts, purple =alt allele counts. (C) Pie chart of differential ref/alt allele counts at ASE1 (chr3:133703666) between Scurred and SHE groups; blue = ref allele counts, purple =alt allele counts. (D) Pie chart of differential ref/alt allele counts at ASE1 (chr3:133705932) between Scurred and SHE groups; blue = ref allele counts, purple = alt allele counts.

3.4 Allele-specific expression (ASE) of KRT1 gene

Four critical allele-specific expression (ASE) loci were identified from 30 ASE (Supplementary Table S4) loci of the KRT1 gene: ASE1 (chr3:133703627), ASE2 (chr3:133703666), ASE3 (chr3:133705599), and ASE4 (chr3:133705932), all located within introns. Analysis of ref/alt allele counts at these four ASE loci in SHE and Scurred groups (Figure 4) revealed a consistent pattern: ref allele counts markedly exceeded alt allele counts in the Scurred group; conversely, alt alleles were completely absent (zero counts) in the SHE group, where only ref alleles were expressed. This systematic bias toward ref allele expression—particularly the complete absence of alt alleles in the horned phenotype—indicates strong negative selection against alt alleles at these loci during horn development. The coordinated ASE pattern across these four intronic regulatory elements suggests their collective role in regulating KRT1 expression and horn type determination.

3.5 Screening for potential functional variants in the KRT1 gene

We analyzed the population differentiation index (Fst) values of variant loci within the KRT1 gene across different sheep breeds (Figure 5A) and identified 20 loci exhibiting high selection signals (Fst > 0.05). This indicates significant population differentiation at these loci between large-horned and small-horned breeds, suggesting their potential association with horn formation or phenotypic variation. Through functional annotation of SNPs across the 199 variant loci in the KRT1 gene region, we prioritized 32 potential functional SNPs (Table 1): SNPs 1–8 are located in the upstream promoter region. SNP 32 is located in the downstream region of the gene. SNPs 13–14, 16–21, 23, 26–27 are located within intronic regions. These potential functional SNPs generally display higher Fst values, further supporting their significant differentiation between large-horned and small-horned populations. This suggests they may play key roles in the transcriptional regulatory network or gene expression regulation governing horn development. Within the exonic regions of the KRT1 gene, a total of 11 SNPs were identified (SNPs 9–12, 15, 22, 24–25, 28–31): six of these are synonymous mutation loci (SNPs 22, 24, 25, 28). These loci exhibit low Fst values, indicating minimal differentiation across populations, likely due to neutral selection. Notably, SNPs 12 and 15, despite being synonymous mutations, show relatively high Fst values. They might be involved in regulation by influencing codon usage bias or mRNA stability. Three SNP loci harbor missense mutations: SNP9 (p. Gly44Ala), SNP10 (p. Gly56Ser), SNP31 (p. Gly139Ser). Although the Fst values for these missense variants do not indicate significant differentiation, they may still exert functional effects on protein structure or function within specific populations, holding potential biological significance. Additionally, SNPs 29 to 31 are located within the 3′ untranslated region 3′ UTR. While their Fst values are not high and their population frequencies are low, their functional roles warrant further validation in larger sample cohorts.

Figure 5

Figure 5. (A) Scatter plot showing the distribution of Fst values for all variant loci across the KRT1 gene region. Loci with Fst > 0.05 (Difference) are marked in blue; loci with Fst < 0.05 (no Difference) are marked in gray. The upstream region (Upstream) is marked in purple; the downstream region (Downstream) is marked in dark blue; exonic regions (Exon) are marked in yellow. (B) Linkage disequilibrium (LD) plot for SNP loci within the KRT1 gene. Darker block colors represent higher degrees of LD (measured as r²).

Table 1

Table 1. Potential SNP functional loci of KRT1 gene.

To further explore the genetic architecture of the KRT1 gene and potential selection signals, we performed linkage disequilibrium (LD) analysis on these SNPs, identifying nine haplotype blocks. These blocks collectively encompass 13 SNPs with Fst > 0.05, specifically: SNP1 (chr3:133701816), SNP2 (chr3:133701864), SNP4 (chr3:133702294), SNP5 (chr3:133702304), SNP15 (chr3:133704641), SNP16 (chr3:133704995), SNP17 (chr3:133705002), SNP18 (chr3:133705243), SNP19 (chr3:133705250), SNP20 (chr3:133705599), SNP21 (chr3:133705611), SNP26 (chr3:133707421), SNP27 (chr3:133707432). In the LD plot (Figure 5B), the blocks representing the pairwise LD relationships among these SNPs show significantly darker coloration, indicating strong LD. These regions likely constitute genetic units closely associated with horn phenotype traits.

3.6 Association analysis between horn length and genotypes

To investigate loci significantly associated with horn length within the KRT1 gene region, we performed association analysis using the average length of the left and right horns from 35 Small-tailed Han sheep individuals as the horn length phenotype. This analysis identified four SNPs significantly associated with horn length: chr3:133703518, chr3:133704023, chr3:133705250, and chr3:133706677 (Supplementary Table S5). At each of these four loci, individuals carrying the wild-type homozygous genotype exhibited significantly greater mean horn length compared to those carrying the heterozygous genotype (p < 0.05) (Figure 6): At chr3:133703518, the wild-type homozygous CC genotype had significantly longer horn than the heterozygous TC genotype (p = 0.0157). At chr3:133704023, the wild-type homozygous AA genotype had significantly longer horn than the heterozygous GA genotype (p = 0.0157). At chr3:133705250, the wild-type homozygous AA genotype had significantly longer horn than the heterozygous GA genotype (p = 0.0115). At chr3:133706677, the wild-type homozygous GG genotype had significantly longer horn than the heterozygous AG genotype (p = 0.0157). These results indicate that the wild-type homozygous genotype is associated with increased horn length, while the heterozygous genotype is associated with decreased horn length at these four loci. This suggests that the wild-type homozygous genotype may be the genotype associated with promoting horn growth. Notably, no individuals carrying the mutant homozygous genotype were detected for any of these four loci in our sampled population. We speculate that under certain conditions, these missing mutant homozygous genotypes may be lethal.

Figure 6

Figure 6. (A) Violin plot of horn length genotypes at chr3:1337303518, with p-value calculated by T-test (p = 0.0157). (B) Violin plot of horn length genotypes at chr3:1337304023, with p-value calculated by T-test (p = 0.0157). (C) Violin plot of horn length genotypes at chr3:1337305250, with p-value calculated by T-test (p = 0.0115). (D) Violin plot of horn length genotypes at chr3:1337306677, with p-value calculated by T-test (p = 0.0157).

4 Discussion

Keratin 1 (KRT1) plays a pivotal role in maintaining the structural integrity of the skin. In the upper epidermis, KRT1 forms heterodimers with KRT10, assembling into intermediate filaments that provide mechanical strength (48). Significantly, homozygous nonsense mutations in KRT1 (c.457C>T, p. Gln153; c.33C>G, p. Tyr11) have been linked to epidermolytic palmoplantar keratoderma (EPPK) with knuckle pads, resulting from nonsense-mediated mRNA decay and loss of KRT1 protein (13). Furthermore, intergenic variation near KRT1 correlates with the migratory capacity of human epidermal keratinocytes (HEKs) during wound healing (49), underscoring the critical function of KRT1 in skin development and homeostasis. Given that horns are skin-derived structures, with growth originating from the dermis and subcutaneous layers as demonstrated by Corner bud tissue transplantation studies in young calf and caprids (50), and considering the skin-specific high expression of KRT1 gene observed across species, we hypothesized that KRT1 may regulate horn development by influencing the growth and differentiation of keratinocytes within the skin.

Our RNA-seq analysis revealed significantly elevated KRT1 gene expression in the scurred group compared to the SHE group in sheep (p = 0.0024). We propose that this KRT1 gene overexpression in scurred individuals contributes to localized thickening of the keratinous horn sheath, potentially inhibiting the longitudinal growth of the underlying horn core. As structural components forming the cytoskeleton of keratinocytes, keratins directly impact epithelial differentiation and cornification–key processes in horn sheath formation (51). This aligns with prior research in humans and model organisms, where KRT1 gene mutations disrupt keratin filament assembly, leading to hyperkeratotic disorders like palmoplantar keratoderma (52), characterized by abnormal epidermal thickening analogous to the keratinization process in horn development. Our findings extend this mechanistic link to ovine horn morphogenesis, suggesting that KRT1 gene overexpression in scurred group accelerates keratinocyte differentiation and matrix hardening, thereby restricting horn size. The consistency of this model among different sheep breeds (the expression of KRT1 gene in hornless related breeds is higher than that in horned breeds) further suggests its evolutionary conservation in the regulation of horn development. Additionally, we observed significantly higher KRT1 gene expression in ram skin compared to ewe skin (p = 0.00691), consistent with the observation that males typically possess larger, more robust horn structures requiring reinforced frontal bones to accommodate growth and withstand increased mechanical stress during physical interactions (53).

Phylogenetic analysis of KRT1 protein sequences across sheep and 17 other species identified several horned-animal-specific amino acid residues (e.g., K312) located within the central α-helical rod domain. This domain is essential for intermediate filament formation (54), and mutations within the H1 and α-helical rod domains of K1/K10 are known to cause bullous congenital ichthyosiform erythroderma (55, 56), highlighting the functional importance of this region in keratinization. Other identified specific residues, often aromatic amino acids, reside in low-complexity, aromatic-rich segments (LARKS) within the head or tail domains. LARKS facilitate molecular interactions crucial for keratin intermediate filament assembly (57), suggesting these loci may represent potential regulatory points for horn development. Protein-protein interaction network analysis (STRING database) confirmed KRT10 as the primary interactor with KRT1. KRT1 (basic, Type II) pairs with KRT10 (acidic, Type I) to form heterodimers (48), which self-assemble into antiparallel, staggered tetramers, ultimately forming intermediate filaments through longitudinal and lateral interactions (48, 58). The KRT1-KRT10 heterodimer interface, stabilized in part by the Cys401K10 disulfide bond, forms a half-staggered antiparallel tetrameric complex (57). Notably, our analysis suggests that the horned-animal-specific residue K312 in KRT1 may form a key hydrogen bond with E262 in KRT10. Mutations at K312 to S/T are predicted to disrupt this bond, indicating that the conservation of K312 is likely critical for the stability of the KRT1/KRT10 complex and, consequently, for proper keratin structure formation.

Utilizing public data, this study identified 20 mutation loci with FST > 0.05 within the KRT1 gene and its flanking 1,000 bp regions, from which 32 potential functional loci were selected. These functional loci, prioritized either due to high FST values or specific selection criteria, are hypothesized to influence horn morphology and size by modulating gene expression or altering protein function. Subsequent linkage disequilibrium (LD) analysis delineated nine haplotype blocks encompassing 13 SNPs with FST > 0.05, suggesting these regions may harbor genetic units critically associated with horn type traits. Association analysis using self-measured average horn length data from Small Tail Han sheep revealed four loci (chr3:133705250, chr3:133704023, chr3:133706677, chr3:133703518) showing significant associations with horn length phenotypic variation. At these loci, the wild-type homozygous genotype was significantly associated with longer horns, while the heterozygous genotype corresponded to shorter horn lengths. This finding parallels observations by Zhao et al. (59). In their study on Gansu Alpine Fine-wool sheep, two SNPs (SNP1: C.-7G/C, SNP2: C.1500G/A) within the KRT71 gene, associated with wool traits, also showed that individuals homozygous for the wild-type allele had a significantly longer mean staple length (MSL) than heterozygous individuals (p < 0.05). Notably, our study work failed to detect individuals homozygous for the mutant allele at these significantly associated loci, leading to the hypothesis that this genotype may confer lethality under certain conditions. Furthermore, this study identified four loci exhibiting allele-specific expression (ASE). At these ASE loci, the reference allele (ref) count was significantly higher in the SHE group compared to the scurred group, while the alternative allele (alt) count was zero in the SHE group. This pattern strongly suggests that these loci may cooperatively regulate horn size development. Hereditary ichthyoses, known to be caused by dominant-negative mutations in keratin genes such as KRT1, KRT2, and KRT10, disrupt epidermal keratinization (60). Drawing an analogy to this established pathological mechanism, we speculate that the specific allelic variations identified in the KRT1 gene in this study (i.e., the FST > 0.05 loci, association loci, and ASE loci) may similarly interfere with epidermal keratinization or related processes, thereby influencing horn morphogenesis and size regulation.

This study provides valuable insights for sheep breeding programs targeting horn phenotypes. Future research should integrate single-cell RNA sequencing (scRNA-seq), epigenomic analyses, and histological examinations for multi-omics profiling. This approach will help identify specific cell types and locate critical regulatory elements of KRT1 gene influencing horn development. Functional validation through CRISPR-Cas9-mediated mutagenesis or in vitro assays will be essential to confirm the transcriptional regulatory effects of prioritized SNPs and elucidate the precise molecular mechanisms by which KRT1 gene governs horn morphogenesis in sheep.

5 Conclusions

This study demonstrates a significant association between the KRT1 gene and horn size in sheep. Expression analysis revealed a relatively higher level of KRT1 expression in the scurred sheep group. Cross-species investigations further confirmed the high tissue specificity of KRT1, with robust expression observed predominantly in the skin and horn tissue of multiple species, including sheep, cattle, pigs, and humans. In-depth analysis of the amino acid sequence encoded by the KRT1 gene and its protein structure identified 11 amino acid loci exhibiting specificity in horned animals. Crucially, an interaction was discovered between lysine at position 312 (K312) of the KRT1 protein and glutamate at position 262 (E262) of the KRT10 protein. This specific interaction loci may play a pivotal role in keratin filament network formation and horn development. Furthermore, this study integrated RNA-seq and whole-genome sequencing (WGS) datasets. Through subsequent analysis of SNP variations in the KRT1 gene, allele-specific expression (ASE), and population differentiation index (Fst), a series of potential functional regulatory loci and genetic variants significantly associated with horn length were successfully identified. These loci likely influence horn size and morphology (horned, polled, or scurred) by modulating gene function or expression levels. Based on these findings, we propose that these loci serve as potential molecular markers for horn type traits (including horn length and presence/absence) in sheep. In conclusion, this research provides novel insights into the biological function of the KRT1 gene in ovine horn development and establishes a crucial theoretical foundation for future molecular breeding strategies targeting horn phenotype in sheep.

Data availability statement

The original data presented in the study are publicly available. The Tibetan sheep RNA-seq data can be found in the NCBI repository under accession number PRJNA1003277. Publicly available datasets analyzed in this study are cited within the article and their sources are provided in the Methods section.

Ethics statement

This study was reviewed and approved by the Animal Welfare and Ethics Committee of the Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS-CAAS; Approval No. IAS-CAAS-AE-03; December 12, 2016). All procedures were strictly conducted in accordance with the institutional guidelines and regulations for the care and use of laboratory animals. Since the Tibetan sheep tissue samples used in this study were obtained through collaborative arrangements for research purposes, written informed consent was not required. The sampling protocol constituted an integral part of the research design and was explicitly approved by the Ethics Committee. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

QM: Writing – original draft, Data curation, Formal analysis, Writing – review & editing. HL: Writing – review & editing, Data curation, Formal analysis. GZ: Writing – review & editing, Data curation, Formal analysis. ZM: Writing – review & editing, Data curation. SZ: Writing – review & editing, Data curation. FL: Writing – review & editing, Data curation. ZP: Writing – review & editing, Methodology. JY: Funding acquisition, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by Proposal for National Natural Science Fund for Excellent Young Scientists Fund Program (Overseas; 2024-HY-01), Biological Breeding-National Science and Technology Major Project (2022ZD04017 and 2023ZD04051), and National Key R&D Program of China (2023YFF1001800 and 2022YFF1000100).

Conflict of interest

The authors declare that the research 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) declare that no Gen AI was 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.1703700/full#supplementary-material

References

1. Clark EL, Bush SJ, McCulloch MEB, Farquhar IL, Young R, Lefevre L, et al. A high resolution atlas of gene expression in the domestic sheep (Ovis aries). PLoS Genet. (2017) 13:e1006997. doi: 10.1371/journal.pgen.1006997

PubMed Abstract | Crossref Full Text | Google Scholar

2. Cappelli J, García AJ, Kotrba R, Gambín Pozo P, Landete-Castillejos T, Gallego L, et al. The bony horncore of the common eland (Taurotragus oryx): composition and mechanical properties of a spiral fighting structure. J Anat. (2018) 232:72–9. doi: 10.1111/joa.12708

PubMed Abstract | Crossref Full Text | Google Scholar

3. Preston BT, Stevenson IR, Pemberton JM, Coltman DW, Wilson K. Overt and covert competition in a promiscuous mammal: the importance of weaponry and testes size to male reproductive success. Proc Biol Sci. (2003) 270:633–40. doi: 10.1098/rspb.2002.2268

PubMed Abstract | Crossref Full Text | Google Scholar

4. Robinson MR, Kruuk LEB. Function of weaponry in females: the use of horns in intrasexual competition for resources in female Soay sheep. Biol Lett. (2007) 3:651–4. doi: 10.1098/rsbl.2007.0278

PubMed Abstract | Crossref Full Text | Google Scholar

5. Case study: prevalence of horns and bruising in feedlot cattle at slaughter. Prof Anim Sci. (2017) 33:135–9. doi: 10.15232/pas.2016-01551

Crossref Full Text | Google Scholar

6. Simon R, Drögemüller C, Lühken G. The complex and diverse genetic architecture of the absence of horns (polledness) in domestic ruminants, including goats and sheep. Genes. (2022) 13:832. doi: 10.3390/genes13050832

PubMed Abstract | Crossref Full Text | Google Scholar

7. Still Brooks KM, Hempstead MN, Anderson JL, Parsons RL, Sutherland MA, Plummer PJ, et al. Characterization of efficacy and animal safety across four caprine disbudding methodologies. Animals. (2021) 11:430. doi: 10.3390/ani11020430

PubMed Abstract | Crossref Full Text | Google Scholar

8. Graduate student literature review: role of pain mitigation on the welfare of dairy calves undergoing disbudding. J Dairy Sci. (2022) 105:6809–19. doi: 10.3168/jds.2021-21349

Crossref Full Text | Google Scholar

9. Chen L, Qiu Q, Jiang Y, Wang K, Lin Z, Li Z, et al. Large-scale ruminant genome sequencing provides insights into their evolution and distinct traits. Science. (2019) 364:eaav6202. doi: 10.1126/science.aav6202

PubMed Abstract | Crossref Full Text | Google Scholar

10. Luan Y, Wu S, Wang M, Pu Y, Zhao Q, Ma Y, et al. Identification of critical genes for ovine horn development based on transcriptome during the embryonic period. Biology. (2023) 12:591. doi: 10.3390/biology12040591

PubMed Abstract | Crossref Full Text | Google Scholar

11. Tian D, Han B, Li X, Liu D, Zhou B, Zhao C, et al. Genetic diversity and selection of Tibetan sheep breeds revealed by whole-genome resequencing. Anim Biosci. (2023) 36:991–1002. doi: 10.5713/ab.22.0432

PubMed Abstract | Crossref Full Text | Google Scholar

12. Guo T, Zhao H, Yuan C, Huang S, Zhou S, Lu Z, et al. Selective sweeps uncovering the genetic basis of horn and adaptability traits on fine-wool sheep in China. Front Genet. (2021) 12:604235. doi: 10.3389/fgene.2021.604235

PubMed Abstract | Crossref Full Text | Google Scholar

13. Mo R, Lin M, Lee M, Yan W, Wang H, Lin Z. Nonsense mutations in KRT1 caused recessive epidermolytic palmoplantar keratoderma with knuckle pads. J Eur Acad Dermatol Venereol. (2022) 36:1857–62. doi: 10.1111/jdv.18189

PubMed Abstract | Crossref Full Text | Google Scholar

14. Kim T, Kim S-C, Lee SE. Two cases of KRT1 mutation-associated epidermolytic ichthyosis without typical epidermolytic hyperkeratosis in the neonatal skin lesions. Pediatr Dermatol. (2023) 40:1149–51. doi: 10.1111/pde.15354

PubMed Abstract | Crossref Full Text | Google Scholar

15. Starr NJ, Khan MH, Edney MK, Trindade GF, Kern S, Pirkl A, et al. Elucidating the molecular landscape of the Stratum corneum. Proc Natl Acad Sci U S A. (2022) 119:e2114380119. doi: 10.1073/pnas.2114380119

PubMed Abstract | Crossref Full Text | Google Scholar

16. Li H, Du X, Li X, Feng P, Chu M, Jin Y, et al. Genetic diversity, tissue-specific expression, and functional analysis of the ATP7A gene in sheep. Front Genet. (2023) 14:1239979. doi: 10.3389/fgene.2023.1239979

PubMed Abstract | Crossref Full Text | Google Scholar

17. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | Crossref Full Text | Google Scholar

18. Shumate A, Wong B, Pertea G, Pertea M. Improved transcriptome assembly using a hybrid of long and short reads with StringTie. PLoS Comput Biol. (2022) 18:e1009730. doi: 10.1371/journal.pcbi.1009730

PubMed Abstract | Crossref Full Text | Google Scholar

19. Zhang G, Chu M, Yang H, Li H, Shi J, Feng P, et al. Expression, polymorphism, and potential functional sites of the BMPR1A gene in the sheep horn. Genes. (2024) 15:376. doi: 10.3390/genes15030376

PubMed Abstract | Crossref Full Text | Google Scholar

20. Teng J, Gao Y, Yin H, Bai Z, Liu S, Zeng H, et al. A compendium of genetic regulatory effects across pig tissues. Nat Genet. (2024) 56:112–23. doi: 10.1038/s41588-023-01585-7

PubMed Abstract | Crossref Full Text | Google Scholar

21. Liu S, Gao Y, Canela-Xandri O, Wang S, Yu Y, Cai W, et al. A multi-tissue atlas of regulatory variants in cattle. Nat Genet. (2022) 54:1438–47. doi: 10.1038/s41588-022-01153-5

PubMed Abstract | Crossref Full Text | Google Scholar

22. GTEx Consortium. Human genomics. The Genotype-Tissue Expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. (2015) 348:648–60. doi: 10.1126/science.1262110

PubMed Abstract | Crossref Full Text | Google Scholar

23. GTEx Consortium. The GTEx consortium atlas of genetic regulatory effects across human tissues. Science. (2020) 369:1318–30. doi: 10.1126/science.aaz1776

PubMed Abstract | Crossref Full Text | Google Scholar

24. Tamura K, Stecher G, Kumar S. MEGA11: molecular evolutionary genetics analysis version 11. Mol Biol Evol. (2021) 38:3022–7. doi: 10.1093/molbev/msab120

PubMed Abstract | Crossref Full Text | Google Scholar

25. CUDA ClustalW: an efficient parallel algorithm for progressive multiple sequence alignment on multi-GPUs. Comput Biol Chem. (2015) 58:62–8. doi: 10.1016/j.compbiolchem.2015.05.004

Crossref Full Text | Google Scholar

26. Ge T, Luo X, Wang Y, Sedlmair M, Cheng Z, Zhao Y, et al. Optimally ordered orthogonal neighbor joining trees for hierarchical cluster analysis. IEEE Trans Vis Comput Graph. (2024) 30:5034–46. doi: 10.1109/TVCG.2023.3284499

PubMed Abstract | Crossref Full Text | Google Scholar

27. Jumper J, Evans R, Pritzel A, Green T, Figurnov M, Ronneberger O, et al. Highly accurate protein structure prediction with AlphaFold. Nature. (2021) 596:583–9. doi: 10.1038/s41586-021-03819-2

PubMed Abstract | Crossref Full Text | Google Scholar

28. Szklarczyk D, Kirsch R, Koutrouli M, Nastou K, Mehryary F, Hachilif R, et al. The STRING database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. (2022) 51:D638–46. doi: 10.1093/nar/gkac1000

PubMed Abstract | Crossref Full Text | Google Scholar

29. Chen E, Pan E, Zhang S. Structure bioinformatics of six human integral transmembrane enzymes and their AlphaFold3 predicted water-soluble QTY analogs: insights into FACE1 and STEA4 binding mechanisms. Pharm Res. (2025) 42:291–305. doi: 10.1007/s11095-025-03822-6

PubMed Abstract | Crossref Full Text | Google Scholar

30. Pan Z, Li S, Liu Q, Wang Z, Zhou Z, Di R, et al. Whole-genome sequences of 89 Chinese sheep suggest role of RXFP2 in the development of unique horn phenotype as response to semi-feralization. Gigascience. (2018) 7:giy019. doi: 10.1093/gigascience/giy019

PubMed Abstract | Crossref Full Text | Google Scholar

31. Li C, Li M, Li X, Ni W, Xu Y, Yao R, et al. Whole-genome resequencing reveals loci associated with thoracic vertebrae number in sheep. Front Genet. (2019) 10:674. doi: 10.3389/fgene.2019.00674

PubMed Abstract | Crossref Full Text | Google Scholar

32. Allais-Bonnet A, Hintermann A, Deloche M-C, Cornette R, Bardou P, Naval-Sanchez M, et al. Analysis of polycerate mutants reveals the evolutionary co-option of HOXD1 for horn patterning in bovidae. Mol Biol Evol. (2021) 38:2260–72. doi: 10.1093/molbev/msab021

PubMed Abstract | Crossref Full Text | Google Scholar

33. Qiao G, Xu P, Guo T, Wu Y, Lu X, Zhang Q, et al. Genetic basis of Dorper sheep (Ovis aries) revealed by long-read de novo genome assembly. Front Genet. (2022) 13:846449. doi: 10.3389/fgene.2022.846449

PubMed Abstract | Crossref Full Text | Google Scholar

34. Luo R, Dai X, Zhang L, Li G, Zheng Z. Genome-wide DNA methylation patterns of muscle and tail-fat in DairyMeade sheep and Mongolian sheep. Animals. (2022) 12:1399. doi: 10.3390/ani12111399

PubMed Abstract | Crossref Full Text | Google Scholar

35. Schultz DT, Haddock SHD, Bredeson JV, Green RE, Simakov O, Rokhsar DS. Ancient gene linkages support ctenophores as sister to other animals. Nature. (2023) 618:110–7. doi: 10.1038/s41586-023-05936-6

PubMed Abstract | Crossref Full Text | Google Scholar

36. Posbergh CJ, Staiger EA, Huson HJ. A stop-gain mutation within MLPH Is responsible for the lilac dilution observed in jacob sheep. Genes. (2020) 11:618. doi: 10.3390/genes11060618

PubMed Abstract | Crossref Full Text | Google Scholar

37. Zhao Y, Zhang X, Li F, Zhang D, Zhang Y, Li X, et al. Whole genome sequencing analysis to identify candidate genes associated with the rib eye muscle area in Hu sheep. Front Genet. (2022) 13:824742. doi: 10.3389/fgene.2022.824742

PubMed Abstract | Crossref Full Text | Google Scholar

38. Guo Y, Liang J, Lv C, Wang Y, Wu G, Ding X, et al. Sequencing reveals population structure and selection signatures for reproductive traits in Yunnan semi-fine wool sheep (Ovis aries). Front Genet. (2022) 13:812753. doi: 10.3389/fgene.2022.812753

PubMed Abstract | Crossref Full Text | Google Scholar

39. Li X, Yang J, Shen M, Xie X-L, Liu G-J, Xu Y-X, et al. Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits. Nat Commun. (2020) 11:2815. doi: 10.1038/s41467-020-16485-1

PubMed Abstract | Crossref Full Text | Google Scholar

40. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics. (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | Crossref Full Text | Google Scholar

41. Li H, Durbin R. Fast and accurate short read alignment with burrows–wheeler transform. Bioinformatics. (2009) 25:1754–60. doi: 10.1093/bioinformatics/btp324

PubMed Abstract | Crossref Full Text | Google Scholar

42. Dong S-S, He W-M, Ji J-J, Zhang C, Guo Y, Yang T-L. LDBlockShow: a fast and convenient tool for visualizing linkage disequilibrium and haplotype blocks based on variant call format files. Brief Bioinform. (2021) 22:bbaa227. doi: 10.1093/bib/bbaa227

PubMed Abstract | Crossref Full Text | Google Scholar

43. Li H, Durbin R. Fast and accurate long-read alignment with burrows–wheeler transform. Bioinformatics. (2010) 26:589–95. doi: 10.1093/bioinformatics/btp698

PubMed Abstract | Crossref Full Text | Google Scholar

44. Van der Auwera GA, Carneiro MO, Hartl C, Poplin R, del Angel G, Levy-Moonshine A, et al. From FastQ data to high-confidence variant calls: the genome analysis toolkit best practices pipeline. Curr Protoc Bioinform. (2013) 43:11.10.1–33. doi: 10.1002/0471250953.bi1110s43

PubMed Abstract | Crossref Full Text | Google Scholar

45. 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

PubMed Abstract | Crossref Full Text | Google Scholar

46. Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly. (2012) 6:80–92. doi: 10.4161/fly.19695

Crossref Full Text | Google Scholar

47. 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

PubMed Abstract | Crossref Full Text | Google Scholar

48. Zhang X, Yin M, Zhang L. Keratin 6, 16 and 17—critical barrier alarmin molecules in skin wounds and psoriasis. Cells. (2019) 8:807. doi: 10.3390/cells8080807

PubMed Abstract | Crossref Full Text | Google Scholar

49. Tao H, Berno AJ, Cox DR, Frazer KA. In vitro human keratinocyte migration rates are associated with SNPs in the KRT1 interval. PLoS ONE. (2007) 2:e697. doi: 10.1371/journal.pone.0000697

PubMed Abstract | Crossref Full Text | Google Scholar

50. Aldersey JE, Sonstegard TS, Williams JL, Bottema CDK. Understanding the effects of the bovine POLLED variants. Anim Genet. (2020) 51:166–76. doi: 10.1111/age.12915

PubMed Abstract | Crossref Full Text | Google Scholar

51. Jaeger K, Sukseree S, Zhong S, Phinney BS, Mlitz V, Buchberger M, et al. Cornification of nail keratinocytes requires autophagy for bulk degradation of intracellular proteins while sparing components of the cytoskeleton. Apoptosis. (2019) 24:62–73. doi: 10.1007/s10495-018-1505-4

PubMed Abstract | Crossref Full Text | Google Scholar

52. Shchagina O, Fedotov V, Markova T, Shatokhina O, Ryzhkova O, Fedotova T, et al. Palmoplantar keratoderma: a molecular genetic analysis of family cases. Int J Mol Sci. (2022) 23:9576. doi: 10.3390/ijms23179576

PubMed Abstract | Crossref Full Text | Google Scholar

53. Gündemir O, Szara T. Morphological patterns of the European bison (Bison bonasus) skull. Sci Rep. (2025) 15:1418. doi: 10.1038/s41598-025-85654-3

PubMed Abstract | Crossref Full Text | Google Scholar

54. Vermeire P-J, Stalmans G, Lilina AV, Fiala J, Novak P, Herrmann H, et al. Molecular interactions driving intermediate filament assembly. Cells. (2021) 10:2457. doi: 10.3390/cells10092457

PubMed Abstract | Crossref Full Text | Google Scholar

55. Suga Y, Duncan KO, Heald PW, Roop DR. A novel helix termination mutation in keratin 10 in annular epidermolytic ichthyosis, a variant of bullous congenital ichthyosiform erythroderma. J Invest Dermatol. (1998) 111:1220–3. doi: 10.1046/j.1523-1747.1998.00451.x

PubMed Abstract | Crossref Full Text | Google Scholar

56. Ishida-Yamamoto A, Takahashi H, Iizuka H. Lessons from disorders of epidermal differentiation-associated keratins. Histol Histopathol. (2002) 17:331–8. doi: 10.14670/HH-17.331

PubMed Abstract | Crossref Full Text | Google Scholar

57. Bunick CG, Milstone LM. The x-ray crystal structure of the keratin 1–keratin 10 helix 2B heterodimer reveals molecular surface properties and biochemical insights into human skin disease. J Invest Dermatol. (2017) 137:142–50. doi: 10.1016/j.jid.2016.08.018

PubMed Abstract | Crossref Full Text | Google Scholar

58. Lomakin IB, Hinbest AJ, Ho M, Eldirany SA, Bunick CG. Crystal structure of keratin 1/10 (C401A) 2B heterodimer demonstrates a proclivity for the C-terminus of helix 2B to form higher order molecular contacts. Yale J Biol Med. (2020) 93:3–17.

Google Scholar

59. Zhao F, He Z, Sun H, Wang J, Liu X, Hao Z, et al. Polymorphism of keratin gene KRT71 and its relationship with wool properties in Gansu alpine fine-wool sheep. Animals. (2025) 15:2028. doi: 10.3390/ani15142028

PubMed Abstract | Crossref Full Text | Google Scholar

60. Hotz A, Oji V, Bourrat E, Jonca N, Mazereeuw-Hautier J, Betz RC, et al. Expanding the clinical and genetic spectrum of KRT1, KRT2 and KRT10 mutations in keratinopathic ichthyosis. Acta Derm Venereol. (2016) 96:473–8. doi: 10.2340/00015555-2299

PubMed Abstract | Crossref Full Text | Google Scholar