A functional polymorphism in the pre‐miR‐146a gene is associated with the risk of nonsyndromic orofacial cleft

microRNAs (miRNAs) are widely involved in craniofacial development, and genetic variants of miRNAs may be associated with the risk of nonsyndromic orofacial cleft (NSOC). Here, we systematically selected five single nucleotide polymorphisms (SNPs) of miRNAs and investigated the associations between these variants and NSOC susceptibility in a two‐stage case–control study including 1,406 NSOC patients and 1,578 controls from the Chinese population. We found that compared with the C allele, the rs2910164 G allele of pre‐miR‐146a was associated with an increased risk of NSOC (additive model: odds ratio [OR] = 1.17, 95% confidence interval [CI]: 1.06–1.30, P = 0.002), including both cleft lip with or without cleft palate (CL/P) and cleft palate only (CPO). Bioinformatic prediction and functional assays revealed that the C allele of rs2910164 was significantly associated with inhibited HEK‐293 and HEPM cell proliferation and decreased abundance of TRAF6. Both miR‐146a and TRAF6 were expressed in the lip tissue samples of NSOC patients, and a moderate inverse correlation was observed between them. Taken together, these results demonstrated that miR‐146a/rs2910164 is associated with susceptibility to NSOC, providing novel insights into the genetic etiology and underlying biology of NSOC.


INTRODUCTION
microRNAs (miRNAs), which are approximately 22 nucleotides in length, are important regulators of gene expression at the transcriptional, epigenetic, and posttranscriptional levels and influence processes such as embryonic development and organogenesis as well as human diseases (Ambros, 2004;Schoen et al., 2017). miR-NAs are widely considered to interact with their mRNA targets via Watson-Crick base pairing, and the functional importance of miRNA targeting events has been widely investigated (Kim et al., 2016).
Accumulating evidence has demonstrated the importance of genetic variations in miRNAs and genes in human diseases (Nicoloso et al., 2010;Ryan, Robles, & Harris, 2010). Single nucleotide polymorphisms (SNPs) in miRNA genes have the potential to affect the biogenesis of miRNAs and/or alter their target specificity (Ha & Kim, 2014).
Examples include miR-196a2 rs11614913 and miR-27a rs11671784, which influence the involvement of these miRNAs in diseases. The SNP miR-196a2 rs11614913 affects the processing of the miRNA precursor to its mature form and reduces its capacity to regulate target genes. miR-196a2 rs11614913 is associated with the risk of breast cancer (Hoffman et al., 2009). Genetic variation of miR-27a (rs11671784) decreases mature miR-27a expression, reduces gastric cancer susceptibility, and plays a role in gastric tumorigenesis (Yang et al., 2014).
Genetic variants of miRNAs have also been shown to be associated with the development of orofacial clefts. The SNP rs7205289, located in pre-miR-140, which is a key modulator of the platelet-derived growth factor (PDGF) signaling pathway during palatogenesis, contributes to nonsyndromic cleft palate susceptibility by influencing the processing of miR-140 (Li, Meng, Jia, Zhu, & Shi, 2010). Nevertheless, the role of SNPs within miRNA genes during craniofacial development has yet to be fully characterized.
In the current study, we exploited bioinformatics data from our previously published study with a validated database , to test the associations between genetic polymorphisms in miRNA genes and the risk of NSOC as well as to elucidate the underlying mechanism.

Human subjects
In the present study, we carried out a two-stage case-control study consisting of 1,406 NSOC patients and 1,578 healthy controls recruited from three affiliated hospitals of Nanjing Medical University since August 2008 (Ma et al., 2014;Yin et al., 2017

SNPs selection and genotyping
We selected the miRNA SNPs based on our previously published criteria . All available human miRNAs in the public miRBase database (version 10.0) were selected. We then performed extensive searches of SNP databases including HapMap (release #24), dbSNP, and Patrocles based on two criteria: (a) SNPs located in the pre-miRNA or mature sequences and (b) minor allele frequency (MAF) ≥5% in the Chinese population. Moreover, the Gibbs binding free energy (ΔG, KJ/mol) for each pair of common and variant alleles was computed using miRanda software. The difference in free energy between the two alleles was calculated as ΔΔG. As the energy parameter ΔΔG may affect the mature miRNA products (Gong et al., 2012) information about the subjects' case/control status. Approximately 5% of the samples were randomly selected to repeat the analysis.

Predicting the potential binding targets of miRNAs
Potential binding targets of miRNAs were predicted according to the following criteria: (1) predicted by three databases (Targetscan, miRDB, miRTarBase) and (2) had a high binding ability score (≥0.95).
Finally, TRAF6 was selected as the target of miR-146a. As a key adaptor molecule in the TLR or cytokine receptor/NF-B signaling pathway, TRAF6 is a direct target of miR-146a through a negative feedback regulation loop (Taganov, Boldin, Chang, & Baltimore, 2006). Iwata et al. reported that the TRAF6/TAK1/p38 signaling pathway, which is activated by TGF-, may regulate craniofacial development, including cleft palate (Iwata et al., 2012). Further identification of the target gene was complemented by its potential role in craniofacial malformations.

Construction of clones
To evaluate the effect of different alleles of rs2910164 on the expression of miR-146a, we amplified a 340-bp DNA fragment from a genomic DNA sample carrying the rs2910164 G allele and cloned it into pEGFP-N3 expression plasmid containing a green fluorescent protein marker (Clontech, Palo Alto, CA, USA). Then, a plasmid containing the C allele was generated from the construct containing the G allele by sitedirected mutagenesis.

Transfection
To assess the effect of rs2910164 on mature miR-146a expression, the HEK-293 and HEPM cells were seeded into 24-well culture plates at 3 × 10 5 cells per wells, and each well was transiently transfected with vector DNA containing either the rs2910164 G or C allele by Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. The pEGFP-N3 vector without an insert was used TA B L E 3 Associations between rs2910164 at miR-146a and risk for NSOC in overall and stratified analyses

Construction of reporter plasmids and dual-luciferase reporter assay
To generate luciferase reporter plasmids, the TRAF6 3 ′ -UTR fragment was inserted at the XbaI site downstream of the luciferase gene in the pGL3-promoter vector (Promega, Madison, WI, USA). All constructs were verified through DNA sequencing.
For the luciferase assay, HEK-293 and HEPM cells were seeded at 3 × 10 5 cells per well in 24-well culture plates and transiently cotransfected with the vector DNA containing either the rs2910164 G or C allele and the reporter plasmids containing the TRAF6 3 ′ -UTR fragment. The luciferase activities of the cellular extracts were quantified 48 hr posttransfection with a dual-luciferase reporter assay system (Promega). The ratio of Firefly luciferase to Renilla luciferase activity was evaluated. Three independent transfections were performed to ensure the reproducibility of results.

Western blot analysis
The total protein concentrations of the lysates were estimated using a BioRad protein assay kit based on the Bradford method. The total protein (20 g) in each sample was separated using 7.5% SDS-PAGE gels and transferred to a polyvinylidene fluoride (PVDF) membrane (Merck Millipore). The membranes were blocked with 5% nonfat milk in PBS/0.05% Tween 20 (PBS-T) for 2 hr at room temperature, incubated overnight at 4 • C, and probed with diluted primary antibodies against TRAF6 (Abcam, ab33915, 1:5000) and -actin (Abcam, 1:1000). Immunoblotted proteins were detected using the Phototope-

HRP Western Blot Detection System (Cell Signaling Technology).
Quantification analyses were carried out on Image J and the intensity of the TRAF6 protein bands was analyzed through densitometry after normalization to the corresponding -actin level.

Experimental animals
Five adult male and ten adult female C57BL/6 mice, purchased from Animal Center of Yangzhou University, were housed in a controlled temperature and a 12 hr/12 hr light/dark cycle environment. Two females and one male each were put together in one cage at 8:00 PM and separated at 8:00 AM the next morning. Embryos were considered as E0.5 on the morning of vaginal plug discovery. Later, E11.5, E12.5, E13.5, E14.5, and E15.5 embryos were collected for histological studies. All procedures were performed in compliance with federal regulations and the protocol followed was approved by the Institutional Animal Care and Use Committee.
Tissue sections were incubated at 4 • C overnight with the antibody.

Statistical and bioinformatics analyses
All statistical analyses were two sided and were performed using PLINK 1.07 software. The chi-square ( 2)   In the combined logistic regression analysis (Table 3), miR-146a/rs2910164 G allele was found to contribute to susceptibility to NSOC in the dominant (OR = 1.21, 95% CI: 1.04-1.41), recessive (OR = 1.27, 95% CI: 1.06-1.53), and additive models (OR = 1.17, 95% CI: 1.06-1.30). Furthermore, stratification analysis revealed that the rs2910164 G allele was a risk factor for CL/P and CPO under an additive model.

miR-146a/rs2910164 is associated with the expression of miR-146a in an allele-specific manner
The structural change in pre-miR-146a was predicted via an online database (https://www.bioguo.org/miRNASNP/) (Supporting Information Figure 1). The rs2910164 G allele was associated with significantly decreased miR-146a expression compared with the level of expression of the C allele in both cell lines (P < 0.001, Figure 1).

miR-146a/rs2910164 exerts different effects on cell proliferation
To explore the impact of miR-146a rs2910164 in the cellular environment, we transfected HEK-293 and HEPM cells with vectors containing different rs2910164 alleles and performed CCK8 assays to evaluate their effects on cell proliferation. A time-dependent increase in proliferation was detected in the G allele from 48 hr to 72 hr F I G U R E 1 Expression of miR-146a in human cell lines. HEK-293 and HEPM cells were transfected to express either the miR-146a G/C allele or control vector. The miR-146a expression levels were estimated by qRT-PCR; n = 3 (colonies) for each group and results are shown as mean values with the SD normalized to GAPDH. *P < 0.05, **P < 0.01, ***P < 0.001.
( Figures 2a and 2b). Flow cytometric analysis demonstrated that the rs2910164 G allele induced a substantial reduction in HEK-293 and HEPM cells in the G1 phase (P < 0.05 in HEK-293 and HEPM cells, Figures 2c and 2d), accompanied by a substantial increase in the number of cells in S phase (P < 0.05 in HEPM cells, Figures 2c and 2d).

TRAF6 might be the target gene of miR-146a
TRAF6 was predicted to be the target gene of miR-146a by three bioinformatics algorithms (Figure 3a). Since the miR-146a/rs2910164 C to G change is associated with a different level of mature miR-146a, it is possible that this variant regulates TRAF6 expression.
Hence, we generated luciferase constructs containing either the miR-146a rs2910164 G or C allele and the TRAF6 3 ′ -UTR plasmid to confirm the bioinformatic prediction. The reporter activity of the G allele was significantly higher compared with the activity of the C allele in the HEK-293 and HEPM cell lines (P < 0.05, Figure 3b). Furthermore, when transfecting miR-146a rs2910164 G/C into these cell lines, the rs2910164 G allele was associated with higher abundance of TRAF6 at the mRNA and protein levels compared with the C allele (P < 0.05, Figure 3c-e).

miR-146a and TRAF6 expression levels in human lip tissues
Both miR-146a and TRAF6 mRNA were detected in NSOC cleft lip tissue samples from 72 patients who underwent surgery. Furthermore, a moderate reverse correlation between miR-146a and TRAF6 was observed by Spearman's correlation (r = 0.243, P = 0.040). However, we did not observe any statistically significant association between the rs2910164 genotype and either miR-146a or TRAF6 mRNA expression in NSOC cleft lip tissue samples.

miR-146a and TRAF6 are expressed during the development of craniofacial structures
Based on the bioinformatic analysis of RNA-Seq count data in the The two cells were transiently transfected with the constructs (miR-146a G/C allele or control vector and the reporter plasmids containing the TRAF6 3 ′ -UTR fragment) and assayed for luciferase activity after 24 hr. n = 3 for each group, luciferase activity was measured three times for each group. TRAF6 mRNA (c) and protein (d) expression levels in human cell lines after transfection of either the miR-146a G/C allele or control vector. n = 3 (colonies) for each group. The data indicate means ± SD from three independent experiments. *P < 0.05, **P < 0.01, ***P < 0.001. (e) The intensity of the TRAF6 protein bands was analyzed by densitometry after normalization to the corresponding -actin level.
(Supporting Information Figure 2), the expression of miR-146a in E10.5-E14.5 mouse embryos was identified in specific embryonic tissues, including those of facial structures.
To further assess the role of TRAF6 during the development of craniofacial structures, we examined its expression in palate development (E11.5-E15.5) using IHC. TRAF6 was found to be expressed in mouse palatal shelves during the fusion of palatal medial edge epithelium (Supporting Information Figure 3). Our findings were consistent with the results from the online functional annotation enrichment analysis with ToppFun (Chen et al., 2009), in which TRAF6 was identified to be expressed in the mandibular arch and neural epithelium overlying medial eminence in embryonic mouse tissues (E9.5-E10.5) (Supporting Information Figure 4).

DISCUSSION
As key regulators of embryogenesis, miRNAs have been demonstrated to play essential roles in developmental processes including the epithelial-mesenchymal transition, cell migration, differentiation, proliferation, and apoptosis (Eberhart et al., 2008;Pauli, Rinn, & Schier, 2011).
SNPs in miRNA genes are likely to affect the expression of miRNA as well as their binding efficiencies with target genes, thereby contributing to susceptibility to common human diseases (Ryan et al., 2010;Schoen et al., 2017 Here we focused on the five common variants in miRNA genes in a two-stage study and identified a potentially functional variant, rs2910164, located in pre-miR-146a, conferred susceptibility to NSOC. NSOCs are commonly categorized as CL/P and CPO, which have historically been analyzed as distinct disorders due to the different developmental origins of the lip and palate (Jiang, Bush, & Lidral, 2006).
Recently, FOXE1 was shown to be associated with both the CL/P and CPO in a GWAS . In the present study, stratification analysis revealed that the rs2910164 G allele was a risk factor for both CL/P and CPO groups.
In vitro experiments demonstrated that the rs2910164 G allele pro- With the bioinformatic prediction of computation alignment, we found that the TRAF6 3 ′ -UTR was in the potential binding site of miR-146a. The luciferase assay in HEK-293 and HEPM cells revealed that the reporter activity of the G allele was significantly higher compared with that of the C allele. Moreover, cell lines transfected with the rs2910164 G allele expressed lower levels of miR-146a and higher levels of TRAF6. Both miR-146a and TRAF6 were expressed in vivo in NSOC cleft lip tissue samples, and a moderate inverse correlation between them was observed. In sum, these findings indicated that miR-146a might specifically interact with the TRAF6 3 ′ -UTR and significantly inhibited the transcription and translation of TRAF6. However, no significant association between the rs2910164 genotypes and either miR-146a or TRAF6 mRNA expression in NSOC cleft lip tissue samples were observed, possibly due to limited sample size.
TRAF6 is involved in various signaling pathways and cell progressions, and acts as an important adaptor protein (King et al., 2006).
The TRAF6/TAK1/p38 MAPK signaling pathway activated through T RI/T RIII is responsible for cranial neural crest cell proliferation defects and the failure of palatal fusion (Iwata et al., 2012). In this study, TRAF6 was found to be expressed in the mouse fusion palate. Therefore, we speculated that the miR-146a rs2910164 G allele was associated with the expression of miR-146a and could influence cell proliferation and modulate TRAF6 expression, leading to the failed fusion of facial processes and finally contributing to the occurrence of NSOC (Supporting Information Figure 5).
Several limitations of the present study should be addressed. First, larger cohorts from different populations will be required to validate the association between rs2910164 and NSOC. Additionally, due to the selection criteria of the analyzed SNPs, potential interesting variants were missed. Furthermore, because numerous different cell types, such as neural crest cells, epithelial cells, and mesenchymal cells, are involved in craniofacial development Tian et al., 2017), more experiments from other cell types are needed to verify our findings.
In conclusion, we found that the base change in rs2910164 results in a modification of the structure of the stem region of pre-miR-146a, which is associated with the risk of NSOC. Our findings will help in understanding the etiology and genetic origin of NSOC.