Knocking-out the human face genes TBX15 and PAX1 in mice alters facial and other physical morphology

DNA variants in or closed to the human TBX15 and PAX1 genes have been repeatedly associated with facial morphology in independent genome-wide association studies, while their functional roles in determining facial morphology remains to be understood. We generated Tbx15 knockout (Tbx15-/-) and Pax1 knockout (Pax1-/-) mice by applying the one-step CRISPR/Cas9 method. A total of 75 adult mice were used for subsequent phenotype analysis, including 38 Tbx15 mice (10 homozygous Tbx15-/-, 18 heterozygous Tbx15+/-, 10 wild-type WT) and 37 Pax1 mice (12 homozygous Pax1-/-, 15 heterozygous Pax1+/-, 10 WT mice). Facial and other physical morphological phenotypes were obtained from three-dimensional (3D) images acquired with the HandySCAN BLACK scanner. Compared to WT mice, the Tbx15-/- mutant mice had significantly shorter faces (P=1.08E-8, R2=0.61) and their ears were in a significantly lower position (P=3.54E-8, R2=0.62) manifesting an “ear dropping” characteristic. Besides these face alternations, Tbx15-/- mutant mice displayed significantly lower weight as well as shorter body and limb length. Pax1-/- mutant mice showed significantly longer noses (P=1.14E-5, R2=0.46) relative to WT mice, but otherwise displayed less obvious morphological alterations than Tbx15-/- mutant mice did. Because the Tbx15 and Pax1 effects on facial morphology we revealed here in mice are largely consistent with previously reported TBX15 and PAX1 face associations in humans, we suggest that the functional role these two genes play on determining the face of mice is similar to the functional impact their human homologues have on the face of humans. Author Summary Several independent genome-wide association studies (GWASs) on human facial morphology highlighted DNA variants in or closed to TBX15 and PAX1 with genome-wide significant association with human facial phenotypes. However, direct evidence on the functional involvement if these genes in the development and determination of facial morphology has not been established as of yet. In the current study, our in vivo gene editing experiments in mice for two well-replicated human face TBX15 and PAX1 genes provide novel evidence on the functional involvement of these two genes in facial and other physical morphology in mice, at least. Tbx15-/- mice showed a shortened facial length and manifesting an ear dropping characteristic, Pax1-/- mice showed an increased nose length. Our geometric morphometrics analysis further indicate that there are significant facial morphology differences between groups (Tbx15-/-and Tbx15+/-, Tbx15-/- and Tbx15+/+, Tbx15+/- and Tbx15+/+, Pax1-/- and Pax1+/+). We provide the first direct functional evidence that two well-known and replicated human face genes, Tbx15 and Pax1, impact facial and other body morphology in mice. The general agreement between our findings in knock-out mice with those from previous GWASs suggests that the functional evidence we established here in mice may also be relevant in humans.


Introduction
Human facial morphology represents a set of highly variable, multidimensional, highly correlated, symmetrical, and strongly heritable complex phenotypes. Unveiling the genetic basis of human facial variation is of fundamental and applied value in developmental biology, evolutionary biology, human genetics, medical genetics, and forensic genetics.
Several independent genome-wide association studies (GWASs) on human facial morphology highlighted DNA variants in or closed to TBX15 and PAX1 with genomewide significant association with human facial phenotypes [1][2][3][4][5][6]. In particular, Xiong et al. found DNA variants near TBX15 and PAX1 significantly associated with human face length and nose width, respectively [1]. Adhikari et al. reported SNPs in TBX15 to be significantly associated with two human ear traits, folding of antihelix and antitragus size [2]. Claes et al. and White et al. showed TBX15 and PAX1 intron variants to be significantly associated with human facial shape [3,6]. Adhikari et al. and Shaffer et al. reported significant association of PAX1 intergenic variants with human nose wing breadth [4,5].
Taking all available genetic knowledge together allows concluding that TBX15 and PAX1 play a role in facial morphology in humans. However, direct evidence on the functional involvement if these genes in the development and determination of facial morphology has not been established as of yet. Although previous in vivo experiments revealed that Tbx15 -/mutant mice exhibited lower weight, glucose intolerance, and obesity on high fat diets [16][17][18], these studies did not explore the effect of Tbx15 on facial morphology in mice. Moreover, no functional work has been reported for Pax1 in mice, or other functional evidence for PAX1 in humans.
In this study, to enhance our functional understanding on the involvement of Tbx15 and Pax1 in facial morphology in mice, we respectively knocked out these genes by applying the one-step CRISPR/Cas9 method [19]. On mutant and wild-type (WT) mice, we acquired 3D images through the HandySCAN BLACK scanner. Based on the collected digital imagery, we quantitatively assessed facial morphology as pairwise Euclidean distances between a set of anatomically meaningful facial landmarks, and also quantified other morphological phenotypes such as weight, body length, fore and hind limb length. Finally, the morphological differences between mutant mice and WT mice were compared using appropriate statistical tests.

Facial phenotypes
For all 75 F2 mice, digital 3D surface data were collected by use of the HandySCAN BLACK scanner. For quality control purposes, we compared the accuracy of the scanner with a millimeter ruler using a cylindrical object. The height and the diameter measured by the ruler (135.71 and 84.55 mm) were highly consistent with those measured by the scanner (135.70 and 84.58 mm), with the maximal difference smaller than 0.03 mm (S3 Fig). A total of 17 facial landmarks were selected for shape analysis.
Among the 17 selected landmarks, 9 overlapped with the 13 facial landmarks in humans as described in the previous face GWAS by Xiong et al (Fig 3A-3B, Fig 4C) [1]. The remaining eight landmarks were all on the mouse ears, and there were no corresponding landmarks in previous GWASs of human face and ear morphology. A total of 136 Euclidean distances connecting all possible pairs obtainable from the 17 facial landmarks were estimated and considered as facial phenotypes in subsequent analyses ( Fig 3A-3B). In all WT mice, the facial phenotypes were largely normally distributed (S4 Table), and no significant differences between males and females were observed (S5 Table). As expected, symmetric facial phenotypes showed higher correlations than non-symmetric ones (S6 Table), further supporting the reliability of the obtained phenotype dataset.
Compared with their WT counterparts, the Tbx15 -/mutant mice had significantly shorter faces, which affected 22 phenotypes as characterized by shorter distances between ear root and nose, mouth landmarks, and their ears had a significantly lower position, manifesting an "ear dropping" characteristic. (12 phenotypes, including L1-L8, Fig 3C and 4D, S7 Table). The most significant face shortening effect in the Tbx15 mutant mice was observed for L6-L15, for which the Tbx15 genotype accounted for 61% of the phenotype variance (0.64 mm per mutant gene). The most significant ear dropping effect in the Tbx15 mutant mice was observed for L2-L5 (0.95 mm per mutant gene), with the Tbx15 genotype explaining 62% of the phenotype variance ( Figure 4D, S7 Table). Principal components (PCs) were derived from 136 facial Euclidean distance.
An unsupervised clustering analysis of PC1 and PC2 could completely separate Tbx15 -/and WT mice ( Fig 5A) and a t-test of PC1 demonstrated significant differences between the heterozygote group and the WT group (P<.001) as well as between the heterozygote group and the Tbx15 -/group (P<.001, Fig 5B). These results overall demonstrated significant facial differences among Tbx15 -/-, Tbx15 +/and WT.
For both genes, the trends of the facial differences observed between mutant and WT mice were largely consistent with the genetic association effects observed in previous human GWASs (Fig 4) Table). The failure of observing an asymmetric effect of Tbx15 in our experiments might be explained by a pronounced effect of gene silencing compared with the potential regulatory effects of the SNPs. For PAX1, previous GWASs on facial shape in humans mainly found effects on nose width (Fig 4B-4C) [1,4,5]. In our experimental mice study, we found that Pax1 mutant mice had a significantly increased nose length (L13-L15, Fig 3D, Fig 4E, S8 Table). No significant changes in nose width were observed in Pax1 -/mice, which may be due to insufficient sample size and facial structure differences between mice and humans.

Differences in other physical phenotypes between mutant and wild-type mice
Next to the effects of Tbx15 and Pax1 gene editing on facial phenotypes in mice, we additionally investigated the effects on other physical morphological phenotypes such as body weight, body length, and length of fore and hind limbs. Compared with WT mice, Tbx15 -/mice had significantly lower body weight (p<0.01), shorter body length (p<0.001), and shorter limb length (P <0.01), (Fig 6). For Pax1, the body weight, body length, limb length of the mutant mice were also reduced compared to WT mice, but not statistically significantly so (S1 Figure).

Protein sequence conservation between mice and humans
The Tbx15 protein sequences are composed of 602 amino acid sequences in both human and mice. An amino acid sequence alignment analysis revealed that the amino acid sequences in both species had very high identity (98.67%, S2A Figure)

Expression of TBX15 and PAX1 in human embryonic cranial neural crest cells
Embryonic cranial neural crest cells (CNCCs) arise during weeks 3-6 of human gestation from the dorsal part of the neural tube ectoderm and migrate into the branchial arches. The later form the embryonic face, consequently establishing the central plan of facial morphology and determining species-specific and individual facial variation [22,23]. Because the is no expression data in CNCC cells of mice, we took the advantage of the available CNCC RNA-seq data [24] and other public RNA-seq datasets from GTEx [25] and ENCODE [26] in humans and investigated preferential expression patterns in CNCCs. Indeed, both TBX15 and PAX1 showed significant preferential expressions in CNCCs compared with that in other cell types (S9 Table).
These results indicated important functional roles of these two genes in facial morphogenesis during the early embryogenesis in humans, and likely also in mice.

Discussion
Our in vivo gene editing experiments in mice for two well-replicated human face TBX15 and PAX1 genes provide novel evidence on the functional involvement of these two genes in facial and other physical morphology in mice, at least. Tbx15 -/mice showed a shortened facial length and manifesting an ear dropping characteristic. In addition, they showed reduced weight, shortened body and limb length which are consistent with the findings from a previous mice study [17]. The TBX15 gene plays a major role in the development of the mesoderm of all vertebrates [27]. The complete inactivation of the mouse Tbx15 in mice and TBX15 mutations in human lead to severe bone deformities [8,11,28]. Notably, previous literature also showed that TBX15 affects the development of head, limbs, vertebrae and ribs by controlling the number of mesenchymal precursor cells and cartilage cells [8], while the Tbx15 effect on the face was not studied before.
Our new findings together with previous ones support our conclusion, that Tbx15 plays an important role in the development of facial and limb morphology in mice. Because the facial shape effects we observed in the Tbx15 -/mutant mice are largely in line with the TBX15 associations effects in humans previously reported in different GWASs, we suggest that the functional impact of Tbx15 on facial shape we revealed here in mice also exists in humans, which needs to be confirmed by future work.
Moreover, Pax1 -/mice showed an increased nose length in our study. Earlier studies revealed that Pax1 is a key transcription factor affecting cartilage development and regulates the expression of cartilage-related genes in the early stages of development [14]. Homozygotes exhibit variably severe morphological alterations of vertebral column, sternum, scapula, skull, and thymus, with reduced adult survival and fertility and some heterozygotes show milder skeletal abnormalities [29][30][31] , while the functional effect of Pax1 on facial morphology was not studied before. Since the nasal bone contains cartilage component, our findings of Pax1 determining nose morphology in mice may be explained by its effect on the development of nose cartilage, which needs to be further explored. In our study, we did observe reduced body weight, body length, and limb length in the Pax1 -/mice compared to WT mice, but the effect was not statistically significant, may be because of small effect and thus insufficient sample size.
Although the facial shape effects we observed in the Pax1 -/mutant mice are not fully consistent with PAX1 associations effects in humans previously reported in different GWASs, as in humans PAX1 associations were seen with nose width, while our in vivo data show that Pax1 affects nose length in mice, in both species this gene is involved in nose-related phenotypes. Therefore, we suggest that the functional effect of Pax1 on facial shape we revealed here in mice is similar to the effect of PAX1 in humans, which needs to be confirmed by future work.
In addition, these two genes may interact with development-related genes to affect the morphological development.Tbx15 was co-expressed with Shox2 genes in multiple tissues in humans and mice [32]. In human, the deletion of the SHOX region on human sex chromosomes caused short stature and Shox2 plays an important role in the face and body structure formation [33]. Homozygous mice (Shox2 +/+ ) had been embryonic lethal (63%). The craniofacial development of homozygous in embryonic was slow, and they showed less body weight and shorter body length [29][30][31]. Tbx15 and Shox2 gene may jointly affect the development of the facial and other morphology in humans and mice. There are protein interactions between Pax1 and Meox-1 in humans and mice [34].In human, Meox-1 protein was a mesodermal transcription factor that plays a key role in somitogenesis and was specifically required for sclerotome development. It is required for maintenance of the sclerotome polarity and formation of the cranio-cervical joints [35,36]. In mice, Meox-1 plays an important role in somatic cell development and limb muscle differentiation. Homozygotes exhibited hemi-vertebrae, rib, vertebral, and cranial-vertebral fusions [29][30][31].In all, Tbx15 and Pax1 may work together with other development-related genes, thus affecting the morphological development in humans and mice, which need to be explored in the future.
In this study, we investigated the function of these two particular human face genes in mice via genes knock-out experiments, because several independent GWASs repeatedly highlighted that intronic and intergenic DNA variants in these two genes demonstrate significant associations with facial shape phenotypes in humans [1][2][3][4][5][6].
However, functional data to explain these statistical associations was missing thus far and is provided for the first time via direct in-vivo experimental evidence in our present study. In this study, we selected the 9-week-old F2 generation mice because by 9-weeks they become sexually mature before engaging complex behaviors such as fertility, lactation, and reproduction, and their morphological characteristics are well developed. Furthermore, our littermate-born design prevented potential noises caused by unobserved environmental factors. In addition, the inclusion of the heterozygote mice in the analysis further expanded the sample size and was helpful in observing potential trend effects. In addition, we adopted reliable 3D measurement methods.
Kaustubh et al analyzed the effect of EDAR on ear morphology by acquiring twodimensional coordinates [2]. Because mice facial phenotypes are consists of facial spatial distances between all combinatorial pairs of the 17 facial landmarks, the threedimensional method can accurately quantify the appearance of mice and obtains mice 3D images with an accuracy of 0.03mm,which has a wide range of applications in animal morphology study.
This study focused on facial morphology phenotypes in mice, which wellcorresponded to the human facial phenotypes investigated in previous GWASs. Future studies may further investigate whether the observed genetic effects can be attributed to bone, fat or muscles.
In conclusion, we provide the first direct functional evidence that two well-known and replicated human face genes, Tbx15 and Pax1, impact facial and other body morphology in mice. TBX15 seems to affect the face globally, while PAX1 is mainly affecting the nose, with such a pattern being similar in both mice and humans. The general agreement between our findings in knock-out mice with those from previous GWASs suggests that the functional evidence we established here in mice may also be relevant in humans. In future research studies, functional effects of other human face genes highlighted in previous GWASs should be investigated to take our increasing knowledge on human face genetics from the statistical associational level to the next level of functional proof.

Samples
The mouse strain was C57BL/6J. As the in vitro transcription template, Cas9 and sgRNAs were amplified and purified with the gene knockout strategy. Fertilized eggs obtained from super-ovulated females (C57BL/6J, 4 weeks) mated with males (C57BL/6J, 7-8 weeks) were microinjected with mixtures of Cas9 mRNA and sgRNA (Fig 1, S1 Table). The injected fertilized eggs were cultured to day 2 and transferred to female mice. Finally, positive mice with a 20 bp frameshift mutation in the second exon of the TBX15 gene were obtained. Positive mice with the second exon of the PAX1 gene completely knocked out were also obtained. The discrimination criterion was the presence of a 20 bp frameshift mutation in the second exon of the TBX15 gene and the second exon of the PAX1 gene. The mice were raised in a pathogen-free environment and bred according to SPF animal breeding standards. The ambient temperature was 20-25 ℃, and humidity 40% -70%. The F0 generation-positive mice were mated with wild-type C57 mice to obtain the heterozygous F1 generation mice.
The heterozygous F1 mice were then mated with a female-to-male ratio of 2:1 to produce the homozygous F2 generation mice. The wild type, heterozygous, and homozygous mice were used for phenotyping when sexually mature at 9 weeks. The use of laboratory animals (SYXK 2019-0022) was licensed by the Beijing Municipal Science and Technology Commission.

Genotyping
Mice tail tissue of 0.3 cm 3 was cut into a 1.5 ml centrifuge tube, then added with 98 μL mouse tail lysate and 2 μL proteinase K before incubated in a metal bath at 55 ℃ for 30 min for full lysis. After that, the proteinase K was inactivated by placing the tube in a 95 °C metal bath for 10 min. The mixture was centrifuged at 12,000 r/min for 5 min, and the supernatant was collected as the DNA template for genotyping. The Tbx15 and Pax1 primers with sequences provided in Table S2 were

Protein sequence analysis
The Tbx15 and Pax1 protein sequences were obtained from the UniProt database [37,38]. Multiple alignments of amino acid sequences were performed using the DNAMAN software. The analysis of the protein conserved domains was performed using PfamScan software [39].

Gene expression analysis
RNA-seq data from the study of Prescott et al., GTEx and ENCODE were used to compare the differences in gene expression levels between these two genes in CNCCs and in other types of cells.

Statistical analyses
(1) 10-week-old heterozygous mice were mated with a female-to-male ratio of 2 : 1, and birth rate of the mice with different genotypes was analyzed using the Student's t-test analysis. The body weight, body length, and limb length of the mice were analyzed with multiple linear regression analysis.  (1) where y i represents a facial characteristic, α represents a fixed effect, ß ADD represents the effect of the genotype on facial phenotype, ß SEX represents the gender effect on facial phenotype, ε i represents the residual error, G i represents the genotype, SEX represents gender. The genotype and sex assignment methods are shown in Formulas (2) and (3), respectively.
(3) To adjust for the multiple testing of multiple phenotypes, the Bonferroni method was used to correct the effective number of independent variables, which was estimated using the Matrix Spectral Decomposition (matSpD) method [40] (4) Data analysis and statistical inference were performed using the R 3.