Analysis of DNA Methylation Profiles in Mandibular Condyle of Chicks With Crossed Beaks Using Whole-Genome Bisulfite Sequencing

Abstract

Crossed beaks have been observed in at least 12 chicken strains around the world, which severely impairs their growth and welfare. To explore the intrinsic factor causing crossed beaks, this study measured the length of bilateral mandibular ramus of affected birds, and investigated the genome-wide DNA methylation profiles of normal and affected sides of mandibular condyle. Results showed that the trait was caused by impaired development of unilateral mandibular ramus, which is extended through calcification of mandibular condyle. The methylation levels in the CG contexts were higher than that of CHG and CHH, with the highest methylation level of gene body region, followed by transcription termination sites and downstream. Subsequently, we identified 1,568 differentially methylated regions and 1,317 differentially methylated genes in CG contexts. Functional annotation analysis of Gene Ontology and Kyoto Encyclopedia of Genes and Genomes showed that these genes were involved in bone mineralization and bone morphogenesis. Furthermore, by combining the WGBS and previous RNA-Seq data, 11 overlapped genes were regulated by both long non-coding RNA and DNA methylation. Among them, FIGNL1 is an important gene in calcification of mandibular condyle. Generally, because the affected genes play key roles in maintaining mandibular calcification, these changes may be pivotal factors of crossed beaks.

Introduction

The beak consists of the maxillary and mandible, which are the main facial feature of birds (Gottschaldt and Lausmann, 1974). Crossed beak is a deformity defined earlier as misalignment of the upper and lower beak (Pomeroy, 1962), and the prevalence ranging from 0.2 to 7.4% was documented in at least 12 chicken strains around the world (Landauer, 1938, 1956; Joller et al., 2018; Hong et al., 2019; Shi et al., 2020). In addition, this trait exists in about 30% of 114 Chinese native chicken strains, according to our survey (Shi et al., 2020). Our previous study also showed that crossed beaks are frequently presented after hatch, and the crossed angle had been more and more serious with age until 56 days (Shi et al., 2020). Generally, chicks with crossed beaks have reduced feed intake (Benkman and Lindholm, 1991), inhibited growth (Chen et al., 2011; Giuseppe et al., 2015), poor performance (Hong et al., 2019), and shorter survival (Bai et al., 2018a, b; Joller et al., 2018), which is a great problem for the birds.

Crossed beak is a complex trait regulated by many genes, and its heritability was estimated to be 0.1 (Bai et al., 2018a). The genetic determinants of the complex trait have been studied at the genomic (Bai et al., 2018a, b; Joller et al., 2018), transcriptional (Bai et al., 2014), and translational levels (Sun et al., 2019). However, the genetic determinants of crossed beaks remain incompletely understood. DNA methylation is an epigenetic regulatory mechanism, which mediates numerous biological processes, such as growth, development and genomic imprinting (Smith and Meissner, 2013). Li C. W. et al. (2015) collected embryos and post-hatched chicks to study the level of global DNA methylation, and found that spatiotemporal specific epigenetic alterations could be critical for the late development of chick embryos and neonates. Xiao et al. (2016) evaluated DNA methylation in mandibular head cartilage in rat, and identified that 440 consistently changed genes in early, middle, and late phases of temporomandibular joint osteoarthritis, and 80% of which were hypomethylated and related to cell cycle regulation. In addition, the genome-wide methylation profile of bone revealed differentially methylated regions (DMRs) in osteoporosis and osteoarthritis, which enriched in genes associated with cell differentiation and skeletal embryogenesis (Delgado-Calle et al., 2013). These studies indicate that DNA methylation plays an important role in bone development. However, little is known about the expression patterns and potential roles of DNA methylation in beak development, especially in the complex genetic disease of crossed beak.

A previous study has shown that crossed beaks from 14 to 70 days of age were characterized with impaired development of unilateral mandibular ramus, and mandibular condyle is the growth center for the mandibular ramus extension (Shi et al., 2020). Moreover, Fukaya et al. (2017) discovered that unilateral IGF-1 injection extensively up-regulated the genes including RUNX2, COL2, and IHH in the mandibular condyle, and induced endochondral growth and a lateral shift of the mandible to the response side. This observation identifies the expression of molecular asymmetry may determines morphological left-right asymmetry in beaks. Based on above, we suspect that the epigenetic regulatory mechanisms of bilateral mandibular condyles of crossed beak chicks may be different, leading to asymmetric calcification of mandible. Therefore, this study systematically analyzed the DNA methylation profiles on both side mandibular condyle of crossed beak chicks using whole-genome bisulfite sequencing (WGBS) technology, aiming to provide new insights into the genetic basis of crossed beaks.

Results

Mandibular Length and Body Weight of Normal and Affected Birds

The affected birds of 7 days of age was caused by asymmetric length of bilateral mandibular ramus (Figure 1A). In particular, the left-side mandibular ramus with left mandibular curvature was shorter than the right-side (P < 0.01; Figure 1B). However, there was no difference among the normal right-side ramus of affected birds and two side ramus of normal birds. Meanwhile, the body weight of affected birds was lower than those of normal ones (P < 0.01; Figure 1C), which indicated that the beak deformity significantly decreased the growth.

FIGURE 1

Shi et al., 2020) with left mandibular curvature. MC means mandibular condyle; MR means mandibular ramus. (B) The length of bilateral mandibular ramus of crossed beaks with left mandibular curvature (n = 16) and normal beaks (n = 16). Results are expressed as means ± standard deviation, means were compared by Student-Newman-Keuls multiple-range tests. (C) Body weight of 7 day-old chicks with crossed beaks (n = 16) and normal beaks (n = 16). The body weight was analyzed using t-test, ** represents p < 0.01.

Genome-Wide DNA Methylation Profiling

The short left-side mandibular condyle of each four affected chicks were mixed in one composite sample and denoted as L group, their corresponding normal right-side condyle was mixed as R group. There are four replicates for each group in total. Global DNA methylation analysis of the four replicates was performed by WGBS with 35 × genome coverage and >99% conversion efficiency. A total of averagely 35.67, and 35.76 clean base were generated for the affected left-side condyle (L) and the normal right-side condyle (R), respectively. After filtering out low-quality data, approximately 71.43 million clean reads with Q30 ranging from 91.40 to 91.92% were generated for each replicate (Table 1). By aligning to chicken genome, the mapped reads were used for subsequent analysis with mapping rate ranging from 81.61 to 82.70%. Detailed quality of sequencing data is shown in Table 1.

All methylated genomic C sites were approximately 4.33% (Table 1). The methylation level of CG, CHH, and CHG (where H is A, C, or T) was significant different. In L group, genome-wide methylated cytosine (mC) levels were 89.66, 1.69, and 8.66% for CG, CHG, and CHH, respectively, and those of R group were 89.56, 1.71, and 8.74% (Supplementary Figure 1).

A violin graph was drawn with points representing different methylation levels. The CG methylation levels were high with wide sections in the violin graph, while CHG and CHH methylation levels were low with narrow sections (Supplementary Figure 2). Chromosome methylation maps for all composite samples were plotted (Supplementary Figure 3). The results showed that most chromosomal cytosine hypermethylation was found in the CG context.

We took the 3,000 bp upstream of a gene as promoter region, made an overlap annotation on CpG islands with methylation levels >0.7 and mC coverage > 5×, but not including the hypermethylated CpG island with C-degree confidence less than 0.1. In two groups, hypermethylated CpG islands were found in the distal intergenic regions (Supplementary Figure 4).

To further compare the genome-wide distribution and the methylation levels of various functional genomic elements, the methylation status of three different regions were analyzed, including upstream, gene body, and downstream regions (Figure 2 and Supplementary Table 1). In the two groups, there was no significant difference among different genetic elements of the three mC contexts. However, the methylation levels in the CG context were higher than those in the CHG and CHH contexts, where the CHH context was hypomethylated except for the transcription start site (TSS), while CHG context was almost completely unmethylated. The DNA methylation levels in the CG context were the highest in gene body region, then followed by transcription termination sites (TTS) and downstream regions, with sites near the TSS showing the lowest level. The methylation levels gradually decreased from the upstream to the TSS and increased from the TSS to the gene body region. In contrast, the DNA methylation levels in the CHH context decreased from the TSS to the gene body region.

FIGURE 2

Characterization of DMRs

In total, 1,568, 7, and 1,153 DMRs were identified in CG, CHG, and CHH contexts, respectively. As compared to the normal side, 1,330 (721 CG, 3 CHG, and 606 CHH) were hypermethylated and 1,398 (847 CG, 4 CHG, and 547 CHH) hypomethylated in the affected side. The DMRs were mostly located at distal intergenic regions, followed by introns, regulatory regions, and exons (Figure 3). In the CG context, only 145, 10, and 21 DMRs were in promoters, 5′UTRs, and 3′UTRs, respectively. In addition, as shown in the heat maps in Figure 4, the results showed a clear separation between the left-side and right-side mandibular condyle of crossed beak chicks. Formation of 53.9 and 47.3% hypomethylated DMRs in CG and CHH contexts, respectively. More detailed information is listed in Supplementary Table 2.

FIGURE 3

FIGURE 4

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis of Differentially Methylated Genes (DMGs)

To explore the change in the methylation status of genes, the GO and KEGG databases were used to annotate (Supplementary Table 3). Because most of the DMGs is CG and CHH methylation context, we focused on CG and CHH methylation for the DMG functional enrichment analysis. The GO enrichment analysis indicated that the 1,317 DMGs in CG methylation were significantly enriched in negative regulation of bone mineralization, bone morphogenesis, osteoblast differentiation, and so on (Figure 5A and Supplementary Table 4). In detail, there were three genes from negative regulation of bone mineralization, i.e., SOX9, AHSG, and BCOR; four genes from bone morphogenesis, i.e., PAX1, MSX1, PAPPA2, and CITED2, and eight genes from osteoblast differentiation, i.e., FIGNL1, FZF9, and so on. Most of the potential target genes were enriched in KEGG pathways of glycosaminoglycan biosynthesis-chondroitin sulfate/dermatan sulfate, glycosaminoglycan biosynthesis-heparan sulfate/heparin and Wnt signaling pathway (Figure 5B and Supplementary Table 5). The DMGs involved in the three pathways are CHST3, CHST13, NDST2, GLCE, DKK1, PRICKLE1, LEF1, AXIN2, MAP3K7, LOC101748851, and so on. In addition, the interaction network of 76 DMGs from three pathways that mentioned above and all GO terms was generated using STRING software. As shown in Figure 6, SOX9, RUNX2, MSX1, AXIN2, DKK1, and LEF1 were identified as hub genes in the interaction network related to the calcification of mandibular condyle.

FIGURE 5

FIGURE 6

Functional enrichment analysis showed that 660 DMGs in CHH methylation were significantly enriched in neural crest cell migration, bone morphogenesis, and so on (Figure 7). More detailed results of the GO and KEGG analyses in CG and CHH methylation are shown in Supplementary Tables 6, 7.

FIGURE 7

Expression Pattern of FIGNL1, MSX1, RUNX2, SOX9, CHST3, and CHST13

The expression of potential genes related with calcification were evaluated using T-test (Figure 8). The quantitative real-time polymerase chain reaction (qRT-PCR) analysis showed that the expression of FIGNL1 and MSX1 in affected-side mandibular condyle were lower than that of normal-side (P < 0.05). However, the expression of RUNX2 and SOX9 were not significantly different between the two sides mandibular condyle of crossed beaks, but the trend showed that those of 0.67 and 0.82 times, respectively lower in the affected-side than the normal-side. Moreover, there was no difference in expression of CHST3 and CHST13 (P > 0.05).

FIGURE 8

Integration Analysis of DMGs, Differentially Expressed Genes (DEGs), and Differentially Expressed Long Non-coding RNAs (DE lncRNAs)

A total of 14 genes were identified as both DEGs and DMGs, using our previous RNA-seq data, which similarly compared between the affected left-side and the normal right-side mandibular condyle of affected chicks with left mandibular curvature (Supplementary Table 8), and 11 candidate genes of them were regulated by both lncRNA and DNA methylation. Moreover, most genes were associated with more than one lncRNA, with methylation differences mainly distributed in the distal intergenic and intron regions. More detailed results on the above genes are listed in Table 2.

To investigate the effect of DNA methylation on gene expression levels, we compared the trend between gene expression and methylation levels using the fragments per kilobase million (FPKM) value for the RNA-seq data and the difference in methylation levels between L and R WGBS data samples (Table 2). The results showed that the DNA methylation level in the promoter regions of FIGNL1 was opposite of that observed for their expression levels. Furthermore, qRT-PCR results showed that the expression levels of the FIGNL1 was down-regulated in the affected left-side mandibular condyle of crossed beak chicks (P < 0.01; Figure 9).

FIGURE 9

Discussion

Within recent years, the sporadic occurrence of crossed beaks has been described by researchers in wildbirds (van Hemert et al., 2012; van Hemert and Handel, 2016) as well as in chickens (Handel et al., 2010; Joller et al., 2018; Hong et al., 2019). These affected chickens are usually normal at hatch, and do not become apparent until 1–2 months old (Landauer, 1938). Meanwhile, the type of mandible deviating laterally from the longitudinal axis of the head was more frequent (Bai, 2017; Joller et al., 2018; Hong et al., 2019). In this study, bilateral mandibular ramus length of crossed beaks chicks was asymmetrical at 7 day of age, where the affected shorter one is shorter than normal, and the other one is similar to that of the normal chicks. These results are similar to our previous study that crossed beaks had short of unilateral mandibular ramus from 14 to 70 day of age (Shi et al., 2020). In addition, many researches indicated that the condyle is essential for mandibular growth, in particular for the enlargement of the ramus (Meikle, 1973; Pirttiniemi et al., 2009). Thus, we suspect that there may be different epigenetic regulatory mechanisms for the growth of affected-side and normal-side mandibular ramus in a crossed beak, and caused the bone growth differently. Meanwhile, the mandibular condyle as the key point of bone development is an ideal sample for studying molecular mechanism of this complex trait.

DNA methylation is an epigenetic regulation form with important roles in gene expression and tissue development (Zhang et al., 2019). Although bone DNA methylation has been analyzed in human (Wu et al., 2019), baboon femora (Housman et al., 2017), and rat (Villagra et al., 2002; Nagaoka et al., 2019), to our knowledge, this is the first systematic comparison of genome-wide DNA methylation profiles of bilateral mandibular condyle in chickens. In the present study, bilateral mandibular condyles of genome-wide methylation patterns were similar in functional genomic regions. However, there were differences among three mC contexts which might be related to differences in the sequences of different genetic elements (Fang et al., 2017). Approximately 4.3% of cytosine sites were methylated, which is lower than the cortical bones in mice from 6.06 to 6.48% (Wei et al., 2020). The highest proportion of CG methylation in this study was similar to that found in other species (Wei et al., 2020) and tissues (Li C. W. et al., 2015; Li J. X. et al., 2015; Zhang et al., 2018; Tan et al., 2020). Among the gene functional regions, TSSs presented the lowest methylation levels, which was consistent with the results in chickens’ liver (Tan et al., 2020) and blood (Zhang et al., 2018), but inconsistent in cortical bones of mice (Wei et al., 2020).

We compared the trend between gene expression and methylation levels for the RNA-seq data and the difference in methylation levels to determine whether DMGs play a determinative role in calcification of mandibular condyle. The results showed that FIGNL1 was regulated by both lncRNA and DNA methylation. The DMRs of this gene is located in the promoter region, and the trend of DNA methylation levels is in contrast with its expression. As shown in a previous study, FIGNL1 is a subfamily member of the ATPases associated with diverse cellular activities protein family, which plays an important role in the inhibition of osteoblast proliferation and the stimulation of osteoblast differentiation (Park et al., 2007). The over-expression of FIGNL1 could reduce the proliferation of calvarial cells, and enhance the mRNA expression of osteocalcin and alkaline phosphatase. Therefore, FIGNL1 is highly expressed in the cells of mineralized tissue and plays a critical function in the formation of hard tissue. In contrast, interference with FIGNL1 significantly increased the proliferation of osteoblasts and decreased the expression of osteocalcin and alkaline phosphatase (Park et al., 2007). At present, there are few studies on epigenetic modifications or chromatin accessibility of FIGNL1. Although, the mechanism of crossed beaks has been studied at the genomic, transcriptional, and translational levels, the genetic determinants of crossed beaks remain incompletely understood. The joint analysis neither results in very consistent conclusions, mostly because the fact that the specific samples varies among studies. With the deepening of the research on the crossed beak, we discovered that the crossed beak was caused by the unilateral short of mandibular ramus, which extends through calcification of mandibular condyle (Shi et al., 2020). Based on the above, the mandibular condyle was studied and the FIGNL1 was identified by integration analysis of RNA-seq and WGBS. The results in this study indicated that, hypermethylation of the promoter regions may inhibit FIGNL1 expression on affected-side mandibular condyle of crossed beaks, and the expression of FIGNL1 was regulated by two lncRNAs. Therefore, FIGNL1 is an important gene in the process of mandibular calcification, and these results provided two candidate lncRNA and methylation marker for a new regulatory mechanism of calcification in mandibular condyle.

In this study, RUNX2, SOX9, and MSX1 are the hub genes identified through the interaction network associated with mandibular condylar calcification. RUNX2 is a transcription factor essential for skeletal development. Osteoblasts are completely absent in RUNX2 knockout mice, which indicates that RUNX2 is an essential transcription factor for osteoblast differentiation (Komori et al., 1997). SOX9 is a master transcription factor that participates in sequential events in chondrogenesis by regulating a series of downstream factors in a stage-specific manner. A previous study indicated that the physiological down-regulation of SOX9 in hypertrophic chondrocyte is associated with up-regulation of osteoblast-associated genes (Lefebvre, 2019). In transgenic mice expressing SOX9, the number of chondrocytes transdifferentiating into osteoblasts was markedly reduced (Lui et al., 2019). Moreover, SOX9 can also physically interact with RUNX2 and may thereby delay the master osteogenic actions of this RUNT-domain transcription factor (Zhou et al., 2006). In this study, the methylation of RUNX2 and SOX9 in distal intergenic and intron regions, respectively, may be involved in the formation of crossed beaks. Recent advances in techniques to study genome-wide methylation patterns have facilitated the identification of significant DNA methylation in intergenic and genebody regions. It is speculated that methylation within these non-promoter regions regulate alternative promoters, RNA processing, transposable elements such as long and short interspersed elements, and non-coding RNAs (Kulis et al., 2013). Intergenic DNA hypomethylation that results from dysfunctional trans-regulatory pathway (Weinberg et al., 2019). In addition, DNA methylation within intergenic regions is a mechanism regulating microRNAs (Pheiffer et al., 2016) and lncRNAs (Bermejo et al., 2019). The expression of RUNX2 and SOX9 were not significantly different between two sides mandibular condyle, however, the trend showed that those of 0.67 and 0.82 times, respectively lower in the affected-side than the normal-side mandibular condyle. Generally, DNA methylation is one of epigenetic mechanism which regulates gene expression. It remains to be further analyzed whether the above genes are regulated by non-coding RNA and histone modification.

MSX1 and PAX1 are two DMGs enriched in bone morphogenesis terms. MSX1 is a homeobox transcriptional factor and involved in limb-pattern formation and craniofacial development, specifically in tooth formation (Satokata and Maas, 1994; Chung et al., 2010; Nassif et al., 2014). Previous studies have reported that the most striking feature of MSX1 mutation is the inhibition of the mandibular basal convexity and absence of endochondral ossification in the mandibular condyle (Satokata et al., 2000; Orestes-Cardoso et al., 2009). Meanwhile, MSX1 is an upstream and downstream regulator for the bone morphogenetic protein BMP2 and BMP4 signaling pathway, respectively (Maxence et al., 2013), which stimulates trabecular bone metabolism and controls the collagen-based mineralization process (Nassif et al., 2014). In this study, methylation of the MSX1 in 3′UTRs regions may be related to the down-regulated of MSX1 in affected-side mandibular condyle of crossed beaks. PAX1 indirectly promotes the early stages of chondrogenic differentiation (Rodrigo, 2003), and PAX1-misexpressing chondrocytes exhibited abnormal cell morphology (Murtaugh et al., 2001; Takimoto et al., 2013). Moreover, Adhikari et al. (2016) and Shaffer et al. (2016) revealed that PAX1 play roles in craniofacial development or face syndromes. In this study, PAX1 had the DMR in exon region, which may be responsible for the calcification. Based on the above, these mentioned genes play important roles in calcification of mandibular condyle, and differentiation and regulation of them through DNA methylation might be one of the mechanisms that determine the difference of mandibular ramus length in crossed beak chicks. Nonetheless, the epigenetic mechanisms involved in the regulation of these genes and genetic regions involved in bone morphogenesis require further study.

This study systematically described the genome-wide DNA methylation patterns of mandibular condyle in chicks for the first time. FIGNL1, several important DMRs/DMGs, and pathways were emphasized to be related with calcification of mandibular condyle in crossed beaks. The results provide valuable data for further understanding the genetic and epigenetic mechanisms of this malformation.

Materials and Methods

Animals and Sample Collection

The study was conducted according to the local ethical guidelines and met the requirement of the institutional animal care and use committee (No. IAS2020-8). As the incidence of crossed beaks did not differ between male and female progeny based on our previous study (Bai, 2017), 32 female Beijing-You chickens including 16 normal beak chicks and 16 affected chicks with left mandibular curvature were used in this study. All birds were incubated contemporally and kept in the same environment without beak-trimming.

As the incidence increased quickly since 7 day of age (Shi et al., 2020), all chicks were weighed at the age of 7 day in this study, and the length of bilateral mandibular ramus were measured from the photo by Digimizer 5.3.4 MedCalc software (Ostend, Belgium). Meanwhile, two sides mandibular condyles of 16 affected chicks were dissected, temporarily frozen in liquid nitrogen and stored at -80°C.

The short left-side condyle of each four affected chicks were mixed in one composite sample and denoted as L group, their corresponding normal right -side condyle was mixed as R group. There are four replicates for each group in total. Then, genomic DNA was isolated from mandibular condyle tissues of each composite sample replicated using the phenol-chloroform method.

Library Preparation

The DNA concentration and quality were determined by NanoDrop (NanoDrop Technologies, Wilmington, DE, United States) and agarose gel electrophoresis before library construction. Four DNA libraries for L and R groups, respectively were constructed. Equal amounts of genomic DNA (2 μg per sample) were fragmented to 400–500 bp by ultrasonication, followed by adenylation and end-repair. The selected fragments were treated with bisulfite and then amplified by PCR to generate the sequencing libraries.

WGBS and Identification of DMRs

The library was sequenced using an IlluminaHiSeqTM2500 platform (Biomarker Technologies, Beijing, China), and the peak signal was transformed into sequence data by base calling, following which the raw reads were quality-filtered to obtain the clean reads. First, 3′ adapter sequence were trimmed. Then, reads with >10% unknown bases (N) and those of low quality (more than 50% of bases with a PHRED score ≤5) were removed. The Q30 and GC content were also calculated.

The clean reads were aligned to the chicken genome (GRCg6a) and the bisulfite mapping of methylation sites was performed using Bismark software. The duplicates were reads that aligned with the same region of the genome, and can estimated the sequencing depth and coverage. The bisulfite conversion rate is the percentage of methylated clean reads to the total number of clean reads in the genome. The binomial distribution test for each C site was used to confirm C site methylation by screening conditions for coverage ≥4 × and false discovery rate (FDR) < 0.05.

To identify DMRs between bilateral mandibular condyle of crossed beak chicks, we referred to a previously reported model (Lister et al., 2011) to estimate the methylation level. All C sites with reads coverage more than 10 × were used for DMR analysis performed in MOABS (Sun et al., 2014). DMRs were defined by the presence of at least three methylation sites in the region, and in which the difference in methylation levels was >0.1 for CHG and CHH context, and >0.2 for CG context, and the P-value from Fisher’s exact test was < 0.05. We annotate the DMR regions using ChIPseeker, and gene overlapped with at least one DMR was defined as DMG (Wang et al., 2020).

Function Enrichment Analysis

GO enrichment analysis of DMGs was implemented by the GOseq R packages based on the Wallenius non-central hypergeometric distribution (Young et al., 2010). KOBAS software was used to analyze the significance of DMGs enrichment in the KEGG pathway (Kanehisa et al., 2016). Pathways with a P-value < 0.05 were considered to be significantly enriched. The STRING database¹ was used to analyze interaction networks of DMGs (Franceschini et al., 2012).

Integration Analysis of DMGs, DEGs, and DE lncRNAs

Many DEGs and DE lncRNAs previously screened between left-side and right-side mandibular condyle of crossed beak chicks with left mandibular curvature using the Illumina platform. For trans target genes, we calculated the Pearson correlation coefficient (>0.9) and significant P-value (< 0.01) for the expression levels of each lncRNA-mRNA pair. For cis target genes, we identified chromosomal co-expressed genes within 100 kbps upstream and downstream of DE lncRNAs. Thereafter, the integrated analysis of the DMGs, DEGs (FDR < 0.05, and | log₂FoldChange| ≥ 1.5), and DE lncRNAs (P-value < 0.05, and | log₂FoldChange| ≥ 1.2) were further integrated analysis.

Validation by Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA from bilateral mandibular condyle of affected birds (n = 4) with left mandibular curvature were reverse transcribed into cDNA using the PrimeScript RT Reagent Kit (TaKaRa, Dalian, China) following the manufacturer’s instructions. qRT-PCR was performed on the ABI QuantStudio 7 Flex Real-Time Detection System (Life Technologies Holdings Pte Ltd., United States) using KAPA SYBR Fast universal qPCR kit (Kapa Biosystems, Boston, United States). Using GAPDH as a reference, relative-expression levels of six genes (FIGNL1, MSX1, RUNX2, SOX9, CHST3, and CHST13) were quantified using the 2^–ΔΔCt method (Livak and Schmittgen, 2001). The primer sequences are listed in Table 3.

Statistical Analysis

Data were analyzed by SAS 9.1 (SAS Institute Inc., Cary, NC). Means of mandibular length were compared by Student-Newman-Keuls multiple-range tests when a significant difference was detected. The body weight, and expression data were analyzed using T-test. The results are presented as the means ± standard deviation. A P-value < 0.05 (^∗) and P-value < 0.01 (^∗∗) implied a statistically significant difference and highly significant difference, respectively.

References

• 1

  AdhikariK.Fuentes-GuajardoM.Quinto-SánchezM.Mendoza-RevillaJ.Camilo Chacón-DuqueJ.Acuña-AlonzoV.et al (2016). A genome-wide association scan implicates DCHS2, RUNX2, GLI3, PAX1 and EDAR in human facial variation.Nat. Commun.7:11616. 10.1038/ncomms11616

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  BaiH. (2017). Genome-wide association studies and genome-wide detection of copy number variations for beak deformity in chickens.Ph. D. thesis.Beijing: Chinese academy of agricultural sciences.

  • Google Scholar
  • View reference in article

• 3

  BaiH.SunY. Y.LiuN.LiuY.XueF. G.LiY. L.et al (2018b). Genome-wide detection of CNVs associated with beak deformity in chickens using high-density 600K SNP arrays.Anim. Genet.49226–236. 10.1111/age.12652

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BaiH.SunY. Y.LiuN.XueF. G.LiY. L.XuS. S.et al (2018a). Single SNP- and pathway-based genome-wide association studies for beak deformity in chickens using high-density 600K SNP arrays.BMC Genomics19:501. 10.1186/s12864-018-4882-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BaiH.ZhuJ.SunY. Y.LiuR. R.LiuN.LiD. L.et al (2014). Identification of genes related to beak deformity of chickens using digital gene expression profiling.PLoS One9:e107050. 10.1371/journal.pone.0107050

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BenkmanC. W.LindholmA. K. (1991). The advantages and evolution of a morphological novelty.Nature349519–520. 10.1038/349519a0

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BermejoJ. L.HuangG.ManoochehriM.MesaK. G.SchickM.SilosR. G. (2019). Long intergenic noncoding RNA 299 methylation in peripheral blood is a biomarker for triple-negative breast cancer.Epigenomics1181–93. 10.2217/epi-2018-0121

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  ChenB. L.HaithK. L.MullensB. A. (2011). Beak condition drives abundance and grooming-mediated competitive asymmetry in a poultry ectoparasite community.Parasitology138748–757. 10.1017/S0031182011000229

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  ChungI.HanJ.IwataJ.ChaiY. (2010). Msx1 and Dlx5 function synergistically to regulate frontal bone development.Genesis48645–655. 10.1002/dvg.20671

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  Delgado-CalleJ.FernándezA. F.SainzJ.ZarrabeitiaM. T.SañudoC.García-RenedoR.et al (2013). Genome-wide profiling of bone reveals differentially methylated regions in osteoporosis and osteoarthritis.Arthritis. Rheum-US65197–205. 10.1002/art.37753

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  FangX.ZhaoZ.YuH.LiG.JiangP.YangY. (2017). Comparative genome-wide methylation analysis of longissimus dorsi muscles between Japanese black (Wagyu) and Chinese Red Steppes cattle.PLoS One12:e0182492. 10.1371/journal.pone.0182492

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  FranceschiniA.SzklarczykD.FrankildS.KuhnM.SimonovicM.RothA.et al (2012). STRING v9.1: protein-protein interaction networks, with increased coverage and integration.Nucleic Acids Res.41808–815. 10.1093/nar/gks1094

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  FukayaS.KanzakiH.MiyamotoY.YamaguchiY.NakamuraY. (2017). Possible alternative treatment for mandibular asymmetry by local unilateral IGF-1 injection into the mandibular condylar cavity: Experimental study in mice.Am. J. Orthod. Dentofac.152820–829. 10.1016/j.ajodo.2017.05.023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  GiuseppeV.MullensB. A.MenchJ. A. (2015). Relationships between beak condition, preening behavior and ectoparasite infestation levels in laying hens.Poult. Sci.941997–2007. 10.3382/ps/pev171

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  GottschaldtK.LausmannS. (1974). The peripheral morphological basis of tactile sensibility in the beak of geese.Cell Tissue Res.153477–496. 10.1007/BF00231542

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  HandelC. M.PajotL. M.MatsuokaS. M.van HemertC.TerenziJ.TalbotS. L.et al (2010). Epizootic of beak deformities among wild birds in alaska: an emerging disease in North America?Auk127882–898. 10.1525/auk.2010.10111

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  HongY.PangY.ZhaoH.ChenS.TanS.XiangH.et al (2019). The morphology of cross-beaks and BMP4 gene expression in Huiyang Bearded chickens.Animals9:1143. 10.3390/ani9121143

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  HousmanG.HavillL. M.QuillenE. E.ComuzzieA. G.StoneA. C. (2017). Assessment of DNA methylation patterns in the bone and cartilage of a nonhuman primate model of osteoarthritis.Cartilage10335–345. 10.1177/1947603518759173

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  JollerS.BertschingerF.KumpE.SpiriA.von RotzA.Schweizer-GorgasD.et al (2018). Crossed beaks in a local Swiss chicken breed.BMC Vet. Res.14:68. 10.1186/s12917-018-1398-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  KanehisaM.SatoY.MorishimaK. (2016). BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences.J. Mol. Biol.428726–731. 10.1016/j.jmb.2015.11.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  KomoriT.YagiH.NomuraS.YamaguchiA.SasakiK.DeguchiK.et al (1997). Targeted disruption of Cbfa1 results in a complete lack of bone formation owing to maturational arrest of osteoblasts.Cell85755–764. 10.1016/S0092-8674(00)80258-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  KulisM.QueirosA. C.BeekmanR.Martin-SuberoJ. I. (2013). Intragenic DNA methylation in transcriptional regulation, normal differentiation and cancer.Biochim. Biophys. Acta18291161–1174. 10.1016/j.bbagrm.2013.08.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LandauerW. (1938). Notes on cross-beak in fowl.J. Genet.3751–68. 10.1007/BF02982143

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  LandauerW. (1956). Hereditary and induced cross-beak of fowl.J. Exp. Zool.13225–38. 10.1002/jez.1401320103

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  LefebvreV. (2019). Roles and regulation of sox transcription factors in skeletogenesis.Curr. Top. Dev. Biol.133171–193. 10.1016/bs.ctdb.2019.01.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  LiC. W.GuoS. S.ZhangM.GaoJ.GuoY. M. (2015). DNA methylation and histone modification patterns during the late embryonic and early postnatal development of chickens.Poult. Sci.94706–721. 10.3382/ps/pev016

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  LiJ. X.LiR. J.WangY.HuX. X.ZhaoY. Q.LiL.et al (2015). Genome-wide DNA methylome variation in two genetically distinct chicken lines using MethylC-seq.BMC Genomics16:851. 10.1186/s12864-015-2098-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  ListerR.PelizzolaM.KidaY. S.HawkinsR. D.NeryJ. R.HonG.et al (2011). Hotspots of aberrant epigenomic reprogramming in human induced pluripotent stem cells.Nature47168–73. 10.1038/nature13843

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  LivakK. J.SchmittgenT. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method.Methods25402–408. 10.1006/meth.2001

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  LuiJ. C.YueS.LeeA.KikaniB.TemnyckyA.BarnesK. M.et al (2019). Persistent Sox9 expression in hypertrophic chondrocytes suppresses transdifferentiation into osteoblasts.Bone125169–177. 10.1016/j.bone.2019.05.027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  MaxenceV. R.KamalB.StefanoM.GiuliaG.PaoloP.SimonettaA.et al (2013). BMP-mediated functional cooperation between Dlx5; Dlx6 and Msx1; Msx2 during mammalian limb development.PLoS One8:e51700. 10.1371/journal.pone.0051700

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  MeikleM. C. (1973). The role of the condyle in the postnatal growth of the mandible.Am. J. Orthod.6450–62. 10.1016/0002-9416(73)90280-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  MurtaughL. C.ZengL.ChyungJ. H.LassarA. B. (2001). The chick transcriptional repressor Nkx3.2 acts downstream of Shh to promote BMP-Dependent axial chondrogenesis.Dev. Cell1411–422. 10.1016/S1534-5807(01)00039-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  NagaokaM.MaedaT.MoriwakiS.NomuraA.KatoY.NiidaS.et al (2019). Petunidin, a B-ring 5′-O-methylated derivative of delphinidin, stimulates osteoblastogenesis and reduces sRANKL-induced bone loss.Int. J. Mol. Sci.20:2795. 10.3390/ijms20112795

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  NassifA.SenussiI.MearyF.LoiodiceS.HottonD.RobertB.et al (2014). Msx1 role in craniofacial bone morphogenesis.Bone6696–104. 10.1016/j.bone.2014.06.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  Orestes-CardosoS.NefussiJ. R.LezotF.OboeufM.PereiraM.MesbahM.et al (2009). Msx1 is a regulator of bone formation during development and postnatal growth: in vivo investigations in a transgenic mouse model.Connect Tissue Res.43153–160. 10.1080/03008200290000547

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  ParkS. J.KimS. J.RheeY.ByunJ. H.KimS. H.KimM. H.et al (2007). Fidgetin-like 1 gene inhibited by basic fibroblast growth factor regulates the proliferation and differentiation of osteoblasts.J. Bone Miner Res.22889–896. 10.1359/jbmr.070311

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  PheifferC.ErasmusR. T.KengneA. P.MatshaT. E. (2016). Differential DNA methylation of microRNAs within promoters, intergenic and intragenic regions of type 2 diabetic, pre-diabetic and non-diabetic individuals.Clin. Biochem.573281–286. 10.1016/j.clinbiochem.2015.11.021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  PirttiniemiP.PeltomakiT.MullerL.LuderH. U. (2009). Abnormal mandibular growth and the condylar cartilage.Eur. J. Orthodont.311–11. 10.1093/ejo/cjn117

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  PomeroyD. E. (1962). Birds with abnormal bills.Br. Birds5549–72.

  • Google Scholar
  • View reference in article

• 41

  RodrigoI. (2003). Pax1 and Pax9 activate Bapx1 to induce chondrogenic differentiation in the sclerotome.Development130473–482. 10.1359/jbmr.081234

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  SatokataI.MaasR. (1994). Msx1 deficient mice exhibit cleft palate and abnormalities of craniofacial and tooth development.Nat. Genet.6348–356. 10.1038/ng0494-348

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  SatokataI.MaL.OhshimaH.BeiM.MaasR. (2000). Msx2 deficiency in mice causes pleiotropic defects in bone growth and ectodermal organ formation.Nat. Genet.24391–395. 10.1038/74231

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  ShafferJ. R.OrlovaE.LeeM. K.LeslieE. J.RaffenspergerZ. D.HeikeC. L.et al (2016). Genome-wide association study reveals multiple loci influencing normal human facial morphology.PLoS Genet.12:e1006149. 10.1371/journal.pgen.1006149

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  ShiL.LiY. L.BaiH.LiD. L.WangP. L.JiangL. L.et al (2020). Phenotype characterization of crossed beaks in Beijing-You chickens based on morphological observation.Poult. Sci.995197–5205. 10.1016/j.psj.2020.07.046

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  SmithZ. D.MeissnerA. (2013). DNA methylation: roles in mammalian development.Nat. Rev. Genet.14204–220. 10.1038/nrg3354

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  SunD.XiY.RodriguezB.ParkH.TongP.MeongM.et al (2014). MOABS: model based analysis of bi-sulfite sequencing data.Genome Biol.15:R38. 10.1186/gb-2014-15-2-r38

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  SunY. Y.LiuN.BaiH.LiY. L.XueF. G.YeJ. H.et al (2019). Differential proteomic analysis to identify proteins associated with beak deformity in chickens1.Poult. Sci.981833–1841. 10.3382/ps/pey519

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  TakimotoA.MohriH.KokubuC.HirakiY.ShukunamiC. (2013). Pax1 acts as a negative regulator of chondrocyte maturation.Exp. Cell Res.3193128–3139. 10.1016/j.yexcr.2013.09.015

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  TanX.LiuR.XingS.ZhangY.LiQ.ZhengM.et al (2020). Genome-wide detection of key genes and epigenetic markers for chicken fatty liver.Int. J. Mol. Sci.21:1800. 10.3390/ijms21051800

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  van HemertC.HandelC. M. (2016). Elements in whole blood of NorthWestern crows (corvus caurinus) in Alaska, USA: no evidence for an association with beak deformities.J. Wildlife Dis.52713–718. 10.7589/2015-10-287

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  van HemertC.HandelC. M.O’BrienD. M. (2012). Stable isotopes identify dietary changes associated with beak deformities in black-capped chickadees (poecile atricapillus).Auk129460–466. 10.1525/auk.2012.12037

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  VillagraA.GutiérrezJ.ParedesR.SierraJ.PuchiM.ImschenetzkyM.et al (2002). Reduced CpG methylation is associated with transcriptional activation of the bone-specific rat osteocalcin gene in osteo-blasts^∗The contents are solely the responsibility of the authors and do not necessarily represent the official views of the National Institutes of Health.J. Cell Biochem.85112–122. 10.1002/jcb.10113.abs

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  WangM.FengS.MaG.et al (2020). Whole-Genome Methylation Analysis Reveals Epigenetic Variation in Cloned and Donor Pigs.Front.Genet.11:23. 10.3389/fgene.2020.00023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  WeiY.FuJ.WuW.WuJ. (2020). Comparative profiles of DNA methylation and differential gene expression in osteocytic areas from aged and young mice.Cell. Biochem. Funct.38721–732. 10.1002/cbf.3539

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  WeinbergD. N.Papollon-CavanaghS.ChenH.YueY.ChenX.RajagopalanK. N.et al (2019). The histone mark H3K36me2 recruits DNMT3A and shapes the intergenic DNA methylation landscape.Nature573281–286. 10.1038/s41586-019-1534-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  WuJ.DuY.SongJ.DangX.LiuR. (2019). Genome-wide DNA methylation profiling of hip articular cartilage identifies dif-ferentially methylated loci associated with osteonecrosis of the femoral head.Bone127296–304. 10.1016/j.bone.2019.06.021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 58

  XiaoJ.MengJ.GanY.LiY.ZhouC.MaX. (2016). DNA methylation profiling in different phases of temporomandibular joint osteoarthritis in rats.Arch. Oral. Biol.68105–115. 10.1016/j.archoralbio.2016.04.006

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 59

  YoungM. D.WakefieldM. J.SmythG. K.OshlackA. (2010). Gene ontology analysis for RNA-seq: accounting for selection bias.Genome Biol.11:R14. 10.1186/gb-2010-11-2-r14

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  ZhangW. Z.ZhangS. X.XuY. Y.MaY. L.ZhangD. X.LiX. Y.et al (2019). The DNA methylation status of Wnt and TGFβ signals is a key factor on functional regulation of skeletal muscle satellite cell development.Front. Genet.10:220. 10.3389/fgene.2019.00220

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  ZhangZ.DuH.BaiL.YangC.JiangX. (2018). Whole genome bisulfite sequencing reveals unique adaptations to high-altitude environments in Tibetan chickens.PLoS One13:e0193597. 10.1371/journal.pone.0193597

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 62

  ZhouG.ZhengQ. Z.EnginF.MunivezE.ChenY. Q.SebaldE.et al (2006). Dominance of SOX9 function over RUNX2 during skeletogenesis.Proc. Natl. Acad. Sci. U S A.10319004–19009. 10.1073/pnas.0605170103

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article