Whole-exome sequencing identified four loci influencing craniofacial morphology in northern Han Chinese

Facial shape differences are one of the most significant phenotypes in humans. It is affected largely by skull shape. However, research into the genetic basis of the craniofacial morphology has rarely been reported. The present study aimed to identify genetic variants influencing craniofacial morphology in northern Han Chinese through whole-exome sequencing (WES). Phenotypic data of the volunteers’ faces and skulls were obtained through three-dimensional CT scan of the skull. A total of 48 phenotypes (35 facial and 13 cranial phenotypes) were used for the bioinformatics analysis. Four genetic loci were identified affecting the craniofacial shapes. The four candidate genes are RGPD3, IGSF3, SLC28A3, and USP40. Four single-nucleotide polymorphism (SNP) site mutations in RGPD3, IGSF3, and USP40 were significantly associated with the skull shape (p < 1×10−6), and three SNP site mutations in RGPD3, IGSF3, and SLC28A3 were significantly associated with the facial shape (p < 1×10−6). The rs62152530 site mutation in the RGPD3 gene may be closely associated with the nasal length, ear length, and alar width. The rs647711 site mutation in the IGSF3 gene may be closely associated with the nasal length, mandibular width, and width between the mental foramina. The rs10868138 site mutation in the SLC28A3 gene may be associated with the nasal length, alar width, width between tragus, and width between the mental foramina. The rs1048603 and rs838543 site mutations in the USP40 gene may be closely associated with the pyriform aperture width. Our findings provide useful genetic information for the determination of face morphology.


Introduction
Human facial shapes are species-specific and have individual specificity. The individual specificity of human facial shapes plays an important role in human personal identity and social activities (Koblenzer 2003;Zhang and Gao 2009). Facial features are affected by genetic environmental factors and their interactions. It has been shown that the heritability of most facial features exceeds 60% (Crouch et al. 2018;Djordjevic et al. 2016). The craniofacial shape has been found to be mainly determined by genetic factors (Relethford 1994;Relethford and Harpending 1994), and the shape of the skull determines the basic shape of the face (Lieberman et al. 2000). The correlation between the skull and facial shapes has also been used in forensic practice by forensic anthropologists for processes such as facial image imposition and craniofacial approximation (Claes et al. 2010;Yoshino et al. 1995). Research into the genetic basis of the skull shape is, therefore, crucial for further studies of the genetic variations of facial shapes. However, there is a lack of research on the molecular genetic basis of facial features.
To date, the published papers have shown whole-genome association analyses of human facial shapes (Adhikari et al. 2016;Cha et al. 2018;Claes et al. 2018;Cole et al. 2016;Lee et al. 2017;Liu et al. 2012;Paternoster et al. 2012;Qiao et al. 2018;Shaffer et al. 2016). Researchers have successfully applied these results to human facial modeling predictions and achieved satisfactory results (Claes et al. 2014).
Wei Wu and Guiying Zhai have contributed equally to this work and should be considered co-first authors.

3
The subjects of the published studies were Europeans, Latin Americans, Africans, and Koreans. Whole-genome association analyses of facial morphology have shown association with PAX3, PRDM16,TP63,C5orf50,COL17A1,DCHS2,RUNX2,GLI3,PAX1,EDAR,MAFB,PAX9,MIPOL1,ALX3,HDAC8,CACNA2D3,SCHIP1,PDE8A,WDR27,SOX9,DHX35,etc. However, the results were not consistent among the above-studied population (except PAX3, which was associated with facial morphology in a few of published papers), indicating that ethnic factors had a large impact on facial features. Therefore, the present study aimed to identify genetic variants influencing craniofacial morphology in northern Han Chinese through whole-exome sequencing (WES) and Sanger Sequencing.

Materials and methods
Study subjects 50 volunteers, including 13 males and 37 females, were recruited from three universities in the Jinzhou area, Liaoning, China. To minimize the impact of nongenetic confounding factors, such as region, ethnicity, age, and BMI, we strictly used the following criteria to select volunteers for this study: (1) no craniofacial surgery, trauma history and craniofacial deformity; (2) no genetic kinship each other; (3) normal Han adults who had been living in Liaoning for at least three generations and were between 20 and 30 years of age; and (4) normal body size with a BMI of 18.5-24. This study was approved by the Ethics Committee of Jinzhou Medical University. All participants provided written informed consent. Venous blood collection and a head CT scan were completed between March 2017 and July 2017 following the standard procedures.
Facial and skull phenotypic measurements using threedimensional (3D) CT scanning of the head were used to collect craniofacial phenotype information and to reconstruct the 3D face and skull images based on the difference in gray values of the acquired CT data. The craniofacial landmarks on the CT were then marked by group discussion as previously described (Evison et al. 2010). The Euclidean distance between the landmarks was calculated to replace the directly measured distance, which avoided the inaccuracies of mapping, rotation, and other issues caused by selection of a coordinate system. The convenience and reliability of this morphometric method have been proven by numerous studies (Lele andRichtsmeier 1991, 1995;Subhash Lele 1993;Starbuck et al. 2015).
The VG Studio MAX 2.2 software (Volume Graphics Co., Ltd., Heidelberg, Germany) was used to perform 3D facial and skull reconstruction using CT images from the 50 subjects. This software was also used simultaneously to measure the 3D coordinates of the facial and skull landmarks.

WES and whole-genome association analysis
WES and whole-genome association analysis were performed by Beijing Novogene Technology Co., Ltd. (Beijing, China, http://www.novog ene.com). The quality of DNA samples was evaluated by agarose gel electrophoresis analysis. The DNA samples that were accurately quantified as greater than 0.6 µg by Qubit were used for library construction. The Agilent's liquid chip capture system was applied to efficiently enrich the DNA containing all human exon regions. High-throughput deep sequencing was subsequently performed on the Illumina platform. The Agilent SureSelect Human All Exon V6 Kit was used for library construction and the capture experiments. After library construction, Qubit 2.0 was used for the preliminary quantification. The Agilent 2100 system was subsequently used to verify the library insert size. The Illumina platform was used for sequencing according to the effective concentration of the library and the data output requirements. High-throughput paired-end sequencing was performed (paired-end 150 bp, PE150).
After WES was completed, a bioinformatics analysis was performed that included three main aspects: sequencing data quality assessment, SNP detection, and whole-genome association analysis. The sequencing data quality control requirements were as follows: the sequencing error rate of each base position was less than 1%, the average Q20 ratio was greater than 90%, the average Q30 ratio was greater than 80%, the average error rate was less than 0.1%, the sequencing reads achieved an alignment rate of 95% or greater, and the read depth of the base at one position reached 10X or more.
The SNP testing was performed as follows: high-quality sequences were aligned with the human reference genome (GRCh37/hg19) to detect variant information in the sample, and the detected variations were analyzed and annotated.

Sanger sequencing
Sanger sequencing validation in this study was performed by the BGI Beijing Branch (http://www.genom ics.cn). Premier 3.0 was used for primer design; the primer sequences are listed in Table 5. The PCR products were purified using Millipore 96 purification plates. The sequencing reaction mixture was as follows: 2 μL of the buffer mix (Bigdye 3.1, 5× sequencing buffer in H 2 O, 2 μL of the purified PCR product, and 1 μL of the primer (5 mmol/L). The sequencing reaction program was as follows: 95 °C 15 s → (95 °C 15 s → 50 °C 5 s → 60 °C 90 s) × 35 cycles → stop reaction. Electrophoresis was used to detect the PCR-amplified products.

Statistical analysis
The GEMMA (Genome-wide Efficient Mixed Model Association algorithm) software was used to construct a linear mixed model and calculate the association between the 48 phenotypes and genotypes (SNP variant loci) (Zhou and Stephens 2012). A mixed linear model was constructed using GEMMA software in R studio: y = Wα + xβ + u + ε, y: n quantitative traits of individuals; W: covariate (gender); x: n individual genotypes; β: slope coefficient, the size of the genotype. A mixed linear model can eliminate confounding results more accurately, and can directly reflect the relationship between craniofacial phenotype and genotype. Combining other studies of GWAS correlation analysis with our stringent requirements for quality control, our p value threshold is set to 1E−6 to avoid the effects of other mixing factors on experimental results. Therefore, the p values of less than 1 × 10 −6 were regarded as statistically different.

Craniofacial phenotypic measurement
Using the facial photographs of 50 individuals, we accessed 35 facial features and 13 cranial features. The craniofacial phenotypic measurements of adult northern Han Chinese are shown in Table 2.

WES and whole-genome association analysis
A total of 23,011 SNP sites were detected by WES. The sequencing quality control is shown in Table 3 and Fig. 2. The association analysis between the 23,011 SNP loci and 48 phenotypes showed that a C > T missense mutation occurred at rs62152530 (NM_001144013) in RGPD3 (2q12). This SNP site mutation is closely associated with 4 facial phenotypes (sl-sn, obs_l-obi_l, sa_r-sba_r, sa_l-sba_l, and al_r-al_l). This mutation causes an Asp > Asn change in the protein amino acid sequence at position 111. A G > C missense mutation occurred at rs647711 (NM_001007237) in IGSF3 (1p13). This SNP site mutation is closely associated with one facial phenotype (sl-sn) and two skull phenotypes (GO_R-GO_L and ML_R-ML_L), and causes an Asp > Glu change in the protein amino acid sequence at position 1020. A T > C missense mutation occurred at rs10868138 (NM_001199638) in SLC28A3 (9q21). This SNP site mutation is closely associated with 4 facial phenotypes (n-sn, sl-sn, t_r-t_l, and al_r-al_l) and one Euclidean distance between the right alare and the left alare sal_r-al_r Euclidean distance from the right subalare to the right alare sal_r-al_l Euclidean distance from the right subalare to the left alare sal_l-al_r Euclidean distance from the left subalare to the right alare sal_l-al_l Euclidean distance from the left subalare to the left alare sal_l-sal_r Euclidean distance from the left subalare to the right subalare go_r-go_l Euclidean distance between the right gonion and the left gonion en_l-ex_l Euclidean distance from the left entocanthion to the left ectocanthion en_r-en_l Euclidean distance between the right entocanthion and the left entocanthion zg_r-zg_l Euclidean distance between the right zygion and the left zygion Morphological facial index The quotient of the height and the width of the morphological surface is multiplied by 100 Zygomatic mandibular index The width between the gonions is divided by the width of the surface and multiplied by 100 Nasal index The nasal width is divided by the nasal height and multiplied by 100 Skull phenotypes used in this study (13) N-RHI Euclidean distance from nasion to rhinion RHI-SP Euclidean distance from rhinion to spinale RHI-SS Euclidean distance from rhinion to subspinale FMO_R-FMO_L Euclidean distance between the right frontomalare and the left frontomalare OR_R-OR_L Euclidean distance between the right orbitale and the left orbitale MF_R-FMO_R Euclidean distance from the right maxillofrontale to the right frontomalare MF_L-FMO_L Euclidean distance from the left maxillofrontale to the left frontomalare MF_R-MF_L Euclidean distance from the right maxillofrontale to the left maxillofrontale NMJ_R-NMJ_L Euclidean distance between the right naso-maxillary junction and the left naso-maxillary junction LP_R-LP_L Euclidean distance between the right lateral piriform and the left lateral piriform ZG_R-ZG_L Euclidean distance between the right and the left zygion GO_R-GO_L Euclidean distance between the right gonion and the left gonion ML_R-ML_L Euclidean distance between the right mentale and the left mentale skull phenotype (ML_R-ML_L). This mutation causes a Tyr > Cys change in the protein amino acid sequence at position 113. G > A missense mutations occurred at the rs1048603 (NM_018218) and rs838543 (NM_018218) sites in USP40 (2q37). These two SNP site mutations are closely associated with one skull phenotype (NMJ_R-NMJ_L). The mutations caused a Cys > Thr change in the protein amino acid sequence at position 1123 and a Thr > Cys change at position 2033 (

Discussion
We identified the SNP loci of four genes which were significantly associated with the craniofacial shape of northern Han Chinese. This result is the first report on the genetics of northern Han Chinese craniofacial shapes, which has significance for understanding craniofacial growth and development and improving the accuracy of human facial modeling based on DNA. Our findings could not replicate the previously published results (Table 5).
Five SNP sites at 2q12, 1p13, 9q21, and 2q37 were found to be associated with craniofacial phenotypes in northern Han Chinese. These SNP sites are located in the four candidate genes RGPD3 (RANBP2-Like and GRIP Domain Containing 3), IGSF3 (Immunoglobulin Superfamily Member 3), SLC28A3 (Solute Carrier Family 28 Member 3), and USP40 (Ubiquitin Specific Peptidase 40).
RGPD3 (Ciccarelli et al. 2005) is located in the Ran-binding protein-related gene cluster on chromosome 2 (cytogenetic location: 2q12.2; genomic coordinates (GRCh37/ hg19) 107,021,136-107,084,832 and (GRCh38/hg38) 106,391,290-106,468,376). This gene is a protein-coding gene, and one of its related pathways is RNA transport. A study of new complex syndromes caused by the absence of 2q12.1-2q14.2 suggests that its absence may lead to a  deformed appearance in child patients, which is manifested as a flat arcus supraorbitalis, wide alveolar ridge, high forehead, low ear, and other facial features (Stoeva et al. 2011). Because the RGPD3 gene is located within the range of 2q12.1-2q14.2, RGPD3 is highly likely to be associated with these special facial features. In this study, we found that a C > T missense mutation occurred at rs62152530. This SNP site mutation is closely associated with the facial nasal length (sl-sn), ear length (obs_l-obi_l, sa_r-sba_r, and sa_l-sba_l), and nasal width (al_r-al_l), which partially overlap with the special facial features caused by the loss of 2q12.1-2q14.2, suggesting that RGPD3 is associated with the morphology of the nose and ears. IGSF3 is located on chromosome 1 (cytogenetic location: 1p13.1; genomic coordinates (GRCh38): 1: 116,574,397-116,667,754) and is a member of the immunoglobulin superfamily. However, studies have not detected the IGSF3 mRNA in peripheral blood leukocytes, the spleen, or the liver, suggesting that IGSF3 is not involved in immune functions (Saupe et al. 1998). Foster et al. revealed that loss of IGSF3 was associated with congenital nasolacrimal duct obstruction. RT-PCR analysis confirmed the presence of IGSF3 transcripts in the embryonic mouse lacrimal gland, and IGSF3 transcripts were detected in the lacrimal canaliculi and lacrimal glands 1 day after mice were born. Congenital lacrimal duct obstruction includes Hasner valve obstruction and congenital lacrimal duct atresia (including bony nasolacrimal duct atresia caused by bone malformation) (Foster et al. 2014). This finding suggests that IGSF3 may be involved in the regulation of craniofacial developmental. In this study, a G > C missense mutation was found at rs647711 in IGSF3, and this SNP mutation was closely associated with the facial nasal length (sl-sn), the width between the gonions, and the width between the mentales (GO_R-GO_L and ML_R-ML_L) of the skull. This result suggests that IGSF3 may not only participate in regulating development of the lacrimal ducts, but may also participate in the growth and development of other structures of the face and skull.
SLC28A3, also known as CNT3 (Ritzel et al. 2001), is a solute carrier family member located on chromosome 9 [cytogenetic location: 9q21.32-q21.33; genomic coordinates (GRCh38): 9: 84,273,368,633], and is involved in the regulation of various cellular processes, including neurotransmission, vascular tone, nucleoside transport and metabolism, and other processes. SLC28A3 has extensive specificity for the transport of pyrimidine and purine nucleosides (Barber et al. 2010). Studies have shown that SLC28A3 is associated with obesity, and it may be associated with the lipid-lowering response of statins (Comuzzie et al. 2012;Guallar et al. 2007). Body fat content has a great influence on facial morphological characteristics, and thus, SLC28A3 may have a certain correlation with facial morphology. In this study, a T > C missense mutation was found at rs10868138 in SLC28A3. This SNP site mutation is closely associated with the facial nasal length, the width between the tragus, the alar width (n-sn, sl-sn, t_r-t_l, and al_r-al_l), and the width between the mental foramina (ML_R-ML_L).
USP40 is located on chromosome 2 [cytogenetic location: 2q37.1; genomic coordinates (GRCh38): 2:233, 1 3 475,518-233,566,788] and belongs to a large family of cysteine proteases. USP40 modifies cellular proteins through ubiquitin by playing the role of a deubiquitinating enzyme (Quesada et al. 2004). USP40 is expressed in many organ tissues, and its biological functions, protein localization, and structural characterization are completely unknown. A few reports have linked this gene to Parkinson's disease or glomerular function (Li et al. 2006;Takagi et al. 2017;Zhao et al. 2012). In this study, we found G > A missense mutations at rs1048603 and rs838543 in the USP40 gene; these two SNP site mutations were closely associated with the width of the pyriform aperture. This result suggests that USP40 has a wide range of biological functions and may participate in the regulation of craniofacial development. Candidate genes for human facial variation and GWAS studies showed a few consistent results. Most of the novel genes identified in each study were different. In fact, racial factors are a major cause because the results come from different populations, such as Europeans, Americans, Africans and East Asians, etc. (Adhikari et al. 2016;Cha et al. 2018;Claes et al. 2018;Cole et al. 2016;Lee et al. 2017;Liu et al. 2012;Paternoster et al. 2012;Qiao et al. 2018;Shaffer et al. 2016). Even in the same race, such as in Europeans, few of the genes or loci identified in each study have been Fig. 3 The whole-exome region showing genome-wide significant associations with facial traits. The SNP site mutations at 1p13 (rs647711), 9q21 (rs10868138), and 2q12 (rs62152530) have stronger associations with the facial morphology of the Han people of northern China. A G > C missense mutation occurred at rs647711 (IGSF3). This SNP site mutation is closely associated with one facial phenotype (nose length) (a). A T > C missense mutation occurred at rs10868138 (SLC28A3). This SNP site mutation is closely associated with four facial phenotypes (nose length, tragus width, nasal width n-sn, sl-sn, t_r-t_l, and al_r-al_l) (b, d, f, j). A C > T missense mutation occurred at rs62152530 (RGPD3). This SNP mutation is closely associated with four facial phenotypes (nose length, ear length, and nasal width) (c, h, i, e). Here, the LocusZoom plot used the linkage disequilibrium (LD) information in the hg19 1000 Genome Nov 2014 ASN (East Asia) database to obtain the r 2 . Since no LD information was available in the region of the database where the rs62152530 site was located, the image of this locus did not have r 2 information Fig. 4 The whole-exome region showing genome-wide significant associations with the skull traits. A G > C missense mutation occurred at rs647711 (IGSF3). This SNP site mutation is closely associated with two skull phenotypes (width between the gonions and width between the mentales) (a, b). A T > C missense mutation occurred at rs10868138 (SLC28A3). This SNP site mutation is closely associated with one skull phenotype (width between the men-tales) (c). Missense mutations of G > A occurred at both rs1048603 (USP40) and rs838543 (USP40). These two SNP site mutations are closely associated with one skull phenotype (pyriform aperture width) (d, e). Here, the LocusZoom plot was generated using the linkage disequilibrium (LD) information in the hg19 1000 Genome Nov 2014 ASN (East Asia) database

Fig. 5
Sanger sequencing validated the mutations at the five SNP sites. A C > T missense mutation occurred at rs62152530 (RGPD3), a G > C missense mutation occurred at rs647711 (IGSF3), a T > C missense mutation occurred at rs10868138 (SLC28A3), and G > A mis-sense mutations occurred at both rs1048603 (USP40) and rs838543 (USP40). The Sanger sequencing results confirmed that the five SNP sites identified by WES, indeed, had mutations replicated (Adhikari et al. 2016;Lee et al. 2017;Paternoster et al. 2012). In addition, during the review process of this paper, Qiao et al. (2018) reported that, in the GWAS study of the discovery panel of Uyghurs, six significant loci were identified, four of which, rs1868752, rs118078182, s60159418 at or near UBASH3B, COL23A1, PCDH7, and rs17868256 were replicated in independent cohorts of Uyghurs or Southern Han Chinese, which are not insistent with our findings. In addition, the environment may be the likely cause of this situation, for example, even if the study participants in each study are Europeans, they are from different countries or from European descent/Americans. The inconsistency between Qiao's and our findings may also be due to ethnic and/or environmental factors in the Uighur or South Han and the Northern Han Chinese. Detection of effective genes or loci affecting craniofacial morphology cannot ignore sample size. The authors acknowledge the limitation of a small sample size; only 50 subjects were selected due to funding limitation. However, the findings through WES and whole-genome association analysis in the Han Chinese are novel and significant in the understanding of the genetic basis of facial development.
Another possibility is that the measurement methods and phenotypic descriptions in each study are different. For phenotypic and genomic association analysis, the accuracy of phenotypic measurement is crucial. In the previous studies, 2D method was often used to measure craniofacial morphology. In the 2D method, even with 3D reconstructions, the linear distance between the two projection points is measured on a plane parallel to the coronal or sagittal plane when measuring the distance between two landmarks. However, the phenotypic accuracy of 2D image measurement is limited. Liu et al. (2012) indicated that 2D image noise, image size, and pixel resolution had an impact on phenotypic measurements, and some marker point coordinates cannot be measured on a 2D image. Cha et al. (2018) also showed that the main limitation of their research was that facial measurements were derived from 2D images rather than 3D images (they subsequently measured more facial features such as angle, ratio, curvature, width, height, etc.). In the present study, we used geometric morphometrics to extract the 3D coordinates of two landmarks and calculated the Euclidean distance between the two points. We believe that this new method improves the accuracy of our study. In addition, to minimize the measurement error, in this study, the two authors performed three geometric measurements on all coordinate points, respectively, using the average of the measured data as the final ones.
Overall, the present study not only identified four loci influencing craniofacial morphology in northern Han Chinese, but also established an approach for analyzing the genetic basis of human facial traits through WES combined with 3D CT scanning. Compared with GWAS, WES display is more accurate and cost-effective, and the sequencing depth can reach 120x and more, which can meet the research of large target and multi-target area. In this study, WES can more directly reflect the effect of protein-coding regions on craniofacial phenotypes.
Taken together, our results show that RGPD3, IGSF3, SLC28A3, and USP40 have significant associations with the facial morphology of the northern Han people and that the associated phenotypes are concentrated in the nose and mandible regions. Although not the same genes as reported, it is consistent with the relevant regions. Of these genes, both IGSF3 and SLC28A3 are significantly associated with facial and cranial morphology, suggesting that genes that affect cranial morphology likely affect facial morphology, as well. Cranial morphology, in turn, should be incorporated into genomics research of facial morphology. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.