article: Nonsyndromic orofacial clefts are common birth defects with multifactorial etiology. The most common type is cleft lip, which occurs with or without cleft palate (nsCLP and nsCLO, respectively). Although genetic components play an important role in nsCLP, the genetic factors that predispose to palate involvement are largely unknown. In this study, we carried out a meta-analysis on genetic and clinical data from three large cohorts and identified strong association between a region on chromosome 15q13 and nsCLP (P = 8. 13×10−14 for rs1258763; relative risk (RR): 1. 46,95% confidence interval (CI): 1. 32–1. 61) ) but not nsCLO (P = 0. 27; RR: 1. 09 (0. 94–1. 27) ). The 5 kb region of strongest association maps downstream of Gremlin-1 (GREM1), which encodes a secreted antagonist of the BMP4 pathway. We show during mouse embryogenesis, Grem1 is expressed in the developing lip and soft palate but not in the hard palate. This is consistent with genotype-phenotype correlations between rs1258763 and a specific nsCLP subphenotype, since a more than two-fold increase in risk was observed in patients displaying clefts of both the lip and soft palate but who had an intact hard palate (RR: 3. 76, CI: 1. 47–9. 61, Pdiff<0. 05). While we did not find lip or palate defects in Grem1-deficient mice, wild type embryonic palatal shelves developed divergent shapes when cultured in the presence of ectopic Grem1 protein (P = 0. 0014). The present study identified a non-coding region at 15q13 as the second, genome-wide significant locus specific for nsCLP, after 13q31. Moreover, our data suggest that the closely located GREM1 gene contributes to a rare clinical nsCLP entity. This entity specifically involves abnormalities of the lip and soft palate, which develop at different time-points and in separate anatomical regions.
  Nonsyndromic cleft lip with or without cleft palate (nsCL/P) is a common human birth defect with a multifactorial etiology, including a strong genetic component [1,2]. Previous studies have identified 16 genetic risk loci for nsCL/P. These studies comprised candidate gene and linkage analyses [3–5], genome-wide association studies (GWAS) with follow-up approaches [6–11], and a meta-analysis [12]. Despite these advances in deciphering the genetic architecture of nsCL/P, a number of additional risk loci still await identification. Some of these as yet unknown susceptibility variants may be detectable in GWAS datasets but have escaped detection at a genome-wide significant level due to low statistical power, which is secondary to limited sample sizes. This suggests that further risk variants for nsCL/P might be identified via one of the following approaches: the combining of available data sets, targeted replication analyses in independent cohorts, and/or the reduction of clinical heterogeneity using detailed subphenotype information. NsCL/P shows considerable phenotypic variability in terms of affected anatomical structures, and can be subdivided into two main forms: nonsyndromic cleft lip only (nsCLO) and clefts involving both the lip and the palate (nsCLP) [2]. This distinction is important in terms of the degree of physical handicap and treatment. Although epidemiological data indicate that these subtypes are determined at least in part by genetic predisposition [13], few data are available concerning the specific genetic factors determining the formation of nsCLP as opposed to nsCLO. To date, one locus (at 13q31) has shown a specific association with nsCLP but not with nsCLO [12,14], while IRF6 has shown a predominant effect in nsCLO [5]. Previous research has implicated the Gremlin-1 (GREM1) locus in human orofacial clefting. This research has involved the investigation of gene networks and, more recently, the finding of association between variants at 15q13 and human nsCL/P [10,15]. However, associations of common variants were not yet significant at the genome-wide level. Also, sequencing studies of GREM1 in nsCL/P patients and controls have been conducted in limited sample sizes only, with inconclusive results: although our group has previously generated some evidence for the role of rare variants in the GREM1 coding and untranslated region [16], the functional relevance of the identified variants remained unclear, and the results of burden analyses varied depending on the test applied. In another sequencing study, no deleterious rare variants were identified in GREM1 [15]. Analyses of Grem1-deficient mouse models have shown that during embryogenesis Grem1 function is crucial for limb development and kidney formation. However, complete loss of Grem1 function causes no obvious craniofacial defects [17,18]. GREM1 acts as a secreted antagonist of various members of the bone morphogenetic protein (BMP) family, which has been shown to play a critical role in both lip and palate development [19,20]. Notably, previous research has indicated a particular role for BMP4, which is involved in facial genesis [21,22]. Moreover, rare mutations within BMP4 have been associated with human clefting [23], and it is established that soluble GREM1 binds with high affinity to BMP4 [24]. Loss-of-function and overexpression studies of the related Bmp antagonist Noggin in mice have previously demonstrated the critical role of restricted Bmp signaling during lip and palate development, as well as for midfacial morphogenesis [25–27]. The present study analyzed the putative risk locus on human chromosome 15q13 in different nsCL/P cohorts, using detailed clinical subphenotype information. The genetic analyses were complemented by functional analyses of the mouse ortholog of the GREM1 gene, which is located adjacent to the associated region. Finally, existing comprehensive genomics data sets were analyzed to annotate the associated region and to establish GREM1 as a plausible candidate gene for nsCL/P at 15q13.
  Regional association statistics for 219 variants at 15q13 (chr15: 32. 95–33. 5 Mb, hg19) were extracted from a previously published meta-analysis (referred to as Ludwig 2012 meta-analysis) [12] (Table 1, S1 Dataset). The top associated variant at 15q13 in this Ludwig 2012 meta-analysis data set was rs1258763 (PnsCL/P_meta_all = 1. 81×10−06, Fig 1 and S1 Dataset), with the A-allele representing the risk allele. As part of the present study, rs1258763 was then genotyped in an independent case-control cohort of mixed ethnicity (replication I, Table 1 and S2 Dataset). After quality control, 580 cases and 1,684 controls remained in the analysis. In this replication I data, rs1258763 showed strong association with nsCL/P (PnsCL/P_rep_I = 4. 34×10−09, relative risk (RR) for the major allele A: 1. 58,95% confidence interval (95% CI): 1. 36–1. 84, S2 Dataset). Additionally, rs1258763 had been previously genotyped in an independent trio sample from the EuroCran cohort [10] (replication II, Table 1). Re-analysis of this data (to exclude individuals that were also part of replication I) revealed PnsCL/P_rep_II = 0. 072 (S2 Dataset), with the A-allele being overtransmitted to affected children. Combined analysis of the Ludwig 2012 meta-analysis and both replication I and II data sets, respectively, totaling 979 nsCL/P cases, 3,002 controls, and 2,061 trios, yielded genome-wide significance (PnsCL/P_comb = 2. 23×10−13), with a strong genetic risk observed for the major allele A (1. 35 (1. 24–1. 46), S2 Dataset). Imputation analysis of the 15q13 region in the Central European cohort that was part of the Ludwig 2012 meta-analysis (see Methods) revealed a 5 kb region of strongest association (Chr15: 33. 050–33. 055 Mb), located intergenically between GREM1 and Formin-1 (FMN1), with rs2600520 being the most strongly associated variant in the imputed data (PnsCL/P_imp = 5. 05×10−07, Fig 2A and S3 Dataset). To determine associations between variants at 15q13 and particular clefting subphenotypes, individuals from the three study cohorts were classified as nsCLP or nsCLO on the basis of available clinical information (S2 Dataset). Variant rs1258763 showed significant association with nsCLP in the genotyped data of the Ludwig 2012 meta-analysis (PnsCLP_meta = 8. 0×10−08), replication I (PnsCLP_rep_I = 1. 28×10−07) and replication II (PnsCLP_rep_II = 0. 035 Fig 3 and S2 Dataset). No significant association was detected for nsCLO (P>0. 2 in all cohorts, S2 Dataset). For each of the three cohorts, the difference in relative risk between nsCLP and nsCLO was statistically significant (P<0. 05). The overall association P-value for nsCLP in the combined analysis was PnsCLP_comb = 8. 13×10−14 (1. 46 (1. 32–1. 61), Fig 3), in contrast to PnsCLO_comb = 0. 27 (1. 09 (0. 94–1. 27) ) for nsCLO. The difference in relative risks between nsCLO and nsCLP in the combined sample was statistically significant in both multiplicative (P = 0. 0015) and general (P = 0. 0033) model. In the subphenotype analysis of the imputed data generated from the Central European case-control cohort rs2600520 remained the most significantly associated variant with nsCLP (PnsCLP_imp = 3. 19×10−06, Fig 2B), while only marginal associations were observed with nsCLO (Fig 2C). For the 5 kb region of strongest association, no compelling evidence for the presence of regulatory elements was found in ENCODE [32], except for a 200 bp region identified as DNAse hypersensitivity site in four cell types, three of which had been derived from skin tissue (chr15: 33,052,666–33,052,875, hg19) [33]. In datasets of relevance to craniofacial development [34,35], no active regulatory elements were detected. According to HaploRegv3 [36], the top associated variant rs2600520 and a second single nucleotide polymorphism (SNP), rs2600519, which is in high linkage disequilibrium (LD) and located four base pairs away, modify a chromatin mark that involves the transcriptional repressor Sin3A. The analysis of blood-sample expression quantitative trait loci (eQTL) yielded 38 SNPs with cis-eQTL effects on GREM1 expression, six of which had false discovery rates below 0. 05 (S1 Table). The strongest eQTL effect was observed for rs17816375, which showed a P-value of 0. 033 in the imputed nsCLP data. Of the 38 SNPs, the SNP with strongest association to nsCLP was rs16958561 (P = 2. 82×10−05). This SNP maps around 6 kb distal to rs2600520, and the two SNPs are in substantial LD (D' = 1 in the European population (1000 Genomes phase 1) ). No trans- or cis-eQTL effect was identified for FMN1, despite the fact that one FMN1-targeting probe was represented on the expression array used [37]. Analysis of chromatin interaction data [38] obtained in human epidermal keratinocytes (NHEK), umbilical vein endothelial cells (HUVEC) and mammary epithelial cells (HMEC) revealed a topologically associated domain (TAD) encompassing the entire GREM1 gene, the associated region and parts of the SCG5 and FMN1 coding regions. Standard annotation of contact domains reveals a subdomain separating the 5' region of FMN1 from its 3' end and GREM1 (S1 Fig). Further analysis of the TAD data using virtual 4C reveals evidence that direct interaction loops are present between the associated region and the GREM1 transcription start site at the given resolution (S2 Fig). X-Gal staining of the craniofacial region of heterozygous mutants of a mouse model of Grem1 deficiency [17] revealed Grem1 expression in the proximal region of both lateral and medial nasal prominences at E11. 5 and E12. 5 (Fig 4A and 4B). In addition, bilateral Grem1 expression was detected in the merging zones of the maxillary and medial nasal prominences during lip development at E12. 5 (Fig 4B). Grem1 was also expressed in the mesenchyme of the developing secondary palate between E12. 5 and E15. 5, with expression being restricted to the posterior palatal shelf region in which the soft palate forms (Fig 4C–4G). Although no craniofacial defects have been reported for Grem1-deficient mice [17], we tested the hypothesis that loss of Grem1 results in alterations of lip and palate morphology without manifesting any profound craniofacial phenotype. However, a histological analysis revealed no developmental abnormalities in the embryonic secondary palate or during lip formation in Grem1-/- mouse mutants (Fig 5A–5D). Next, E13. 5 secondary palate shelves were cultured in the presence of exogenous Grem1 protein to assess whether this treatment could affect their development. Significant growth (approximately 25% increase) of the palatal shelves was observed in both Grem1-treated and control groups during the 48 hour culture period, with no significant difference in the final size of the shelves observed between groups (P = 0. 27, S2 Table). However, Grem1-treated palatal shelves showed a more rounded shape, with retraction of their medial edges (Fig 5E–5G). Measurement of the area between both palatal shelves at 0 and 48 hours revealed a stable area in the control group but a significant increase in the Grem1-treated group (P = 0. 0014 for difference between control and Grem1-treated pairs; Fig 5H, S3 Table). Informed by the murine expression pattern of Grem1, we tested for involvement of the soft palate in human data sets. Twenty-one patients with a cleft lip and a cleft of the soft palate in the presence of an intact hard palate (nsCLPsoft) were selected (see Methods). In this group, analysis of rs1258763 revealed P = 0. 03 and an RR of 3. 76 (95% CI: 1. 47–9. 61) in nsCLPsoft, representing a two-fold increase of effect size in comparison to 320 patients with complete cleft lip and palate (nsCLPhard+soft, Fig 3 and Table 2). While only a limited number of patients with this rare nsCLPsoft phenotype could be recruited this association P-value passed the statistical significance threshold after permutation-based correction (P = 0. 04). No association was found in a group of 45 patients with cleft of the soft palate only and an intact lip (PnsCPO_soft = 0. 94, Table 2), or in 115 patients with a submucous cleft of the soft palate (PnsCPO_submuc = 0. 85, Table 2).
  The present study established that the 15q13 region contains a genetic risk factor for nsCL/P. This is the second locus after 13q31 to show a particularly strong association with nsCLP and not with nsCLO. Moreover, our results suggest this risk factor to be involved in nsCLPsoft, and our genetic and functional analyses provide strong support for the hypothesis that GREM1 is an effector gene that contributes to this rare clinical subphenotype. In aggregate, our data suggest that GREM1 plays a specific role not only in the development of the lip, but also during formation of the soft palate. These findings provide a framework for further functional analyses of GREM1 in human cell systems and for inducible modulation of Grem1 expression during the development of specific craniofacial regions in model organisms. These analyses may broaden our understanding of the processes that regulate facial morphogenesis, and help to decipher the molecular mechanisms underlying the manifestation of specific nsCL/P subphenotypes in humans.
  Genotyping of rs1258763 in the nsCL/P replication sample and nsCPOsoft individuals was performed using the Sequenom MassArray system (Agena Bioscience, San Diego, USA). nsCLPsoft-patients from the Dutch cohort and the nsCPOsubmuc patients were genotyped using ABI3130XL sequencing and BigDye v3. 1. The replication II sample was analyzed using the transmission disequilibrium test. In the replication I cohort, genotype frequencies in the cases and controls of each of the three subsamples (Bonn, Mexico, Yemen) were compared using the Cochran-Armitage trend test. To combine all results for rs1258763, effect estimates for the different studies were combined in an inverse variance-weighted fixed-effects meta-analysis. To test for heterogeneity of the genotypic RRs between the nsCLO and nsCLP phenotypes, a heterogeneity likelihood-ratio test was applied using a general and a multiplicative model. For each of these analyses, this resulted in an asymptotic Chi2 null distribution with two degrees of freedom. For the analysis of rs1258763 in the nsCL/P subphenotypes nsCLPhard+soft, nsCLPsoft, nsCPOsoft, and nsCPOsubmuc, a common set of 1,318 controls from the Central European case-control cohort was used. To determine whether a statistically significant difference was present between the results for nsCLPhard+soft and nsCLPsoft and to account for the limited number of individuals in the nsCLPsoft group, an empirical P-value was determined using 100,000 permutations (larger test statistics were assigned a value of 1, equal test statistics were assigned a value of 0. 5). Genotypes from the Central European case-control cohort that was part of the meta-analysis [10] were imputed using IMPUTE2 [59] and 2,184 alleles of the 1000genomes project. Dosage values were included in a logistic regression model (using SNPTEST and—method expected) in which we included as covariates the first five components eik (which were obtained from MDS analysis [60], to account for population stratification. Only variants showing a SNPTEST info-score (detailed description: https: //mathgen. stats. ox. ac. uk/genetics_software/snptest/snptest. v2. pdf) of at least 0. 4 and a minor allele frequency of at least 1% were retained, as accuracy of imputation is compromised for low frequency variants. The resolution of the 15q13 region (250 kb, chr15: 32. 95 Mb—33. 2 Mb) increased from 106 variants in the genotyped data to 1,042 variants in the imputed data. Relative risks were calculated using the beta values, representing the logarithm of the RR. Please note that inclusion of the genome-wide data from Beaty et al. 2010 in the imputation analysis was not possible at the time of the present study, due to insufficient coverage of that approach in the dbGaP approval used for the Ludwig 2012 meta-analysis. Potential regulatory effects of the top associated region on either GREM1 or FMN1 were assessed using previously published eQTL and HiC-datasets. Since there is a complete absence of eQTL studies of fetal craniofacial samples, eQTLs obtained in large datasets from other tissues were used. Here we used a recently published and comprehensive study of blood samples by Westra et al. [37]. Both FMN1 and GREM1 probes are present on the Illumina arrays used in the Westra et al. study. To identify altered transcription factor binding we used the v3 version of the HaploReg tool, which was developed to annotate disease-associated genetic variants located in non-coding regions [36]. In particular, information on transcription factor binding altered by nsCL/P risk alleles was used. These analyses were complemented by recently published HiC-data illustrating chromatin formation and loops at 15q13 using the Juicebox tool [38] and virtual 4C as provided at http: //promoter. bx. psu. edu/hi-c/virtual4c. php using the same data. Grem1LacZ mutant mice [17] were kept on a C57Bl/6 genetic background and embryos were staged using mid-day on the day of vaginal plug detection as embryonic day 0. 5 (E0. 5). Genotyping of embryos was carried out using allele-specific PCR. The wild type allele was detected as a 435bp product using primers 5’-TGCAATTGTGTCAGGAGCCA-3’ and 5’-ACTGGGTCTGCTCAGAGTCA-3’. The lacZ allele was detected as an approximately 550bp product using primers 5’-TGCAATTGTGTCAGGAGCCA-3’ and 5’-GGGAACAAACGGATTGACCG-3’. Following embryo dissection, mandibles and tongues were removed to facilitate the penetration of fixatives and staining solutions. Expression of Grem1 was visualized by X-Gal staining of heterozygous Grem1LacZ whole mount embryos using standard procedures, as described elsewhere [61]. To assess tissue-specific domains of Grem1 expression, 15μm cryosections of whole mount stained embryos were prepared and mounted onto glass slides. Images were obtained using an Axiocam color camera mounted on a Zeiss Stemi SV11. Images were processed using the softwares Axiovision AC (release 4. 4) and Photoshop C4 (version 11). Organ culture experiments of secondary palatal shelves were carried out using conditions essentially similar to those described elsewhere [61,62]. Briefly, palatal shelves were dissected at E13. 5 and cultured for 48 hours on sterilized filter papers (Millipore) placed on a metal grid located on the central well of a culture dish. With the aboral side facing downwards, the palatal shelves were cultured in Dulbecco’s modified eagles medium (DMEM), supplemented with 10% (v/v) fetal bovine serum, 1% (v/v) penicillin/streptomycin (5000 units/ml, 5000ug/ml), 1% glutamine (200mg/ml) and 1% (v/v) ascorbic acid (50mg/ml). Grem1 recombinant protein (Peprotech (USA) ) was used at a final concentration of 10μg/ml medium. The organ rudiments were cultured at 37°C using 5% CO2 for 48 hours and culture medium was replaced after 24 hours. Images were obtained using a Axiocam mounted onto a Zeiss Stemi SV6 stereomicroscope. Areas between palatal shelves were measured using imagej (http: //imagej. nih. gov/ij/). To calculate P-values, arbitrary units quantifying the area between the palatal shelves were analyzed using a one-tailed t-test. For the craniofacial phenotype analysis of Grem1-deficient mouse embryos, processing of heads for embedding in paraffin and generation of Hematoxilin/Eosin stained sections were performed as described elsewhere [61]. For the purposes of the present study, the following databases were accessed: https: //www. broadinstitute. org/mpg/snap/ [63], http: //genenetwork. nl/bloodeqtlbrowser/ [37], http: //www. emouseatlas. org/gxdb/dbImage/segment1/2906/2906. html [41], http: //www. aidenlab. org/juicebox/ [38].
summary: Clefts of the lip and palate are common birth defects, and require long-term multidisciplinary management. Their etiology involves genetic factors and environmental influences and/or a combination of both, however, these interactions are poorly defined. Moreover, although clefts of the lip may or may not involve the palate, the determinants predisposing to specific subphenotypes are largely unknown. Here we demonstrate that variations in the non-coding region near the GREM1 gene show a highly significant association with a particular phenotype in which cleft lip and cleft palate co-occur (nsCLP; P = 8. 13×10−14). Our data suggest that the risk is even higher for patients who have a cleft lip and a cleft of the soft palate, but not of the hard palate. Interestingly, this subphenotype corresponds to the expression of the mouse Grem1 gene, which is found in the developing lip and soft palate but not in the hard palate. While Grem1-deficient mice display no lip or palate defects, we demonstrate that ectopic Grem1 protein alters palatal shelve morphogenesis. Together, our results identify a region near GREM1 as the second, genome-wide significant risk locus for nsCLP, and suggest that deregulated GREM1 expression during craniofacial development may contribute to this common birth defect.
keywords: medicine and health sciences
  cleft lip
  palate
  developmental biology
  mathematics
  statistics (mathematics)
  cleft palate
  mammalian genomics
  embryos
  morphogenesis
  digestive system
  research and analysis methods
  embryology
  cleft lip and palate
  birth defects
  mathematical and statistical techniques
  otorhinolaryngology
  congenital disorders
  genetic loci
  animal genomics
  mouth
  anatomy
  meta-analysis
  genetics
  biology and life sciences
  physical sciences
  genomics
  statistical methods
title: Meta-analysis Reveals Genome-Wide Significance at 15q13 for Nonsyndromic Clefting of Both the Lip and the Palate, and Functional Analyses Implicate GREM1 As a Plausible Causative Gene