Genome-wide association study identiﬁes a new susceptibility locus for cleft lip with or without a cleft palate

Nonsyndromic cleft lip with or without a cleft palate (NSCL/P) is among the most common human congenital birth defects and imposes a substantial physical and ﬁnancial burden on affected individuals. Here, we conduct a case–control-based GWAS followed by two rounds of replication; we include six independent cohorts from China to elucidate the genetic architecture of NSCL/P in Chinese populations. Using this combined analysis, we identify a new locus at 16p13.3 associated with NSCL/P: rs8049367 between CREBBP and ADCY9 (odds ratio ¼ 0.74, P ¼ 8.98 (cid:2) 10 (cid:3) 12 ). We conﬁrm that the reported loci at 1q32.2, 10q25.3, 17p13.1 and 20q12 are also involved in NSCL/P development in Chinese populations. Our results provide additional evidence that the rs2235371-related haplotype at 1q32.2

N onsyndromic cleft lip with or without cleft palate (NSCL/ P) is among the most common human congenital birth defects 1 . Despite improvements in surgical correction techniques, this serious malformation has lifelong implications for health and imposes a substantial financial burden on affected families 2 . Previous genetic and epidemiological studies clarified that the aetiology of NSCL/P is multifactorial, with both genetic and environmental factors contributing to the phenotype [3][4][5] .
To date, four multi-stage, well-designed genome-wide association studies (GWAS) of NSCL/P have been conducted. In addition to the susceptibility gene IRF6, which was previously identified in candidate studies before the GWAS era 6,7 , Birnbaum et al. 8 identified a key locus for NSCL/P at 8q24.21 using 686 cases and 1,337 control subjects of European ancestry. An extension of this study, which involved 401 cases and 1,323 controls in addition to independent replication in European trios, identified two additional relevant loci at 10q25.3 and 17q22 ref. 9. Beaty et al. 10 applied a case-parent design to explore genomewide transmission patterns in 825 trios of European ancestry and 1,038 trios of Asian ancestry and validated the identified singlenucleotide polymorphisms (SNPs) with 331 Asian and 51 South Asian families, 1,149 European families and 434 South/Central American families, which resulted in the identification of two new susceptibility loci at 1p22.1 and 20q12. Subsequent validation of Beaty's study identified three other loci at 1p36.13, 3q12.3 and 17p13 that are significantly associated with NSCL/P 11 . By combining the discovery phase of previous GWAS 9,10 , Ludwig et al. 12 reported five additional loci (1p36.13, 2p21, 3p11.1, 8q21.3 and 13q31.1) that are associated with NSCL/P. These observations not only targeted multiple regions influencing the risk of NSCL/P but also suggested that European and Asian populations may have different causal variants in the same locus as a result of different genetic backgrounds 2 . Although these studies [8][9][10]12 were well designed, they mainly focused on susceptibility loci in people of European ancestry. Some of them involved individuals from Asia 10,12 , but they recruited samples from different countries and only some of them included individuals from China.
Thus, to further elucidate the genetic architecture of NSCL/P in Chinese individuals and to identify new susceptibility loci for NSCL/P, we conduct a GWAS followed by two rounds of replication in a Chinese population. Using this combined analysis, we identify a new locus at 16p13.3 associated with NSCL/P: rs8049367 between CREBBP and ADCY9 (odds ratio ¼ 0.74, P ¼ 8.98 Â 10 À 12 ). We confirm that the reported loci at 1q32.2, 10q25.3, 17p13.1 and 20q12 are also involved in NSCL/P development in Chinese populations.

Results
Genome-wide association analysis. In the initial discovery phase, we performed a GWAS scan in two cohorts (Huaxi and Nanjing, Supplementary Fig. 1; Supplementary Table 1). After quality control ( Supplementary Fig. 2), 842,556 autosomal SNPs were analysed (see Methods) in 858 cases and 1,248 controls (504 cases and 455 controls from the Huaxi cohort (GWAS I), and 354 cases and 793 controls from the Nanjing cohort (GWAS II)). We performed logistic regression in each cohort with adjustments for sex and the significant principal components of the corresponding cohort to test the additive effect of minor alleles of each SNP. Meta-analysis was used to combine the results of the two cohorts in the GWA scan. Because the controls in the Nanjing cohort were derived from imputed data from a previous study and did not include newborn infants, we applied a random-effect model in the meta-analysis when heterogeneity existed (P for Cochran's Q statistic o0.05) to avoid false-positive results caused by putative population stratification or bias derived from age-based differences and imputation error. Quantile-quantile plots of genome-wide association test statistics revealed minimal inflation; thus, substantial cryptic population substructure and differing genotypic variants between the cases and controls were unlikely after meta-analysis of the two cohorts (l Huaxi ¼ 1.03, l Nanjing ¼ 1.10 and l meta ¼ 1.06, Supplementary Fig. 3).
Thirty-four SNPs exhibited consistently significant associations on analysis during the discovery phase (see Methods). Among these SNPs, 18 representative SNPs were selected for subsequent evaluation, and the other 16 were ignored because of their high linkage disequilibrium (LD) (r 2 40.8) with at least one of the 18 selected SNPs (see Methods, Supplementary Table 2). Coupled with an additional 11 SNPs identified in previous studies (see Methods, Table 1), we included a total of 29 SNPs in the first validation conducted in Guangzhou (validation Ia), Huaxi (validation Ib) and Nanjing (validation Ic) (see Methods, Fig. 1 and Supplementary Fig. 4).
Replication and combined analysis. In the first round of replication, sixteen SNPs exhibited evidence of association with NSCL/P (Meta P value of validation I o0.05, Supplementary Table 3), and we then evaluated these SNPs in the Hui population (validation II), a racial minority in China. In the combined zFor SNPs without data in initial GWAS scan, we selected tag SNPs that were in LD with published SNPs (r 2 40.8).
yP value from meta-analysis of two cohorts in the initial discovery phase (logistic regression).
Differential expression analysis. To gain further insight into the possible involvement of genes at the identified loci in the development of NSCL/P, we examined the expression pattern of nine genes that were present within 50 kb of the identified SNPs or reported by previous studies by comparing microarray expression data from NSCL/P (n ¼ 7) and control (n ¼ 6) dental pulp stem cell cultures. We found that ADCY9 and NTN1 were significantly upregulated in the dental pulp stem cultures of NSCL/P patients (P Wilcoxon o0.05, Supplementary Table 6).

Discussion
In this GWAS, four recognized loci were confirmed (1q32.2, 10q25.3, 17p13.1 and 20q12) and one new locus were identified (16p13.3) in the Chinese populations. The newly identified SNP rs8049367 at 16p13.3 is located 50 kb upstream of CREBBP and 32 kb downstream of ADCY9 (Fig. 2a). CREBBP is a transcription factor that has been reported to play essential roles in embryonic development via co-activation with the protein P300 (refs 13,14). Several studies have reported a causal relationship between deletions in the CREBBP-ADCY9 region and Rubinstein-Taybi syndrome [15][16][17][18][19][20] , which is characterized by distinctive facial features, including a cleft lip and palate 21 . Expression analysis revealed that expression of ADCY9 but not CREBBP was significantly different in dental pulp stem cell cultures from NSCL/P patients compared with the dental pulp stem cell cultures from normal controls (P Wilcoxon ¼ 4.66 Â 10 À 3 and 6.28 Â 10 À 1 , respectively, Supplementary Table 6). Thus, ADCY9 represents a potential susceptibility gene for NSCL/P.  We further analysed the biological potential of the identified variant but did not observe any signals of enhancers or other regulatory elements in the regions nearby rs8049367 and its closely linked SNPs ( Supplementary Fig. 4b). Notably, we found that rs2262251, which was in strong LD with rs8049367 in people of Asian ancestry (r 2 ¼ 0.96), is located in the exon of a newly identified long non-coding RNA (lncRNA), RP11-462G12.2 (GENCODE version 19). Therefore, we predicted an interaction between the lncRNA and two nearby genes (CREBBP and ADCY9) and found that both lncRNA-gene pairs exhibited high  Table 7), which suggested potential interactions. Together, these data suggest that the genetic variants in the CREBBP-ADCY9 region influence the development of NSCL/P, possibly by regulating non-coding RNA. Additional tissue-specific expression quantitative trait loci studies are warranted to further illuminate the relationship between the variants and CREBBP-RP11-462G12.2-ADCY9. Before the GWAS era, IRF6 had been consistently identified as a susceptibility gene for nonsyndromic CL/P across studies 1 . To date, three SNPs (rs642961, rs861020 and rs10863790) near IRF6 were reported to be associated with NSCL/P in previous GWAS 8,12 , and rs642961 was found to be a causative variant 6 in people of European ancestry. In our study, six SNPs (rs2235371, rs1044516, rs742214, rs596731, rs674433 and rs12117075) were consistently associated with NSCL/P in all three stages (Supplementary Table 4). rs2235371 was in strong LD with rs10863790, and rs674433 was in strong LD with rs642961 and rs861020 (r 2 Chinese 40.8). The other SNPs in the region identified in our study revealed low to moderate LD with rs2235371 (r 2 Chinese ¼ 0.22-0.77; Fig. 2b; Supplementary Table 8). A conditional analysis of the combined all cohorts demonstrated that the effects of all identified SNPs were abolished or substantially reduced by rs2235371 (Supplementary Table 8), suggesting that these variants may represent the same causal variants. Interestingly, the SNP rs2235371, a missense variant with a functional Combined Annotation Dependent Depletion (CADD) score of 32, is located in a region with a regulatory element signal in hESH1-derived mesenchymal cells according to the NIH's Roadmap Epigenomics Project ( Supplementary Fig. 5). The CCCTC-binding factor (CTCF) Chromatin Interaction Analysis by Paired-End Tag Sequencing (ChIA-PET) data in the Encode Project revealed an interaction signal between regulatory regions and DNA sequences containing rs2235371 in multiple cell lines, indicating that rs2235371 may be involved in chromatin activities. Thus, the effect of the variant rs2235371 may be complicated, and additional studies are warranted to clarify the underlying mechanisms. SNP rs1880646 near NTN1 at 17p13.1 was originally reported in a GWAS by Beaty et al. 10,11 and in their follow-up study. In our study, rs1880646 was not associated with NSCL/P in the GWAS stage (P meta ¼ 0.76), suggesting that rs1880646 may not be an appropriate proxy of causal variants at 17p13.1 in the Chinese population. Rather, our newly identified SNP, rs4791774, in the intron of NTN1 was in strong LD with SNP rs4791331, identified by Beaty et al. in Asian populations (r 2 ¼ 1), but not with rs1880646 (r 2 ¼ 0.12) (Fig. 2c). Motif analyses by HaploReg indicated that the G allele of rs4791774 could greatly alter the affinity of MYC (score: 0.5-12.5). NTN1 encodes the protein NETRIN 1, which plays a role in the developing the nervous system by promoting both axonal outgrowth and axonal guidance in pathfinding [22][23][24][25] . In addition, NTN1 was upregulated in dental pulp stem cell cultures from NSCL/P patients, suggesting the importance of NTN1 in the development of NSCL/P.
One of the two reported SNPs, rs7078160 at 10q25.3, is located in the intron of KIAA1598 and is downstream of VAX1 (Fig. 2d). VAX1 encodes a transcriptional regulator with a DNA-binding homeobox domain. Mouse knockouts for VAX1 develop a cleft palate, and this gene is expressed widely in developing craniofacial structures 26 . Nonetheless, recent target sequencing in the region of the gene failed to identify causal variants 27 , indicating that further attention should be devoted to the regulatory region of VAX1. Another SNP, rs13041247 at 20q12, maps 45 kb downstream of MAFB (Fig. 2e), which encodes a basic leucine zipper transcription factor. In mice, MAFB expression was demonstrated to be high in the epithelium of the palatal shelves and in the medial edge epithelium during palatal fusion 10 .
Previous work has revealed evidence of heterogeneity between European and Asian populations in the identified susceptibility loci of NSCL/P. The SNP rs987525 at 8q24.21, which has recently been demonstrated to be located in the region near a remote MYC-regulated enhancer 28 , was not validated in either Beaty's Asian trios 10 or in our study; the derived allele is so rare in Asian populations that studies of it cannot achieve sufficient power. We screened 500-kb regions around the SNP in our imputation results; however, we did not find any SNPs with Po1.0 Â 10 À 5 (Supplementary Fig. 6b). A similar situation occurred at rs2235371 and rs10863790 near IRF6, although these SNPs are Asian specific. In addition, our study is the first to demonstrate that the effect of rs642961, which was thought to be a causal variant in European populations, can be weakened by rs2235371, suggesting that the key variants may be located on rs2235371 in related haplotypes. In contrast to the SNPs at 8q24.21 and 1q32.2, the newly identified SNP rs8049367 tended to be shared throughout the population and to undergo positive selection; the frequency of the C allele is 68% in East Asians in contrast with only 27% in Africans and 37% in Europeans (F st ¼ 0.113). We noted that the trend of the frequency of the C allele in people of Asian, European and African ancestries is consistent with the prevalence of NSCL/P 2 ; thus, rs8049367 may help to explain the variable prevalence of NSCL/P in populations of different ancestries.
In this study, five loci-1q32.2, 10q25.3, 16p13.3, 17p13.1 and 20q12-were consistently identified as being associated with NSCL/P and reached genome-wide significance after combined analysis, and 16p13.3 was newly identified as an NSCL/P susceptibility region. To the best of our knowledge, this is the first GWAS performed in the Chinese population and the largest in Asia. The results of our study provide new insight into the aetiology of NSCL/P and indicate potential targets for further clinical translational research. However, the power of the study was limited because of the cost-efficient multiple stages design. Joint analysis on published GWAS is warranted to identify more loci associated with NSCL/P.

Methods
Subjects. In this study, we performed a three-stage case-control analysis, including an initial discovery phase and two validation phases (Supplementary Table 1; Supplementary Fig. 1). The discovery phase comprised 858 NSCL/P cases and 1,248 controls derived from two independent cohorts: (i) the Huaxi cohort (504 cases and 455 newborn controls) and (ii) the Nanjing cohort (354 cases and 793 controls). The first validation comprised 1,319 NSCL/P cases and 1,526 newborn controls from three independent cohorts: (i) the Guangzhou cohort (Ia: 487 cases and 496 controls); (ii) the Huaxi cohort (Ib: 503 cases and 508 controls); and (iii) the Nanjing cohort (Ic: 329 cases and 522 controls). The second validation was derived from the Ningxia Hui population, a minority nationality in China, and included 344 cases and 348 newborn controls. All participating patients were interviewed and clinically assessed by an experienced oral surgeon based on detailed diagnostic information from medical records and physical examinations to ensure that individuals with syndromic orofacial clefts or other major congenital defects were excluded. In addition, general characteristics, including age, sex, ethnicity, and birth defects, were documented. At recruitment, informed written consent was obtained from each of the participants or their guardians, and this study was approved by the institutional review boards of general hospital of Ningxia Medical University, west China hospital of stomatology Sichuan University, Nanjing Medical University and Guangdong women and children hospital.
Genotyping and quality control in GWAS. All samples from the Huaxi cohort and NSCL/P cases from the Nanjing cohort were genotyped using Affymetrix Axiom Genome-Wide CHB1 and CHB2 arrays by the CapitalBio corporation (1,280,786 SNPs); the controls from the Nanjing cohort were from a previous study 29 and genotyped using an Affymetrix Genome-Wide Human SNP Array 6.0 (905,119 SNPs). Considering that there were only approximately half a million SNPs overlapping between the two platforms, we used imputed data as a control for the Nanjing cohort. We performed a systematic quality control analysis on the raw data in each cohort before conducting the association analysis (see the workflow diagram in Supplementary Fig. 2). First, we performed basic quality NATURE COMMUNICATIONS | DOI: 10.1038/ncomms7414 ARTICLE control on genotyping data. SNPs were excluded if they (i) did not map to autosomal chromosomes, (ii) had a call rate of o95%, (iii) had a minor allele frequency (MAF) of o0.05 or (iv) had a genotype distribution in the controls that deviated from the Hardy-Weinberg equilibrium (Po1.0 Â 10 À 5 ). Samples were excluded if they (i) had overall successful genotyping call rates o95%, (ii) had sex discrepancies between the records and the genetically inferred data or (iii) had outlying autosomal heterozygosity rates (46 s.d. from the mean). Following the basic quality control, we extracted best-guess genotype data for SNPs with imputation quality info 40.8 and MAF 40.05 ref. 30 of sex-matched individuals from previous studies 29 and combined them with the genotype data of the Nanjing cases. Next, we removed individuals who had unexpected duplicates or probable relatives (all PI_HAT 40.25). Finally, we detected population outliers using a method based on principal component analysis. Common autosomal SNPs (MAF40.25) with a low LD (r 2 o0.05) in each cohort were employed to identify population outliers in the samples that had passed quality control, with Han Chinese from Bejing (CHB, N ¼ 97), Southern Han Chinese (N ¼ 100) and Japanese from Tokyo (N ¼ 89) from the 1000 Genomes Project (the Phase I integrated variant set release, v3) as the internal controls ( Supplementary Fig. 3a). The Principal component analysis (PCA) suggested that the cases and controls were genetically matched and homogeneous (Supplementary Fig. 3) and that the genomic control inflation factor (l) was 1.03 in the Huaxi cohort and 1.10 in the Nanjing cohort. After quality control procedures had been performed, 858 cases and 1,248 controls (504 cases and 455 controls in the Huaxi cohort; 354 cases and 793 controls in the Nanjing cohort) with 842,556 overlapping SNPs were included in the final analysis.
SNP selection and genotyping in the replication study. We performed a metaanalysis on the results of the two GWAS cohorts, and the following criteria were used to select SNPs for validation from the discovery phase: (i) P meta r1.0 Â 10 À 5 , (ii) P Huaxi r0.01 and P Nanjing r0.01, (iii) clear genotyping clusters and (iv) the SNP with the lowest P value was selected when multiple SNPs were observed in strong LD (r 2 Z0.8, Supplementary Table 2). SNPs identified in previous GWAS with P meta o0.05 or highly correlated SNPs with P meta o0.05 in the discovery phase were also included (see Table 2). rs2235371 and rs674433 met our GWAS criteria and were in high LD with published SNPs; we validated them rather than rs10863790 and rs642961(rs861020), respectively. To better explore the causal variants at 1q32.2, we included rs10863790 in the validation. As a result, 29 SNPs remained in the first validation (Supplementary Table 4). Genotyping was performed using the iPLEX MassARRAY platform (Sequenom Inc.) for two-phase validations. The primers and probes are available upon request. The laboratory technicians who performed the genotyping experiments were blinded to the case or control status.
Functional annotation and different expression analysis. To investigate the underlying mechanisms of the lead SNPs identified in our studies, we integrated chromatin biofeature annotations with 1,000 Genomes genotyping data using the Bioconductor R package FunciSNP 31 . The following resources were employed to filter the correlated SNPs lying within putative regulatory elements: (1) Roadmap: Chromatin Primary Core Marks Segmentation by HMM from H1 cell lines and hESH1-derived mesenchymal cells; (2) Fantom: an active, in vivo-transcribed enhancers atlas 32 ; (3) Enhancer identified in a previously published paper 33 . Regional functional plot were generated by the UCSC genome browser (http:// genome.ucsc.edu/). Correlated SNPs (r 2 40.8) within a window size of 1 Mb around the newly identified SNP were used for a FunciSNP analysis. HaploReg V2 (http://www.broadinstitute.org/mammals/haploreg/haploreg.php) was used to predict regulatory motif alteration. lncPro (http://bioinfo.bjmu.edu.cn/lncpro/) was applied to predict interactions between lncRNA and nearby proteins. Processed microarray expression data from NSCL/P (n ¼ 7) and control (n ¼ 6) dental pulp stem cell cultures were downloaded from EMBL-EBI (E-GEOD-42589). Expression measurements were performed using the Affymetrix Human Gene 1.0 ST array, which interrogates 28,869 transcripts, followed by RNA labelling and hybridization protocols as recommended by the manufacturer 34 . A Wilcoxon rank sum test was applied to assess differences in expression levels.
Statistical analysis. We used PLINK 1.07 (ref. 35) and R 2.15.1 (http://www. r-project.org/) for general statistical analysis. The population structure was evaluated by PCA using the EIGENSTRAT 3.0 software package 36 . Associations between SNP genotypes and disease status were assessed in an additive model using logistic regression model (1 degree of freedom) with adjustment(s) for sex (and significant principle components in GWAS). A meta-analysis-based method was applied to combine the results from different cohorts. A fixed-effect model with inverse variance weighting was used when there was no indication of heterogeneity (P for Cochran's Q statistic 40.05); otherwise, a random-effect model for the corresponding SNPs was adopted. The genomic inflation factor for meta-analysis was estimated on SNPs without heterogeneity. A Manhattan plot of -log 10 P was generated using the ggplot2 package 37 in R 2.15.1. Non-genotyped SNPs were imputed in the GWAS discovery samples using Shapeit v2 (http://www.shapeit.fr/, phasing step) and IMPUTE2 (http://mathgen.stats.ox.ac.uk/impute/ impute_v2.html, imputation step) based on haplotype information from the 1,000 Genomes Project (http://www.1000genomes.org) (the Phase I integrated variant set release, v3, across 286 Asians). Regional plots were generated using LocusZoom 1.3 Standalone 38 . The odds ratios and 95% confidence intervals were calculated using an additive model in logistic regression analyses, except where otherwise specified.