Targeted sequencing in FGF/FGFR genes and association analysis of variants for mandibular prognathism.

To identify variants of the genes in fibroblast growth factors/fibroblast growth factor receptors (FGF/FGFR) signal pathway that predispose to mandibular prognathism (MP) in the general Chinese population systematically.Targeted sequencing of the FGF/FGFR genes was conducted in 176 MP individuals and 155 class I malocclusion controls. The associations of common and rare variants with MP as a categorical phenotype and also continuous malocclusion phenotypes generated by principal component (PC) analysis were analyzed.One common variant, rs372127537, located in the 3'-untranslated region of FGF7 gene, was significantly related to PC1 (P  =  4.22 × 10), which explained 23.23% of the overall phenotypic variation observed and corresponded to vertical discrepancies ranging from short anterior face height to long anterior face height, after Bonferroni correction. Also, 15 other variants were associated with PC1-4, although not significant after multiple corrections (P < .05). We also identified 3 variants: rs13317 in FGFR1, rs149242678 in FGF20, and rs79176051 FGF12 associated with MP (P < .05). With respect to rare variant analysis, variants within the FGF12 gene showed significant association with MP (P  =  .001).Association between FGF/FGFR signaling pathway and MP has been identified. We found a previously unreported SNP in FGF7 significantly related to increased facial height. Also, rare variants within the FGF12 were associated with MP. Our results provide new clues for genetic mechanisms of MP and shed light on strategies for evaluating rare variants that underlie complex traits. Future studies with larger sample sizes and more comprehensive genome coverage, and also in other population are required to replicate these findings.

Introduction

Mandibular prognathism (MP) is a complex maxillofacial disorder, which imposes significant aesthetic and functional burden on affected individuals worldwide. The prevalence of MP varies greatly relative to the population examined. Epidemiology studies suggested that the prevalence of MP is highest in individuals of Asian descent (approximately up to 15%) and lowest in individuals of Caucasian descent (about 1%).[ Genetic components play an important role in the pathogenesis of MP. However, the inheritance pattern of MP is still controversial. Both monogenic (dominant or recessive) and polygenic inheritance manners have been suggested by different studies.[ With the evidence accumulated, it is now accepted by most researchers that MP is a polygenic disorder with both environmental and genetic risk factors attributed to its etiology.[ Although genetic linkage analysis and association studies have identified many genes and loci associated with MP,[ the genes underlying the risk of MP in the general population remain elusive, leaving some impetus to search for new candidate genes.

As characterized by overgrowth of the lower jaw with or without undergrowth of the upper jaw,[ MP is deemed to be caused by abnormal craniofacial morphogenesis. The genetic mechanisms of craniofacial development has been elucidated, with FGF, bone morphogenetic protein, sonic hedgehog, and many other signal pathways playing critical roles.[ In the past decade, the role of FGF signaling in craniofacial development has been extensively investigated. It has been shown that FGF signaling exerts an inductive impact on facial primordia formation. Fgfr1 and Fgfr2 spread widely in facial primordia, whereas the FGF ligands are present in restricted region: Fgf8, Fgf9, and Fgf10 are highly expressed at nasal pits, whereas expression of Fgf3, Fgf15, and Fgf17 is confined to the medial side of the nasal pits.[ Exogenous Fgf2 and Fgf4 could give rise to increased length of the cartilage rod formed in the frontonasal and mandibular mesenchyme.[ Inactivation of Fgf8 in the ectoderm of the first branchial in zebrafish results in almost complete loss of first-arch derived structures, including the mandible.[ In addition, FGF signaling also has a vital function in craniofacial skeletogeny. It is expressed in both endochondral and intramembranous bones and involved in modulating their development and growth.[

Given the crucial roles of FGF signaling in craniofacial development, it is not surprising that mutations in the FGF pathways are involved in various congenital bone diseases. It has been reported that gain-of-function mutations in Fgfr1 and Fgfr2 lead to craniosynostosis syndromes, such as Apert, Crouzon, and Pfeiffer syndrome, all of which often manifest MP phenotype.[ Recently, we have identified a novel mutation in FGF23, c.35C>A, strongly associated with MP.[ Moreover, 2 single-nucleotide polymorphisms (SNPs), rs2162540 and rs11200014, of FGFR2 are suggested to increase the risk for classes II and III skeletal malocclusion.[ According to these clues, it is reasonable to hypothesize that the variants of genes in FGF/FGFR signal pathway play a significant role in MP pathogenesis. The purpose of this study was to identify variants of the genes in FGF/FGFR signal pathway that predispose to MP in the general Chinese population systematically.

Materials and methods

Study subjects

In all, 176 subjects with MP (mean age 23.36 ± 9.83 years; 79 males) and 155 subjects with normal skeletal class I (mean age 23.01 ± 6.30 years; 60 males), who were seeking orthodontic treatment at the affiliated Stomatology Hospital of Tongji University, were recruited from January 2013 to September 2014 (Table 1). All participants were unrelated and were of Han Chinese ancestry. Written informed consent (including the release for dental records) was acquired from each participant or the parental guardians in case of the minors. This study was approved by the Human Ethics Committee of Tongji University and was conducted according to Declaration of Helsinki principles.

Table 1

Demographical characteristics of the cases and controls.

All individuals were first diagnosed by lateral cephalograms, in conjunction with orthodontic study models or visual inspection. Digital lateral cephalograms were taken by a dental X-ray equipment (Veraviewepocs X550, Kyoto, Japan), using a standardized technique with the patients’ jaws in centric occlusion. Individuals were diagnosed as MP if they had an ANB angle (Point A-Nasion-Point B) of the centric jaw relationship less than 0.0 degrees[ and a negative Wits appraisal greater than −2.0 mm.[ And the criteria of normal skeletal class I was defined as follows: an ANB angle range from 0.3 to 4.8 degree along with a Wits appraisal between −1.3 and 2.4 mm. Individuals who had previous orthodontic treatment, severe facial trauma, congenital abnormalities (eg, cleft lip and palate), or general physical disease (eg, endocrine diseases) were excluded from this study.

Approximately 5 mL of ethylenediaminetetraacetic acid-anticoagulated peripheral blood were obtained from each individual, and genomic DNA was extracted using a QIAamp DNA Blood Kit (QIAGENE GmbH, Hilden, Germany). All the samples were stored at <−80°C until analysis.

Cephalometric analysis

Pretreatment cephalometric tracing was performed using NemoCeph NX software (version 6.0, Nemotec, Madrid, Spain) by 2 independent orthodontists. All the lateral cephalograms were traced twice by each rater at least 2 weeks apart. Sixty-one cephalometric parameters digitized with 27 skeletal landmarks and 9 soft landmarks were used for phenotyping. Inter-rater and intrarater reliability was tested by intraclass correlation method as described previously.[ A Procrustes routine was implemented to wipe out variants irrelevant to shape.[ After that, the Procrustes residuals was employed for principal component (PC) analysis to determine the most significant aspects in the data.

Targeted region sequencing and data analysis

The coding and flanking regions of 26 genes in the FGF/FGFR signaling pathway, approximately 91.3 kb, were selected and sequenced in this study. The targeted regions (shown in Supplemental Table S1) were captured according to the standard procedures using a customized NimbleGen capture array (Roche-NimbleGen Inc. Madison, WI) and then sequenced on the Illumina Hiseq2000 platform (Illumina Inc, San Diego, CA). The raw reads were aligned to the human reference genome (hg19) using Burrows-Wheeler Alignment tool v0.7.1 (http://maq.sourceforge.net), producing binary sequence alignment/map files containing various mapping information. The duplicate reads were then removed using Picard v1.137 (https://github.com/broadinstitute/picard/releases) and realigned using Genome Analysis Toolkit v3.4–46. The coverage, mean quality, and global depth of aligned reads were calculated by perl scripts based on the pileup files generated by SAMtools v1.2. Only the bases with ≥20 base quality were taken into account. Variants were called using SNPTools and annotated using the ANNOVAR software package (http://www.openbioinformatics.org/annovar/). Indels (insertion/deletion) were verified manually. To access reproducibility, 6 random selected samples (3 cases and 3 controls) were analyzed in duplicate and Sanger sequencing of the positive SNPs in this study was also conducted.

Statistical analysis

For common variants (minor allele frequency [MAF] ≥1%), the allelic and genotypic distributions of the case and control groups were compared using the Pearson chi-square test. Fisher exact test was used when the expected count was <5. Logistic regression analysis was used to identify the effects of the variants on MP with odds ratio (OR) and 95% confidence interval (CI). Linear regression analysis adjusting for age and sex was performed to identify the associations between the variants and PCs, explaining more than 5% of the total variance in the cephalometric data. For rare variants (MAF <1%), the cumulative variants within each gene region in cases and controls were compared by Fisher exact test. All statistical analyses were performed using SAS, version 9.2 (SAS institute Inc., Cary, NC), and 2-tailed P values <0.05 were considered statistically significant.

Results

Targeted sequencing data

The average sequencing coverage was 67× (interquartile range 43–87×). The concordant of the calling variants in duplicate samples were more than 99%. The variants with calling rates less than 95% or inconsistent with Hardy–Weinberg equilibrium (P < .01) in the control group were removed. Based on the MAF of the SNPs tested in the control group, the retained variants were then classified into 2 groups: common variants (MAF ≥0.01) and rare variants (MAF <0.01). Overall, 460 variants were identified across all sequenced individuals in the targeted regions, including 139 common variants and 321 rare variants. Among the 139 common variants identified, 43 of the variants conducted 15 haplotype blocks with the fact that the variants in the same block were in almost complete linkage disequilibrium (LD, r2 > 0.90, D’  =  1). In this case, only 1 variant from each block was chosen for further analysis randomly and the final number of common variants to be analyzed was 111.

Association analysis of common variants

Among the 111 common variants, we only presented those significantly associated with MP and PCs which explaining more than 5% of the facial variation. The results of Sanger sequencing of the positive SNPs in this study were in accordance with the raw results. No common variants within the 26 genes were significantly associated with MP after Bonferroni correction for multiple testing (cut-off P value  =  .05/111  =  4.50 × 10–4). Only 3 SNPs reached nominal significance (P < .05), including rs13317 in the 3’-untranslated region (UTR) of FGFR1, rs149242678 in the 5’-UTR of FGF20, and rs79176051 in the intron of FGF12. As shown in Table 2, the genotypic and allelic frequencies at rs13317, rs149242678, and rs79176051 were significantly different between case and control groups. The C allele of rs149242678 decreased the risk of MP, whereas the C allele of rs13317 and the G allele of rs79176051 increased the risk of MP, with ORs of 4.379 (1.234–15.536), 0.666 (0.470–0.944), and 0.369 (0.168–0.812), respectively.

Table 2

The associations between SNPs identified in FGF/FGFR signaling pathway and MP.

The results of the principal component analysis revealed 5 PCs accounted for 73.4% of the total variance, and each of them represented 23.2%, 19.8%, 13.6%, 10.0%, and 6.7% of the total variance, respectively (Fig. 1). As shown in Table 3, 1 common variant (rs372127537; P  =  4.22 × 10–4) in the 3’-UTR of FGF7 gene was significantly associated with PC1, which depicted vertical discrepancies ranging from short anterior face height to long anterior face height, after Bonferroni correction. This variant was also associated with PC2 (P  =  4.96 × 10–3), which captured mandibular shapes ranging from a larger mandibular body, a large ramus height to a small mandibular body, a small ramus height, and PC4 (P  =  1.35 × 10–2), which captured inclination of mandibular incisor, although the associations were not significant after multiple corrections. We also detected 3 other SNPs associated with PC1, 2 other SNPs associated with PC2, and 2 other SNPs associated with PC4 with nominal significance. In addition, 6 SNPs associated with PC3, which captured the inclination of maxillary incisor, and 3 SNPs associated with PC5, which captured horizontal discrepancies of the jaws with respect to anterior skull base, were also identified at the P < .05 significance level.

Figure 1

Examples of individuals at opposite extremes of the distributions for each of the 5 principal components. PC1 explained 23.2% of the observed facial variation and corresponded to variations within the anterior facial height. PC2 explained 19.8% of the observed facial variation and corresponded to differences in the ramus height (Co-Gn in mm) and mandibular body length (Ar-Go in mm). PC3 corresponded to characteristics of incisor inclination and explained 13.6% of the observed facial variation. PC4 corresponded to the mandibular incisor angulation and explained 10% of the observed facial variation. PC5 corresponded to the position of the mandibular and accounted for 6.7% of the variation.

Table 3

The associations of common SNPs with continuous malocclusion phenotypes.

Association analysis of rare variants

Table 4 displays the association results of rare variants with MP for the 26 targeted genes. Compared with the controls, the MP group had more rare variants in FGF12 (P  =  .001). However, when restricted to nonsynomymous variants, none of the targeted genes showed association with MP.

Table 4

Burden test of rare variants.

Discussion

It is widely held that genetic components play an important role in MP. So far, numerous chromosomal loci implicated in MP pathogenesis has been reported, and also a host of genes that predispose to MP, such as EPB41, MATN1, COL2A1, MYO1H, TGFB3, LTBP2, ADAMTS1, DUSP6, FGFR2, and FGF23.[ Most of these studies were based on family linkage study. However, the polygenic nature of MP makes it possible to study its genetic mechanism by case-control design.

With the evidence accumulated, it is now accepted by most researchers that MP is a polygenic disorder with both environmental and genetic risk factors attributed to its etiology.[ Although genetic linkage analysis and association studies have identified many genes and loci associated with MP, the genes underlying the risk of MP in the general population remain elusive, leaving some impetus to search for new candidate genes. In the current study, we aimed to identify the association between variants in the FGF/FGFR signaling pathway and MP in MP cases and controls using target sequencing strategy and have found some novel variants in these genes associated with MP.

By analyzing common variants in the coding and flanking regions of 26 selected FGF/FGFR genes in 176 MP cases and 155 controls, we found that 1 SNP, rs372127537, was significantly associated with PC1 after Bonferroni correction. This SNP was also correlated with PC2 and PC4 with nominal significant difference. This SNP was located in the 3’-UTR of FGF7 gene and may influence the gene expression.[ FGF7 is a member of FGF family, which is known as a mediator of epithelial-mesenchymal tissue interactions in several organs.[ It may act directly to induce the formation of an additional apical ectodermal ridge in the ectoderm of the dorsal midline,[ and the apical ectodermal ridge is indispensable for limb outgrowth proceeds.[ Furthermore, FGF7 is also expressed in perichondrium of growth plate during bone formation.[ Perlecan, a prominent component of human cartilage, is a receptor for FGF7. Their interaction initiates cell signaling and subsequent down-line effects on cell proliferation and differentiation, thus coordinates chondrogenesis and angiogenesis during skeletal development.[ And condylar cartilage are key to the regulation of mandibular growth.[ Our results indicated that rs372127537 was significantly associated with an increased anterior and posterior height of mandible. This SNP may play an important role in the increase of facial height by affecting growth of condyle cartilage, hence mandible. The potential function of this SNP and the FGF7 gene in MP development needs to be studied further.

We also detected 3 SNPs in FGFR1, FGF12, and FGF20, respectively, associated with MP, with nominal significant difference. FGFR1 has an extensive function during craniofacial morphogenesis, and is almost involved in all the structures, including craniomaxillofacial skeleton, muscle, palate, tooth, and submandibular salivary gland.[ It modulates osteoblast differentiation as a positive regulator for skeleton formation.[ Mutation in FGFR1 leads to craniosynostosis syndrome (Pfeiffer syndrome), in which relative MP due to maxillary hypoplasia is a common finding.[ Loss-of-function variants in FGFR1 are also determined in patients with combined pituitary hormone deficiencies, which can be associated with complex phenotypes such as cranial/facial midline defects.[ Moreover, ColII-FGFR31ach transgenic mice exhibit shortened long bones and a domed-shaped skull, probably owing to craniofacial hypoplasia.[ And according to HaploReg v.4.1, rs13317 is on protein-binding region of CCAAT enhancer-binding protein β (CEBPB), which is demonstrated as an important determinant of osteoblast function and bone mass.[ p20C/EBP β (a dominant negative inhibitor of Cebpb) transgenic mouse exhibits significant bone volume reduction of mandible.[ Therefore, the SNP rs13317 in FGFR1 is presumed to interfere normal craniofacial shaping and result in MP. On the contrary, FGF20 is shown to be a major downstream effector of Eda, and affects Eda-regulated characteristics of tooth morphogenesis, including the number, size, and shape of teeth.[ And fgf20a is demonstrated to directly affect suture and skull development in zebrafish, and fgf20a deficiency also causes craniofacial defects similar to Albertson syndrome.[ The identified SNP rs149242678 was located in the 5’-UTR of FGF20 gene, and related to CCCTC-binding factor (CTCF) protein binding. CTCF is a heritable component of an epigenetic system regulating the interplay between DNA methylation, higher-order chromatin structure, and lineage-specific gene expression.[ Future work is needed to gain an insight into the role of these variants in MP pathogenesis.

With respect to FGF12, which has high sequence identity with the FGF family, but does not activate FGFRs, except for detecting common variant associated with MP, we also identified more rare variants in FGF12 in MP individuals compared with controls. Emerging researches have figured out that the common variants with given complex diseases are unable to explain fully of their genetic etiology. Thus, the rare variants are suggested to give rise to explaining a slice of these genetic diseases. Multiple rare variants have been examined affect complex traits strongly, especially the extremes of a disease.[ A recent study reveals significant association of cleft lip with variants in FGF12. Cleft lip can cause other dentofacial malformations in humans such as MP. Also, dysfunction of FGF12 is found to contribute to skeletal growth and development failure of grade II and III Kashin-Beck disease.[ There are no previous association studies with FGF12 and MP, and it is also the first time that we identified rare variants associated with MP in FGF12. However, as none of these effects survived multiple testing correction, further independent replication of these findings will be required.

In a previous research, we have identified a novel heterozygous mutation in FGF23, c.35C>A, associated with MP, in a Chinese pedigree by family linkage analysis.[ However, no significant association was found between variants of FGF23 and MP in this case-control association study. These suggests us that the genetic mechanisms of MP are complex,[ the main effects of the identified variants are considered to confer relatively small increments in risk, and explain only a small proportion of the heritability.[ The mutation in FGF23, c.35C>A, may only explain the development of MP in a small number of family samples. The contribution of FGF23C35A to MP may be small, and not the susceptibility gene of the group of MP individuals studied in this study.

Identification of the variants to MP in FGF/FGFR signal pathway is the first step to reveal the genes contributing to this disease. What necessary to be done is to detect the genetic architecture across more candidate genes, and to test the variants in large-scale individuals and verify the function of the genes pertinent to MP. Knowledge of these important genes will be helpful in elucidating the mechanisms of MP, and in improving diagnosis or even treatment by simple intervention strategies.

In this study, we have identified 1 common variant in FGF7 significantly associated with PC1, which demonstrated vertical discrepancies ranging from short anterior face height to long anterior face height. We also identified 3 other variants associated with MP and 15 other variants associated with PC1–4, although not significant after multiple corrections. Moreover, the rare variant within FGF12 showed significant association with MP. Future studies with larger sample sizes and more comprehensive genome coverage, and also in other population are required to replicate these findings, and further functional studies are also warranted.

Acknowledgments

We are indebted to participants for providing an excellence study sample, and we appreciate the dentists who took clinical examination.