A novel deletion in KRT75L4 mediates the frizzle trait in a Chinese indigenous chicken

Background Highly diversified in morphology and structure, feathers have evolved into various forms. Frizzle feathers, which result from a developmental defect of the feather, are observed in several domestic chicken breeds. The frizzle phenotype is consistent with incomplete dominance of a major gene, but the molecular mechanisms that underlie this phenotype remain obscure. Kirin, a Chinese indigenous chicken breed that originated in the Guangdong province, is famous for its frizzle feathers. The KRT75 gene is considered as the dominant gene responsible for the frizzle trait in several chicken breeds, but this is not the case in the Kirin breed. Thus, the objective of our study was to investigate the genomic region and mutation responsible for this phenotype in this particular breed. Results A resource population was produced by crossing Kirin and Huaixiang chickens to produce F1 and F2 generations. DNA samples from 75 frizzle feather and normal feather individuals were sequenced with double-digest genotyping by sequencing (dd-GBS). After the detection of 525,561 high-quality variants, a genome-wide association analysis was carried out and the gene responsible for the frizzle phenotype was localized within the type II α-keratin cluster on chromosome 33. Sanger sequencing was used to screen for mutations in the exons of five genes of this type II α-keratin cluster. A 15-bp deletion in exon 3 of KRT75L4 that showed complete segregation with the frizzle phenotype was detected within the F2 population. Transcriptome sequencing demonstrated that KRT75L4 was expressed but that the transcript was shorter in Kirin than in Huaixiang chickens. In addition, by using Sanger sequencing, we were able to confirm that the deletion was in complete linkage with frizzle feathers. Conclusions A deletion in the KRT75L4 gene is responsible for the frizzle feather phenotype in the Kirin chicken. The identification of this mutation, which causes a developmental defect of avian integument appendages, will improve our understanding of the mechanisms that are involved in feather formation. Electronic supplementary material The online version of this article (10.1186/s12711-018-0441-7) contains supplementary material, which is available to authorized users.


Background
Birds display a high degree of diversity in feather structure. In chickens, the frizzle trait is a varietal characteristic, and birds that carry this trait have all their contour feathers curling outward and upward [1,2]. The rachises of frizzle and normal feathers are different i.e., frizzle feathers turn outward relative to the skin, whereas normal body and flight feathers have a dorsal orientation. In frizzle chickens, rectrices and remiges are less affected but have an irregular appearance, and display other modifications, such as thickening of the barbs and barbules, alteration of the hooklets and in some cases, other structural abnormalities [3,4]. These modifications of the feather structure reduce the insulating effect of feathers [5][6][7].
To date, the frizzle mutation has been reported to occur in a single autosomal gene, denoted as the F gene (see [5]), and shows incomplete dominance inheritance [4,5]. Although the curling of feathers of the frizzle phenotype is visible in heterozygous individuals, it is even more pronounced in homozygous individuals, but frizzle feathers are more ornamental (the curve is less severe and feathers look like daisy petals) in heterozygous than in homozygous chickens. Because of their enhanced heat dissipation under tropical conditions, both homozygous and heterozygous frizzle chickens have a relatively higher meat and egg production than wild type individuals [8,9].
Kirin chicken (KRC) is a Chinese indigenous breed of chicken that is noted for its frizzle feathers. The KRC breed originates from the Guangdong province, and is raised exclusively in this region with a warm and humid subtropical climate. Adult KRC have a distinct disorientation of feathers ( Fig. 1a-d). In general, the frizzle phenotype is not observed in the down of frizzle dwarf chicken at hatch but becomes more obvious in the second-generation bilateral feathers [10]. In contrast, in KRC, Tao et al. [11] observed the presence of frizzle downy feathers in chicks at hatch (Fig. 1e, f ), which suggests a different mechanism for the frizzle trait in this breed. Ng et al. [10] demonstrated that a deletion at the junction between exon 5 and intron 5 of the KRT75 gene causes the frizzle feather phenotype in frizzle dwarf chicken, but this deletion is absent in KRC, which suggests that another gene is involved in the formation of feathers [10][11][12] and that the molecular mechanism responsible for frizzle feathers in KRC is different.
To dissect the genetic basis of the frizzle phenotype in KRC, we crossed homozygous frizzle KRC with normal feathered Huaixiang chickens (HX). A double-digest genotyping by sequencing (dd-GBS) based [13] genome-wide association analysis on individuals from the F 2 generation, allowed us to identify a candidate region for the F gene, in KRC. Then, we used Sanger sequencing to ascertain the causative mutation in several candidate genes. We also compared the transcriptome profiles of KRC and HX. Finally, we validated our finding in other breeds, e.g. Houdan (H), Luhua (L) and HX. Fig. 1 Comparison between frizzle and normal phenotypes. a-d Adult Kirin (frizzle feathers) and Huaixiang (normal feathers) chickens. Adult frizzle feathers curve away from the body, which gives a highly divergent appearance from normal feathers. e One-day-old Kirin and f Huaixiang chickens. The downy feathers appear frizzled, which contrasts with feathers from previously reported breeds

Animals and dd-GBS sequencing
A three-generation resource population was produced by crossing purebred KRC with frizzle feathers and HX with normal feathers. The F 1 generation was derived by mating one KRC male with 10 HX females (KRC × HX), and by mating one HX male with 10 KRC females (HX × KRC). Three male and eight female F 1 individuals from KRC × HX and three male and eight female F 1 individuals from HX×KRC were mated, respectively, to produce the F 2 generation. We conducted a χ 2 -test to study inheritance patterns of the frizzle trait by using the R package software [14].
Eighteen F 1 and 57 F 2 homozygous frizzle feather and normal feather individuals were sampled from 12 families [see Additional file 1: Table S1]. Genomic DNA was extracted following the instructions of the AxyPrep ™ DNA Gel Extraction kit, double-digested with different restriction enzyme pairs and sequenced. We found that digestion with the enzyme combination, Pst I and Msp I, was appropriate and highly consistent with in silico digestion-site prediction. Thus, these two restriction enzymes were used to digest total DNA. Adapters marked with specific barcodes were ligated to the restriction-digested overhanging sequences by T4 ligase, following the dd-GBS protocol [13]. DNA samples from the above 75 individuals, marked by different barcode adapters, were mixed and purified according to the manufacturer's instructions. DNA fragments that carried ligated adapters were amplified with primers that had complementary sequences for each adapter. The constructed 2 × 150 bp libraries were sequenced by the Illumina HiSeq 2000 sequencing system in Shanghai Personal Biotechnology Co., Ltd.
All raw sequences are deposited at NCBI under the BioProject registration number PRJNA445355. Removal of adapter sequences from the reads was done by using AdapterRemoval v2 [15] and the trimmed reads were mapped to the reference genome of Gallus gallus (gal-Gal5) by the Burrows-Wheeler Alignment (BWA) tool v0.7.12 [16] in order to identify variants. High-confidence single nucleotide polymorphisms (SNPs) were called by using the GATK software [17], and the outputs were transformed into PED and MAP files format, which were the standard input files of PLINK-1.9 [18].

Genome-wide association analysis
The dd-GBS based SNPs were filtered out based on the following criteria: (1) minor allele frequency lower than 0.05; (2) SNP missing rate across all individuals higher than 0.2; (3) P value for Hardy-Weinberg equilibrium less than 1 × 10e −40 ; and (4) Mendelian errors based on parental genotypes. GWAS on 57 F 2 samples were performed by using the case-control model in PLINK.

Sanger sequencing
In order to detect all possible causative mutations in five keratin genes (KRT5, KRT6A, KRT75, KRT75L2 and KRT75L4), we used exon sequencing. Primers for these candidate genes were designed by using primer premier 5.0 and are provided in Additional file 2: Table S2. DNA obtained from blood samples of five homozygous normal and five frizzle feather individuals was pooled, PCR-amplified and sequenced, respectively, with three replicates for each pool.
To validate that the mutation was in exon 3 of the KRT75L4 gene, blood was sampled from individuals that were randomly selected from KRC, HX, H and L chicken. Heterozygous frizzle samples (F) were obtained from crosses between KRC and other breeds with normal feather breeds. DNA was extracted and pooled from 49 H, 34

Transcriptome analysis of frizzle and normal Kirin chicken
We performed the transcriptome analysis of follicle tissues from three KRC and three HX individuals at the embryonic age of E13. Under an RNase-free laminar flow cabinet, feather follicles with skin were removed surgically from the back of each individual, sprayed with RNAsafety, and immediately frozen in liquid nitrogen. Total RNA was extracted with Trizol reagent, and mRNA sequencing was carried out using poly (A)-enriched mRNA. Each of these mRNA samples was converted first into barcoded RNA-seq libraries, and then sequenced on an Illumina HiSeq 2000 sequencer.
The FASTQ sequence of each library was submitted to NCBI (PRJNA445349). Using the Tophat2 [19] program with default parameters, stock-specific transcripts were mapped to the chicken reference genome. The normalized gene expression levels were measured in fragments per kb of exon per million fragments mapped (FPKM) using Cufflinks. Differentially expressed genes (DEG) were identified using the DESeq package [20] based on a P value less than 0.05 and |log2(K/H)| higher or equal to 1, where K and H denote the expression level of each gene in the KRC and HX groups, respectively.

Resource population from crossbreeding between KRC and HX
Reciprocal matings between homozygous KRC and HX resulted in a heterozygous F 1 , and, as previously reported [3][4][5], the frizzle phenotype was incompletely penetrant in this population. Penetrance was similar regardless of whether the F gene was inherited from a cockerel or a hen, thus ruling out sex-linked inheritance. In order to determine the gene underlying the frizzle trait, six males and 16 females were selected from the F 1 generation to produce 486 F 2 offspring. Of these, 113 individuals were homozygous frizzle, 241 were heterozygous frizzle and 132 were homozygous normal feathered. The χ 2 -test (df = 2) showed that the frizzle phenotype is consistent with Mendelian inheritance (P = 0.22) and incompletely dominant in KRC.

QTL mapping by genome-wide association analysis
dd-GBS sequencing applied to 75 F 1 and F 2 individuals generated 1464 million reads for an average of 3× sequencing depth. After filtering, these reads were aligned to the reference genome, generating about 100,000 effective fragments, which is equivalent to 30 Mb and represents 2.8% of the entire genome with an average depth of 54×. A total of 1,948,658 raw variants were identified, loaded into PLINK 1.9 [18], and filtered according to the criteria presented in the Methods section, which yielded 525,561 variants. The distribution of these variants across the chromosomes was relatively even [see Additional file 3: Figure S1]. Eighteen of the 75 F 1 individuals were used as references to remove variants with Mendelian errors.
We applied a case-control model to calculate the association of variants with the frizzle phenotype across the whole genome (see Fig. 2). SNPs that showed a significant association and their nearest gene are listed in Table 1. A majority of the most significant SNPs were located within a region of 130 kb (P = 9.46 × 10 −26 ), which corresponds to intron 5 of the keratin 75-like 4 gene (KRT75L4). The other most significant SNPs were located in the vicinity of the KRT75L2, KRT5, KRT7 and KRT6A genes, which encompass a region between 1.19 and 1.52 Mb on chromosome 33 and are putative candidates for the F locus. In summary, the candidate QTL overlaps (1.19-1.39 Mb) with the cluster of α-keratin genes on chromosome 33 [21]. To our knowledge, the other detected genes, such as FAIM2, AQP5 and GALNT6 are not involved in feather formation.

Identification of a mutation in KRT75L4 in the KRC
The candidate region identified by GWAS lies within a genomic interval on chromosome 33, which covers the cluster of type II α-keratin genes. Although in frizzle dwarf chicken, the causative mutation for the frizzle phenotype is localized in the KRT75 gene [10], previous results indicate that this is not the case in KRC [11]. To detect the putative causative mutations in KRC, we PCR-amplified and sequenced all the exons of five keratin genes (KRT5, KRT6A, KRT75, KRT75L2, and KRT75L4) located within the region detected by GWAS. The KRT75L4 gene has three isoforms (ENSGALT00000086819, ENSGALT00000090244 and ENSGALT00000047324) and the transcript is associated with two gene entries (ENSGALG00000044875 and ENSGALG00000043689). We analyzed the DNA extracted from 10 homozygous F 2 individuals, i.e. five with normal and five with frizzle feathers and detected 90 synonymous substitutions, five missense substitutions, two deletions and three insertions (indels). Among these five indels, only one was located within the coding sequence of KRT75L4, i.e. a 15-bp deletion (Fig. 3), which can alter the sequence of the corresponding protein. This deletion showed complete segregation with the frizzle phenotype in the F 2 samples that were sequenced; and thus, was named F deletion.
All the other variants were excluded because at least one wild type chicken had the same genotype as that in the chickens with frizzle feathers.

Whole-transcriptome sequencing
Whole-transcriptome sequencing produced an average of 49 million reads for each sample. After the trimming and filtering steps, 37.4 to 53.6 million high-quality sequence reads from each of the individuals remained and were aligned to the reference genome. On average, 88% of the reads were successfully aligned. Of the 94 significant differentially expressed genes (DEG) based on a |log 2 (K/H)| higher than 1 and a P less than 0.05 [see Additional file 4: Table S4], 50 DEG were up-regulated and 44 were down-regulated in KRC. However, none of these DEG were located in the 1.19-1.52 Mb region detected by GWAS on chromosome 33 or in the type II α-keratin cluster.
The level of expression of KRT75L4 was similar in the KRC and HX groups. However, the sequence of the KRT75L4 transcript (displayed by the integrative genomics viewer (IGV) [22,23]) differed between the KRC and HX groups (Fig. 4). Due to the F deletion in one exon (exon 3 of ENSGALT00000090244.1 or exon 4 of ENSGALT00000086819.1), the KRT75L4 transcript was comparatively shorter in KRC than in HX, leading to the loss of Val-Leu-Pro-Ala-Ser in the translation product, which thus may have a different function compared to the normal protein. The transcripts of the other keratin genes on chromosome 33 were similar for both KRC and HX, which is concordant with the results from the analysis of the KRC's genome sequence. Thus, our findings indicate that KRT75L4 is probably the F gene.

Validation in other breeds
To validate furthermore the association detected between KRT75L4 and the frizzle feather phenotype, we PCRamplified and sequenced the F deletion and its flanking sequences located in KRT75L4 exon 3 from DNA obtained from other breeds of chickens. The F deletion was not found in HX, H and L chickens, which are breeds with normal feathers (Table 2). These results suggest that even if the F deletion is not the causative variant, it is strongly linked to the phenotype of frizzle feathers.

Discussion
Heat stress can be a major concern for chickens that are raised in warm climates [24,25]. A reduction in feather coverage increases heat dissipation and body heat irradiation, and such adaptation to high environmental temperatures has been well-documented [26][27][28]. Adomako et al. [8] reported favorable effects of the F locus on production traits, such as egg weight, egg mass, body weight and productivity index. Thus, identification of the F gene is important for breeding programs that are designed to enhance the productive performance of chickens maintained in hot humid climates.
Chr33: 1294968- Feathers consist mainly of two types of keratin proteins: α-and β-keratins, which are encoded by multigene families. β-keratins are found only in reptiles and birds, whereas α-keratins exist in all vertebrates [29]. Prum and Williamson [30] attributed the alteration of feather shape to mutations in keratin genes. Cellular and biochemical evidence indicates that α-keratin may have a key role in the early formation of rachides, barbs, and barbules [31]. In addition, the human and mouse homologs of KRT75, KRT6A and KRT6B are expressed in distinct parts of the hair follicle [32][33][34]. The morphological and structural diversity of feathers probably results from different combinations of α-and β-keratin genes in different intrafeather tissues [21].
Analyses of mouse and human cells showed that some α-keratins are involved in intracellular signaling pathways, and that mutations in the corresponding genes may interfere in the construction of the cell cytoskeleton [35,36]. Although type II keratin genes are located within a single cluster on human chromosome 12 and mouse chromosome 15 [37], studies on avian keratin genes have been challenging because of their conservation and duplications. In spite of considerable efforts to establish a correspondence between individual α-keratin chicken genes and their homologs in mammals [21], they are not included in the most recent chicken genome assembly and annotation release. Ng et al. [10] reported that the loss of the authentic splice site at the exon5/intron 5 junction of KRT75 resulted in a 69-bp in-frame deletion within its coding region and played a key role in frizzle feather development. The gene referred to as KRT75 [10] is actually KRT6A according to RefSeq. Furthermore, KRT75L4 lies within the noncoding region of KRT6A, along with KRT75L1, KRT75L2, KRT5 and other predicted loci. The official full name of KRT75L4 is keratin, type II cytoskeletal 75-like 4, and encodes an uncharacterized protein, with an annotation score of 1 [38]. Our results of transcriptome analysis, combined with other evidence, show that the gene is expressed in chicken tissues, although its function is unknown.