In this study, the genotyping data of single-nucleotide polymorphisms (SNPs) at SLC23A2 in 1,047 patients with NSCLO and 2,255 normal controls were used from two previous GWASs (17, 18). The human subject study protocols were reviewed and approved by the institutional review board (IRB) of West China Hospital of Stomatology, Sichuan University, in 2016 (WCHSIRB-D-2016-012R1) and conformed to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) guidelines. Written informed consent was obtained from recruited individuals of consenting age and from parents on behalf of their participating children.

Considering the important role of oral epithelium in facial morphogenesis and its known association with NSCL/P, a human oral epithelial-derived cell line (GMSM-K, kindly gifted by Dr. Zhang from Peking University) was selected for functional analysis in our study (19). GMSM-K was cultured in Dulbecco's modified Eagle's medium (DMEM) supplemented with 10% fetal bovine serum (FBS; Gibco, United States) and 1% penicillin-streptomycin solution (Gibco, United States). Small interfering RNA (siRNA) targeting SLC23A2 (NM_005116.6) and negative control siRNA were both designed and synthesized by GenePharma (Shanghai, China). Information related to the siRNA can be found in the Supplementary Materials. Following the manufacturer's instructions, GMSM-K cells were seeded in a 6-well culture dish at a density of 1 × 10⁵ per well. When the cells reached a confluence of 70%–90%, siRNA was transfected into the GMSM-K cells using Lipofectamine 3000 (Thermo Scientific, United States). After 6 h, we replaced the Lipofectamine 3000-containing media with fresh complete media for further culturing. Following 48 or 72 h of transfection, the cells were collected to perform further examinations. The effective duration of siRNA-mediated gene silencing is 5–7 days.

GMSM-K cells were seeded on a 6-well plate, rinsed with phosphate-buffered saline (PBS), and fixed with ice-cold methyl alcohol for 5 min. Next, the cells were permeabilized with 0.25% Triton X-100 for 5 min, washed with PBS twice, and blocked with 2.5% bovine serum albumin in PBS for 1 h. Antibodies against SLC23A2 (Novus, NBP2-13319) were diluted 150-fold with PBS and incubated at 4°C overnight.

Cytosolic ROS were detected by staining the GMSM-K cells with 10 μM 2’,7’-dichlorodihydrofluorescein diacetate (DCFH-DA) (Sigama,Germany) in serum-free medium for 30 min at 37℃. The cells were washed twice with PBS. Under a fluorescence microscope, a fluorescein isothiocyanate (FITC) filter was used to observe fluorescence. Superoxide dismutase (SOD) activity and glutathione/glutathione disulfide (GSH/GSSG) were detected via a specific kit (Beyotime Biotechnology, China) following the manufacturer's instructions.

GMSM-K cells were seeded into 96-well plates at a density of 2 × 10⁴ cells/100 μL. At 21, 45, and 69 h after transfection, a mixture of 10 μL of Cell Counting Kit-8 (CCK-8) (APExBIO, United States) and 90μL DMEM was added to each well after removing the original medium and the cells were further incubated at 37℃ for 3 h. The optical density (OD) of the mixture was measured at a wavelength of 450 nm. The experiments were repeated three times and five parallel holes were set in each experiment.

GMSM-K cells were seeded into 6-well plates at a density of 2 × 10⁵ cells/2 mL. After transfection, five horizontal lines were scored on the bottom of the plates and the cells were scratched perpendicularly to the horizontal line using a sterile 20-µL pipette. The cells were then washed three times with PBS, placed under a microscope, and the medium was changed to 0.1% FBS experimental medium. The scratches were observed after incubation at 37°C for 0, 24, and 48 h.

The apoptosis rate was evaluated using the Annexin V-PE/7-amino-actinomycin D (7-AAD) Apoptosis Detection Kit (Vazyme, China) according to the instructions from the manufacturer. The cells were seeded into 6-well tissue culture plates (2 × 10⁵ cells/well). Following treatment, the cells were collected, washed with PBS, and resuspended in 500 μL binding buffer. Then, 5 μL Annexin V-PE and 5 μL Annexin V-PE were added to the buffer and incubated at room temperature for 10 min in the dark. Cells were analyzed using flow cytometry (Thermo Fisher, United States) within 1 h. Flow Cytometry Standard (FCS) files were downloaded and analyzed using FlowJo software (version 10.4).

GMSM-K cells were transfected with a siRNA-negative control or siRNAs-SLC23A2 for 48 h. Then, RNA was extracted from the cells and RNA sequencing (RNA-seq) was performed using the BGISEQ-500 platform (BGI, China). Three biological replicates were included within each group. Differential gene expression analysis was performed using the DESeq2 method (|log2| ≥ 0.8, q-value ≤ 0.05), and Gene Ontology (GO)/Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed using ChiPlot (https://www.chiplot.online/).

RNA was extracted 48 h after transfection using RNA-easyTM Isolation Reagent (Vazyme, China), and was then reverse-transcribed to cDNA using a PrimeScript™ RT reagent Kit (Takara Biotechnology, China). Real-time quantitative PCR (RT-qPCR) was performed using TB Green® Premix Ex Taq™ (Takara Biotechnology, China) on a LightCycler 480 System (Roche, Switzerland). All the experiments were performed in triplicate, each with three technical replicates. The results were calculated using the 2^−ΔΔCt equation, normalizing values to GAPDH within each sample. The primers used are shown in Supplementary Table 1.

A gene-environment interaction model in the GMSM-K cell line was established by knocking down SLC23A2 and adding Sin-1 (20) and AA simultaneously. There were the following six groups: Negative control (NC), NC + Sin-1 (Aladdin, China), NC + Sin-1 + AA (Sigama, Germany), siRNA (si), si + Sin-1, and si + Sin-1 + AA. The detection methods of oxidative stress and apoptosis are the same as above.

All the animal experiments performed were approved by the Animal Ethical and Welfare Committee of Nanjing University with ID IACUC-D2310004 (2023.10.8). First, we used the previously verified morpholino (MO) technology to construct the slc23a2-knockdown zebrafish model. MO targeting at slc23a2 (slc23a2-MO) (5'-GCACTGAATATGAAAAGATTGTACT-3’) was designed and produced by Gene Tools (United States). According to the preliminary experiment, the final concentration of slc23a2-MO was 2 ng/μL. Injections were carried out at the single cell stage, and after 8 h, unfertilized eggs and dead eggs were removed and replaced with fresh medium. At 48 h postfertilization (pf), the zebrafish embryos were collected, incubated with different concentrations of 3-Morpholinosydnonimine(Sin-1)(exogenous ROS mimics) and AA, treated with a 20.5 μM DCFH-DA probe (a chemically reduced form of fluorescein used as an indicator for ROS) and 5 μg/mL acridine orange (AO; an indicator for apoptosis) and then incubated in the dark at 28.5°C for 1 h. The embryos were then drenched with water three times, anaesthetized with 0.02% tricaine, and photographed under a fluorescence microscope with FITC filters.

The chi-square test and 95% confidence interval (95% CI) for the odds ratios were used to compare the allele frequency between the cases and controls. Each SNP was assessed using the Hardy–Weinberg equilibrium (HWE) and the minor allele frequency (MAF) was calculated. Moreover, the difference in allelic and genotypic frequencies of each SNP between the cases and normal controls was calculated using PLINK software (21). Pairwise linkage disequilibrium (LD), which shows both D′ and R², was computed for all the SNPs using the Haploview program (http://www.broad.mit.edu/haploview/haploview). The results are shown as mean ± SD. A statistical analysis of the in vivo and in vitro experiments was performed using an unpaired two-tailed t-test in GraphPad Prism 8 software.

A total of 306 common SNPs (MAF ≥0.01, with call rates >95%) that passed the HWE threshold (p > 0.05) were used in the association analysis (Supplementary Table 2). Both the allelic (Supplementary Table 3) and genotypic (Supplementary Table 4) association analyses indicated that 10 SNPs located at SLC23A2 were significantly correlated with NSCLO, and were both adjusted for multiple corrections (p = 0.05/306). The pairwise LD results showed that rs6053029 was tightly linked to other SNPs and had the lowest p-value (9.44E-17) in NSCLO (Figure 2). Our previous RNA-seq results revealed that the expression of SLC23A2 in the lip tissues was higher than that in the palate tissues (22) (Supplementary Figure 1).

We first detected the expression of SLC23A2 in GMSM-K cells and found that SLC23A2 was expressed in the cell membrane (Figure 3A). Then, siRNA was used to construct an SLC23A2-knockdown GMSM-K cell model. The siRNA sequence information is as follows: F:GAGCCAUCCUGUCUUUAGATT, R:UCUAAAGACAGGAUGGCUCTT. The qPCR results showed that si-SLC23A2 effectively reduced the transcription level of the SLC23A2 gene (Figures 3B,C). Studies have shown that knockout of Slc23a2 can increase the level of oxidative stress in the embryonic tissues of mice, so we first examined the effect of SLC23A2 knockdown on ROS levels. The data showed that the intracellular ROS levels in the GMSM-K cells significantly increased after SLC23A2 gene knockdown (Figures 4A,B).

The CCK-8 assay showed that there was no significant change in the proliferation level of the GMSM-K cells after knocking down the SLC23A2 gene (Figures 4C,D). The flow cytometry results showed that the early apoptosis levels of the GMSM-K cells increased significantly when the SLC23A2 gene was knocked down (Figures 4E,F). However, cell migration was not significantly affected (Figures 4G,H).

To investigate the potential role of SLC23A2 in the etiology of NSCLO, we performed RNA sequencing of GMSM-K cells with or without SLC23A2 knockdown. Three biological replicates were set up in each group. The differential gene expression analysis identified 342 differentially expressed genes (DEGs) in total, including 64 upregulated genes and 278 downregulated genes (Figure 4I). KEGG analysis showed that a series of biological processes were enriched, including the PI3K-Akt signaling pathway, regulation of the cell cycle, cell senescence, and the Wnt signaling pathway (Figure 4J).

In order to further clarify the relationship between the DEGs enriched in GMSM-K cells and NSOC, we extracted the genotype data of differential genes from two previously published GWASs and conducted an association analysis with each NSOC subtype. The results showed that IGFBP2 (rs9341191), ITGB4 (rs820392, rs820390, rs820389, rs820387, rs866581, rs820388, and rs1008177), LAMC3 (rs3780275), LFNG (rs10261289 and rs375386359), NRARP (rs34679617), COL9A3 (rs2294995), 18 SNPs at TLE2, and 23 SNPs at FLT4 were statistically significant (Supplementary Table 5). The genes that were statistically significant in the association analysis were verified by real-time fluorescent quantitative PCR, and the results showed that the expression of COL9A3 in the PI3K-Akt signaling pathway was statistically different between the control group and the knockdown group (Figure 4L).

This is the first time Sin-1 has been used in GMSM-K cells as an ROS mimicry drug. Therefore, we set the gradient according to the concentration of other cells in other studies, and determined the concentration by taking ROS detection and cell survival rate into account. The results from the DCFH-DA fluorescent probe showed that the effect of Sin-1 on the ROS levels in GMSM-K cells was dose-dependent. When the Sin-1 concentration was 400 , 600 , 800 , or 1,000 μM, there was no significant change in cell survival rate. Therefore, the intermediate concentration of Sin-1 of 800 μM was selected for the following experiments (Supplementary Figure 2A).

The concentration of AA changes its antioxidant effect. Therefore, we added different concentrations of AA to cells stimulated with 800 μM of Sin-1 and observed the altered ROS levels. The data showed that the antioxidant effect of AA was dose-dependent before 250 μM, with no significant increase in antioxidant activity over 250 μM. Thus, an AA concentration of 250 μMwas selected (Supplementary Figure 2B).

In order to explore the effect of AA antagonism on ROS levels in cells with normal and impaired SLC23A2 gene function, we set up six groups for verification. The results showed that the normal SLC23A2 gene function group (NC) had significantly increased ROS levels after Sin-1 stimulation (NC + 800 μM Sin-1). After incubation with AA (NC + 800 μM Sin-1 + 250 μM AA), the ROS level was lower than that after stimulation with Sin-1 alone, and there was no significant difference between the NC and NC + 800 μM Sin-1 + 250 μM AA groups. In the SLC23A2 gene impaired group (si), the ROS level increased significantly after Sin-1 stimulation (si + 800 μM Sin-1), and the ROS level decreased slightly after coincubation with AA (si + 800 μMSin-1 + 250 μM AA) compared with that after pure Sin-1 stimulation (Figures 5A,B). The GSH/GSSG results showed that GSH/GSSG decreased significantly after Sin-1 stimulation in both the NC and si groups, while GSH/GSSG increased significantly after AA supplementation in the NC + 800 μM Sin-1 group and GSH/GSSG did not change significantly after AA supplementation in the si + 800 μM Sin-1 group (Figure 5C). The SOD activity detection results showed that the SOD activity of the NC and si groups significantly decreased after Sin-1 stimulation, while the SOD activity of the NC + 800 μM Sin-1 group significantly increased after AA supplementation and the SOD activity of the si + 800 μM Sin-1 group was not significantly changed after AA supplementation (Figure 5D).

After Sin-1 stimulation, apoptosis increased in both the NC and si groups. Apoptosis was significantly decreased after AA supplementation in the NC + 800 μM Sin-1 group, while there was no significant change in the si + 800 μM Sin-1 group (Figure 5E).

To validate our findings in vivo, we detected the ROS and apoptosis levels in slc23a2-MO zebrafish via DCFH-DA and an AO fluorescent probe. Compared with the wild type (WT) group, the craniofacial ROS and apoptosis levels of the slc23a2-MO group were significantly increased. We consistently observed that 227 μM AA supplementation significantly decreased the ROS level and inhibited apoptosis induced by 400 μM Sin-1 stimulation in the WT group; however, there was no significant change in ROS and apoptosis levels after 227 μM AA supplementation in the slc23a2-MO group (Figures 6A,B).