Allelic and genotypic association analysis on 18 candidate genes
Of the 18 candidate genes associated with muscle development, 399 SNPs were successfully genotyped. The results of HWE Tests and MAF of these SNPs are listed in Supplementary Table 1. HWE test was used to ensure the representation of each SNP in the population and 344 qualified SNPs with HWE p-value more than 0.05 were included in discovery phase and performed allelic and genotypic association analysis.
Allelic association analysis was aimed to assess the association between alleles on candidate SNPs and NSOC subtypes, including non-syndromic cleft lip only (NSCLO), non-syndromic cleft lip and palate (NSCLP) and non-syndromic cleft palate only (NSCPO) and the first two subtypes can be combined as non-syndromic cleft lip with or without cleft palate (NSCL/P). NSCLO was furtherly divided into bilateral cleft lip (BCL) and unilateral cleft lip (UCL) subgroups, including left cleft lip (LCL) and right cleft lip (RCL), and performed allelic association analysis accordingly. SNPs with p-value<0.05 were listed in Supplementary Table 2 with detailed information of the results, including p-value, odds ratio (OR) and 95% confidence interval (CI). After Bonferroni correction, the p-value threshold was set as 0.0001453 (0.05/344 SNPs).
Results indicated that all significant SNPs passed Bonferroni correction were located at ALX4 gene (Supplementary Table 2). Allele C of rs11037948 is associated with NSCPO (p = 8.181E-07, OR = 0.6203, 95%CI: 0.513–0.7499), BCL (p = 1.874E-07, OR = 0.5981, 95%CI: 0.493–0.7257) and UCL (p = 0.00001108, OR = 0.6732, 95%CI: 0.5643–0.8031), including RCL (p = 7.043E-08, OR = 0.5933, 95%CI: 0.4908–0.7174) and LCL (p = 0.00001336, OR = 0.6696, 95%CI: 0.559–0.8021). Allele T of rs3861063 is associated with BCL (p = 3.083E-07, OR = 0.6116, 95%CI: 0.5066–0.7383) and UCL (p = 4.825E-07, OR = 0.6397, 95%CI: 0.5376–0.7613), including RCL (p = 3.759E-07, OR = 0.6195, 95%CI: 0.515–0.7452) and LCL (p = 3.361E-07, OR = 0.6296, 95%CI: 0.5271–0.7521). Besides, allele C of rs10838258 is associated with UCL (p = 3.134E-09, OR = 0.5679, 95%CI: 0.471–0.6848), more specifically, with RCL (p = 1.07E-12, OR = 0.4781, 95%CI: 0.3902–0.5857).
Genotypic association analysis was conducted to assess whether specific genotypes of SNPs could influence NSOC risk. The results showed that rs17447140 on GATA3 gene was associated with NSCLP (p = 9.46E-07), NSCPO (p = 9.46E-07) and every subgroup of NSCL, including BCL (p = 4.05E-05), LCL (p = 7.40E-05), RCL (p = 2.90E-05) and UCL (p = 6.57E-05). Six SNPs on ETS1 showed significant association with NSCLP, including rs4245078 (p = 2.53E-05), rs4245079 (p = 9.49E-05), rs4262739 (p = 7.57E-05), rs4245080 (p = 0.0001302), rs7115781 (p = 7.54E-05) and rs7130041 (p = 1.03E-04). Every other SNP passed Bonferroni correction was located at ALX4 gene and a total of 7 SNPS reached the significant p-value in at least one different subtypes of NSOC or sub-phenotypes of NSCLO, including rs3861063, rs11037948, rs10838258, rs10742700, rs879238, rs7102098 and rs10838257 (Supplementary Table 3).
Replication association analysis on SNPs located at ALX4
To further verify the association between ALX4 and NSOC, we chose three SNPs (rs12420503, rs879239 and rs3861063) on ALX4 gene with HWE p-value > 0.05 (Supplementary Table 4) to perform replication analysis in another group of case-control subjects. Allelic association results between SNPs and subtypes of NSOC were showed in Fig. 1 and detailed information including p-value, odds ratio and 95% CI were listed in Supplementary Table 5. Allele C of rs879239 was associated with unilateral cleft lip with or without cleft palate (UCL/P) (p = 0.04042, OR = 1.12, 95%CI: 1.005–1.249), left cleft lip with or without cleft palate (LCL/P) (p = 0.008197, OR = 1.183, 95%CI: 1.044–1.341), UCL (p = 0.03863, OR = 1.132, 95%CI: 1.006–1.274) and LCL (p = 0.02445, OR = 1.168, 95%CI: 1.02–1.338). Allele T of rs3861063 is associated with microform cleft lip with or without cleft palate (MCL/P) (p = 0.03624, OR = 0.7902, 95%CI: 0.6338–0.9854) and microform cleft lip (MCL) (p = 0.02814, OR = 0.7759, 95%CI: 0.6183–0.9736). Allele C of rs10838259 is only significantly associated with BCL (p = 0.03624, OR = 1.355, 95%CI: 1.019–1.801).
The results of genotypic association analysis were presented in Fig. 2 and Supplementary Table 6. rs10838259 is associated with NSCPO (p = 0.03705), bilateral cleft lip with or without cleft palate (BCL/P) (p = 0.03263) and BCL (p = 0.006806). Besides, rs879239 is associated with LCL/P (p = 0.02925) and rs12420503 is associated with complete cleft lip (CCL) (p = 0.01023).
ALX4 was differentially up-regulated in NSCLO patients
We collected six MSCLO patients as case group and two healthy subjects with lip trauma as control group and performed RNA sequencing using their lip tissues. ALX4 was differentially up-regulated in NSCLO case group compared with control group. There were 11 differentially expressed miRNAs, of which 2 were up-regulated, including has-miR-455-3p. Three databases, including TargetScan, RNAhybrid and miRanda, predicted ALX4 as the target gene of has-miR-455-3p.
To validate the accuracy of the sequencing data, we performed RT-qPCR in NSCLO cases and controls and verified the up-regulation of ALX4 and has-miR-455-3p (Fig. 3). Dual-luciferase reporter assay results showed that compared with transfection of wild-type ALX4-3’ untranslated region (UTR) alone, the co-transfection of has-miR-455-3p with wild-type ALX4-3’ UTR significantly reduced the luciferase activity. The luciferase activity showed no significant change with the co-transfection of has-miR-455-3P with mutant ALX4-3’ UTR (Fig. 4). This indicated that ALX4 was the target gene of has-miR-455-3P.