Study on LOC426217 as a candidate gene for beak deformity in chicken

Background The beak deformity (crossed beaks) was found in some indigenous chickens of China, such as Beijing-You (BJY), Qingyuan Partridge, and Huxu Chickens. Birds with deformed beaks have reduced feed intake and drinking, impeded growth rate, and poor production performance. Beak deformity reduces the economy of poultry industry and affects animal welfare as well. The genetic basis of this malformation remains incompletely understood. LOC426217, also named claw keratin-like, was the most up-regulated gene in the deformed beaks from a previous digital gene expression (DGE) analysis and was selected as an important candidate gene for further analysis. Results In the present study, quantitative real-time PCR (qRT-PCR) was firstly performed to determine the expression pattern of LOC426217 gene in deformed and normal beaks to verify the DGE results. Tissue-specific expression profile of this gene in 14 tissues was also determined using qRT-PCR. The LOC426217 was amplified from the genomic DNA of 171 deformed and 164 normal beaks, and sequenced to detect the single nucleotide polymorphisms (SNPs). The results showed that LOC426217 was significantly high-expressed in the deformed beaks, which was in good agreement with the DGE results. This gene was specifically high-expressed in beaks than other tissues. Eight SNPs were detected in LOC426217: -62G > T, 24 T > C, 36G > C, 192A > T, 204C > T, 222 T > C, 285G > T, and 363 T > C. Genotype frequency of G-62 T, T24C, G36C, T222C, and T363C loci was significant different between deformed and normal beaks. Haplotype analysis revealed one block with SNPs T24C and G36C, and one block with SNPs A192T, C204T, T222C, and G285T in normal birds, while the block with SNPs G36C and A192T in deformed ones. Conclusions It was concluded from these results that the over-expression of LOC426217 in the beak maybe related to the malformation. The polymorphisms of LOC426217 gene were associated with the beak deformity trait where the SNPs of G-62 T, T24C, G36C, T222C, and T363C loci maybe used as markers. The specific haplotype block in deformed birds may be a potential linkage marker for this trait. Electronic supplementary material The online version of this article (doi:10.1186/s12863-016-0353-x) contains supplementary material, which is available to authorized users.


Background
The beak is an external structure of birds, consisting of the upper and lower mandibles covered with a thin keratinized layer of epidermis [1]. It is used for many important activities such as feeding, drinking, fighting, and preening. In addition to striking morphological differences between species, beak deformities of different forms (noticeably elongated, crossed, bent at right angles) have been documented in many wild birds [2][3][4][5][6][7]. Frequencies of 1 % to 3 % of beak deformity (normally a crossed beak) were found in various indigenous chickens of China, such as Beijing-You (BJY) (studied here), Silkies, Qingyuan Partridge, and Huxu Chickens. Chickens with deformed beaks have reduced feed intake and growth rate. Therefore, beak deformity represents an economic as well as an animal welfare problem in poultry industry. According to our observations in a BJY population, in the absence of known environmental factors contributing to the malformation, birds with deformed beaks present consistently in each generation and cannot be eliminated from a population simply on the basis of the phenotype. This indicated the genetic effects underlying this trait. Studies have been performed to identify the teratogenic genes or molecular genetic background of beak deformity. Previously recognized genetic factors associated with beak deformity include some knwon genes such as fibroblast growth factor 8 (FGF8) [8], bone morphogenetic protein 4 (BMP4) [9][10][11], calmodulin (CaM) [12], and ALX homeobox 1 (ALX1) [13]. The over-expression of homeobox A1 (HOXA1) and homeobox D3 (HOXD3) may result in beak deformity in chicks [14].
Sets of differently expressed genes in the deformed and normal beaks have been detected using digital gene expression (DGE) profile analysis based on highsequencing technology. Of these genes, LOC426217, also known as claw keratin-like gene, was the most upregulated in the deformed beak (log2-Ratio (deformed/ normal) = 10.91) [15]. Located on GGA 25, LOC426217 is a member of the keratin family, containing 417 base pairs with only one exon (Fig. 1). Keratin is crucial for maintaining normal cell morphology involved in the cytoskeleton remodeling keratin filaments and cytoskeletal signaling pathways. Change of its structure results in dysmorphic cells [16]. The cytoskeleton is a complex of intracellular proteins that contribute to shape, support, and movement of cells [17]. Up to now, less study was reported about this gene in chickens. Although highlighted in the DGE analysis, further study of this gene is still needed to verify its roles in beak malformation.
In the present study, qRT-PCR was used to detect the relative expression of LOC426217 gene in the deformed and normal beaks, to verify the results of DGE profiling. Tissue expression profile of this gene was also determined in 14 tissues of the birds. Eventually, LOC426217 was amplified and sequenced to seeking the SNPs and haplotypes related with the beak malformation.

Animals and samples collection
The Institutional Animal Care and Use Committee at Institute of Animal Science, Chinese Academy of Agricultural Sciences (IAS, CAAS) approved all procedures involving the use of animals. All efforts were made to minimize the suffering of animals following the animal care guidelines [18]. The animals used in this study came from a pure-line stock of a local breed (Beijing-You) kept by IAS, CAAS (Beijing, China). They were incubated contemporaneously and housed under the same conditions.
The lower mandibles of the beaks were collected from 18 BJY chickens of 56 days of age: 9 with crossed beaks and 9 with normal beaks. Total RNA of the lower mandibles of the crossed and normal beaks above was collected for the verification of DGE profiling results using quantitative real-time PCR (qRT-PCR).
Three normal birds of 56 days of age were killed by stunning and exsanguination. Tissues samples including bursa of fabricius, beak, brain, breast, feather, heart, kidney, thigh, liver, lung, skin, small intestine, stomach, and testicle (50-100 mg) were rapidly collected and snapfroze in liquid nitrogen and storage at -80°C. The RNA of these samples was used to determine the tissue expression profile of LOC426217.
Blood samples were collected from the brachial vein by venipuncture from 171 beak-deformed birds (deformed) and 164 normal ones (control). Based on the case-control study design, we selected these birds according to the phenotype of the birds without family structure. The beak-deformity birds were collected from two generations. The normal birds were selected randomly from the same generation. DNA was isolated from the blood samples and stored at -20°C for the detection of SNPs located in LOC426217 gene.

DNA and RNA extraction and reverse transcription (RT)
Genomic DNA (gDNA) was extracted from blood samples using phenol-chloroform. Total RNA was isolated at 4°C using the Trizol reagent (Invitrogen, USA). Any residual gDNA and protein were removed with Dnase I   To validate the DGE results, qRT-PCR was performed to determine the expression of LOC426217 in the 9 deformed and 9 normal beaks using the ABI 7500 Realtime Detection System (Applied Biosystems, USA) and TaKaRa DRR018A reagents. Each 20 μL PCR mixture contained10 μL of SYBR Premix Ex Taq™ II, 0.8 μL (10 pM) of each primer (Table 2), 0.4 μL of ROX Reference Dye II (50×), 2 μL of cDNA (100 ng) and 6 μL of ddH 2 O. After an initial denaturing for 30 s at 95°C, there were 40 cycles of amplification (95°C for 5 s and 60°C for 32 s), followed by thermal denaturing to generate melting curves to verify amplification specificity. βactin was amplified in the same plates as endogenous control. Samples were assayed in triplicate for standard curves. PCR efficiency of the LOC426217 gene and βactin was consistent. cDNA from normal beaks served as a standard control for tissue-specific expression profile study. The amplification efficiency of transcripts of interest and the internal standard (β-actin) were consistent. Dissociation curves verified that amplification was specific.

SNPs filter and genotyping
PCR amplification product of LOC426217 gene was then directly sequenced by BGI company (Beijing, China)   Table 3 Genotypes of the birds used for qRT-PCR analysis for validation of Digital Gene Expression results

Trait
No. SNP Loci G-62 T  T24C  G36C  A192T  C204T  T222C  G285T  T363C   Deformed  1  TT  TC  GC  AT  CT  TT  GT  CC   2 T Normal  10  TT  TC  GC  AT  CC  TC  GT  TT   11  TT  CC  CC  AA  TT  TT  GG  TT   12  TT  TC  GC  AA  CT  TC  GG  TT   13  TT  TC  GC  AA  TT  TT  GG  TT   14  TT  TC  GC  AT  CT  TT  GT  CC   15  TT  TC  GC  AT  CT  TT  GT  CC   16  TT  TC  GC  AA  CT  TC  GG  TT   17  TT  CC  CC  AA  CT  TC  GG  TT   18  TT  TC  GC  AT  CC  TC  GT  TT using Sanger sequencing methods [19]. The SNPs and amino acids were determined and filtered by the software DNAStar (Version 5.01).

Statistical analysis
The relative abundance of transcripts was calculated from 2 −ΔΔCT [20]. All data presented graphically are means ± SEM. The significance level was P < 0.05 or P < 0.01. Student's t-tests were used to evaluate the relative expression differences of LOC426217 between the RNA samples of deformed and normal beaks. The ANOVA procedure of SAS 8.0 was used to assess the differences expression of LOC426217 in all the tissues. All SNPs were checked for Hardy-Weinberg Equilibrium (HWE) in both groups (P > 0.05 means equilibrium). Allele frequency, genotype frequency, and polymorphism information content (PIC, PIC < 0.25: low polymorphism; 0.25 < PIC < 0.5: moderate polymorphism; PIC > 0.5: high polymorphism) were calculated by PopGene (Version 1.31). Chi-square tests were used to evaluate the genotype frequency differences between deformed and normal beaks. Benjamini & Hochberg method was used for the Bonferroni correction [21]. Linkage Disequilibrium (LD) pattern for the SNPs genotyped was plotted using Haploview (Version 4.2). The sliding window method was used to generate different haplotypes between two groups [22,23].

Verification for LOC426217 gene of DGE results
To verify the previous DGE analysis, where LOC426217 was the most up-regulated in the deformed beaks, qRT-PCR was used to estimate the expression of this gene in 9 deformed and 9 normal beaks. As shown in Fig. 2, the relative expression of LOC426217 in deformed beaks was significantly higher than that in normal ones (P < 0.05). This was in good agreement with the DGE analysis.

Tissue expression profile of LOC426217
As shown in Fig. 3, LOC426217 gene was hardly expressed in brain, heart, bursa, or small intestine. The relative expression in beak was significantly higher than  (20) Note: χ 2 value means the test values of different genotypes to Hardy-Weinberg equilibrium. * (P < 0.05) means the loci were not in agreement with the Hardy-Weinberg equilibrium that in other tissues (P < 0.01). The results revealed that this gene maybe specifically expressed in beak tissue. In order to seek the structural changes of this gene, DNA sequencing was carried out subsequently.

DNA sequencing and SNPs detection
LOC426217 gene of 171 beak-deformed and 164 normal birds was amplified and sequenced. Eight SNPs were detected, including one locus in promoter region: -62G > T, and seven loci in coding region: 24 T > C, 36G > C, 192A > T, 204C > T, 222 T > C, 285G > T, and 363 T > C (Fig. 4). The loci in coding region were synonymous mutations resulting with no amino acid changing. These SNPs were presented in both beak-deformed and normal birds. The 18 birds used in the qRT-PCR analysis for DGE verification were also sequenced and their genotypes were shown in Table 3 and Additional file 1: Table S1.

Genetic diversity analysis
Genetic diversity was analyzed by PopGene (Version 1.31). All SNPs were firstly checked for HWE. In the normal birds (n = 164), T24C, G36C, A192T, T222C, G285T, and T363C were not in agreement with the HWE (P < 0.05), while in the birds with a deformed beak (n = 171), T24C, G285T, and T363C were not in agreement with the HWE. Two loci (G-62 T and C204T) of normal beaks and two loci (G-62 T and T363C) of deformed ones were low polymorphism (PIC < 0.25) (Tables 4 and 5).

Genotype frequency differences in the beak-deformed and normal birds
Chi-square tests were used to evaluate the genotype frequency differences of SNP loci in LOC426217 gene between two groups with a Bonferroni correction for the P-values. As shown in Table 6, the genotype frequency of G-62 T locus showed significant difference between two groups (P < 0.05), while the genotype frequencies of T24C, G36C, T222C and T363C loci showed highly significant differences (P < 0.01). There was no significant difference for the rest SNPs (P > 0.05).

Haplotype analysis
In the normal birds, two haplotype blocks of LOC426217 were identified: one block with SNPs T24C and G36C, and one block with SNPs A192T, C204T, T222C, and G285T. In the birds with a deformed beak, one block with SNP G36C and A192T was identified (Table 7 and Fig. 5). The haplotype blocks of deformed birds were different from those of the normal ones.

Discussion
The molecular genetic mechanism underlying beak deformity is likely to be very complex. Beak deformities of different forms (noticeably elongated, crossed, bent at right angles) have been documented in many wild birds. The molecular mechanism of beak deformity trait is not clear yet. For the wild birds, it is difficult to obtain the individuals for genetic study. Beak deformity was also found in various indigenous chickens and it is easy to collect individuals. This made the chicken a perfect model for the genetic study of this defect. Based on the previous DGE profiling and bioinformatics analyses, we identified a cluster of differentially expressed genes (DEGs) in the deformed and normal beaks. Some of the DEGs were quite extreme, especially the LOC426217 gene in the present study.

Validation and tissue expression profile
In order to validate the reliability and accuracy of DGE results, qRT-PCR was performed. The results showed that the expression of LOC426217 in the deformed beaks was significantly higher than that in the normal beaks, which was in great agreement with DGE analysis. The result of tissue expression profile also revealed that this gene was specifically expressed in beak tissue. LOC426217 is a member of the keratin family [24]. Keratin is a key gene family for maintaining normal cell morphology [16]. Variation of keratin structure can lead to beak deformity [25]. It is also an intermediate filament protein that has essential functions in maintaining the structural integrity of epidermis and its appendages [26], presumably including the beak. In addition, keratin is the main composition of the chicken beak. This may be the main reason for its high expression in beak tissue. According to our early observation of beak anatomy, the lower mandibles of the beaks were abnormal/asymmetry. It was reported that avian keratin disorder could result in gross over growth of the rhamphotheca [7]. The beak deformity caused by the excessive growth of one side of the lower mandibles maybe a result of abnormal high expression of LOC426217.

SNPs and haplotypes associated with beak deformity
Based on the case-control study design, the investigated birds used for sequencing were selected from the pure line BJY population according to the beak phenotype of the birds with no family structures.
Since the percentage of beak-deformity birds in the chicken population was 1 %-3 %, it was not easy to collect enough birds for gene sequencing. The beakdeformity birds were collected from two generations (about 5000 birds in each generation) and the normal birds were selected randomly from the same generation.
Using direct sequencing, eight SNPs located in the LOC426217 gene were detected, including one locus in promoter region and seven loci in coding region. These loci were found in both birds with deformed and normal beaks. Four loci (T24C, G36C, T222C, and T363C) were Note: χ 2 value means the tested values of different genotypes between two groups. * means significant difference between two groups (P < 0.05). ** means highly significant difference between two groups (P < 0.01). Corrected P-value using Benjamini & Hochberg method [21] was used for the Bonferroni correction  [27,28]. In this current study, one SNP (-G62T) was detected in the promoter region and the genotype frequency was significant different between two groups. This mutation could possibly affect the start, expression time and level of this gene. As shown in Table 3, the genotype of 9 normal birds were all TT, In the normal birds, two haplotype blocks were identified: one block with SNPs T24C and G36C, and one block with SNPs A192T, C204T, T222C, and G285T. In the birds with a deformed beak, one block with SNP G36C and A192T was identified while 3 out of the 9 beak-deformed were TG. Further validation study with larger sample size is needed to validate the mRNA expression profile of different genotypes. Similarly, the mutations in coding region are also very important. As we all know, exon is a nucleotide sequence in DNA that carries the code for the final mRNA molecule and thus defines the amino acid sequence during protein synthesis. If the base pairs changed in coding region, the structure and function of protein could possibly change, such as missense mutation or silent mutation [29][30][31][32]. In this study, all the loci in coding region were synonymous mutations with no amino acid and protein changing. Synonymous mutations, which do not alter the protein sequence, have also been shown to affect protein function [33,34] and play key roles in human diseases [35]. Therefore, the mutations detected here may also affect this gene function. Sometimes, the phenotype altering or diseases occurring was not associated with one single SNP changing but related to several SNPs. Global patterns of DNA sequence variation (haplotypes) defined by common SNPs have important implications for identifying disease and traits [36,37]. A previous study showed that the ALX1 haplotype has contributed to diversification of beak shapes among the Darwin's finches and, thereby, to an expanded utilization of food resources [13]. In this present study, different haplotypes were analyzed between two groups. It indicated that the beak deformity trait might be related to the specific haplotype block which was only existed in beak-deformed chickens. To sum up, SNPs and haplotypes described here were interesting and worthy of further study.

Conclusions
To the best of our knowledge, this is the first time that LOC426217 was studied as an important candidate gene for beak deformity in birds. The over-expression of LOC426217 may be the cause of beaks malformation. The genotype frequency of SNPs at G-62 T, T24C, G36C, T222C, and T363C loci showed significant differences between deformed and normal birds, and might be used as candidate SNP markers for this trait. The specific haplotype block in the deformed group could be served as a potential linkage marker for this trait. Further functional verification studies like over-expression or RNA interfere of LOC426217 during the embryonic development are required to reveal its roles in the beak malformation.