{'fraction_non_alphanumeric': 0.04052251820948603, 'fraction_numerical': 0.020038983688404063, 'mean_word_length': 4.754427390791027, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 9, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Preliminary study on the molecular features of mutation in multiple primary oral cancer by whole exome sequencing

Multiple primary cancers (MPCs) refer to cancers that occur simultaneously or metachronously in the same individual. The incidence of MPC has increased recently, as the survival time of malignant tumor patients has been greatly prolonged. It is difficult to differentiate MPC from primary cancers (PCs) in the same anatomical region from the clinical manifestation alone. However, their biological behaviors appear to be distinct. In this study, we show that the prognosis of multiple primary oral cancers (MP-OCs) is worse than primary oral cancers (P-OCs). To better understand the molecular mechanisms of MP-OC, we used whole exome sequencing (WES) to analyze samples from 9 patients with MP-OC and 21 patients with P-OC. We found more somatic mutations in MP-OC than in P-OC. MP-OC had more complicated mutation signatures, which were associated with age-related and Apolipoprotein B mRNA Editing Catalytic Polypeptide-like (APOBEC) activity-related signatures. Tumor mutational burden (TMB) and mutant-allele tumor heterogeneity (MATH) of MP-OC trended higher compared to P-OC. KEGG and GO analysis showed the differential pathways of MP-OC versus P-OC. In addition, MP-OC took amplification, not loss, as the main pattern of copy number variation (CNV), while P-OC took both. Lastly, we did not find significantly different mutant germline genes, but MSH-6 mutation may be a potential MP-OC driver. In short, our preliminary results show that MP-OC and P-OC have different molecular characteristics.


Introduction
Multiple primary cancers (MPCs) are cancers that occur simultaneously or metachronously in the same individual. Metastasis of initial primary cancers should be carefully excluded (1). As early as 1889, Billroth first reported a case of multiple primary gastric cancer after epithelial carcinoma of the outer ear (2). In 1932, Warren and Gates formally determined the diagnostic criteria of multiple primary cancers (3). MPC was once regarded as a rare disease (4). With the improvement of early screening, diagnostic techniques, treatment methods, postoperative follow-up, and supportive treatment, survival following malignant tumors has been greatly prolonged. Meanwhile, the incidence rate of MPC has significantly increased (5)(6)(7). Without a clear medical history, there are few obvious differences between MPC and primary cancers (PCs) in terms of their clinical manifestations. However, Kai-Hsiung et al. reported that the 5-year overall survival and 5year disease-free survival were worse in individuals with multiple primary lung cancers than those with a primary lung cancer (8). Qi-Wen et al. showed that compared with matched cases, esophageal squamous cell carcinoma accompanied with synchronous MPC was related to significantly impaired survival and an elevated risk of locoregional disease progression following concurrent chemoradiotherapy (9). It thus appears that the underlying biological behavior of these two clinical groups might be distinct. A better understanding of the clinical features and molecular mechanisms of multiple primary cancers may be beneficial to cancer management. However, current studies on multiple primary cancers had mostly focused on its clinicopathological features (10). There are only a handful of studies that reveal the differences in the molecular mechanisms of MPC.
Gene mutations, including single-nucleotide polymorphism (SNP), insertion-deletion (InDel), and copy number variation (CNV), play an important role in the occurrence and progression of tumors. For example, the RB1 gene mutation results in retinoblastoma (11). Detection of the mutation helps to reveal the molecular mechanisms of MPCs and to promote prognostication and treatment planning. Whole exome sequencing (WES) is the most common method used to detect gene mutations. Tumor mutational burden (TMB) and microsatellite instability (MSI), which have been approved as pan-cancer biomarkers by the FDA, can also be obtained by WES (12,13). Patients with high TMB (TMB-H) or high MSI (MSI-H) benefit from immunotherapy.
In this study, we found that the prognosis of individuals with multiple primary oral cancers (MP-OCs) was worse than those with primary oral cancers (P-OCs) (8). We used WES to analyze 9 patients with MP-OC and 21 patients with P-OC, to better understand the underlying molecular mechanisms. The mutant germline genes of two groups were also analyzed.

Dataset and bioinformatics analysis
The SEER database is currently the largest publicly available cancer database, covering approximately 28% of the US population (14). We chose SEER Research Data, 9 Registries, Nov 2020Sub (1975 to calculate the proportion of MP-OC in oral cancer (OC). Then, we analyzed the survival of MP-OC and P-OC. SEER*Stat software (version 8.3.92; National Cancer Institute, Bethesda, MD, USA) was used to analyze the proportion of MP-OC from 1975 to 2018.
The range of OCs included tongue, floor of mouth, and gum and other mouth. According to the International Classification of Diseases for Oncology, third edition (ICD-O-3), the codes for tongue cancer are C019-C029, those for floor of mouth cancer are C040-C049, and those for gum and other mouth cancers are C030-C039, C050-C059, and C060-C069.
Data on sex, site of cancer, histopathological type, survival time, and cause of death were collected. OC patients with a history of malignancies were classified as MP-OC, and the remaining cases were classified as P-OC. We excluded patients with ambiguous information on survival, which was labeled with "N/A not seq 0-59".

Patient cohort
A total of 317 patients with OC from the Hospital of Stomatology of Sun Yat-sen University (Guangzhou, China) were enrolled from October 2014 to September 2016 to further confirm the different prognoses of MP-OC and P-OC. Of these, 16 were MP-OC patients. The inclusion criteria of MP-OC followed the SEER criteria (15): multiple neoplasms were malignant with a definite pathological diagnosis; multiple cancers come from different tissue origins, or the edges of the cancer lesions are at least 2 cm apart from each other; if multiple cancers were of the same site and tissue origin, their times of diagnosis had to be at least 5 years apart. Recurrence and metastasis were carefully excluded. The endpoint of follow-up was 31 December 2020.
A total of 64 patients were enrolled for WES and subsequent analyses, including 20 with MP-OC and 44 with P-OC. Fresh cancerous tissues were collected from 30 patients for somatic genome mutation analysis, including 9 with MP-OC and 21 with P-OC. Blood or normal mucosa was taken as controls. We also collected blood or mucosal samples from all remaining participants. These samples were sequenced for germline gene mutations. The samples were stored in liquid nitrogen or at −80°C until WES. The DNA libraries were sequenced on the Illumina sequencing platform by Gene Denovo Honour Biotechnology Co., Ltd (Guangzhou, China).

Nucleic acid preparation
DNA was extracted using the Cetyltrimethyl Ammonium Bromide (CTAB) method to generate libraries for Next-Generation sequencing. NEBNext dsDNA Fragmentase (NEB, Ipswich, MA, USA) was used to fragment genomic DNA ends followed by DNA end repair. Then, NEBNext adaptor (NEB, Ipswich, MA, USA) was used to dA-tail and ligate these DNA fragments. Biotinylated RNA library baits and magnetic beads were mixed with the barcoded library to select targeted regions with the SureSelect Human All Exome V6 Kit (Agilent Technologies, Palo Alto, Calif.). Next, we further amplified the captured sequences for 150-bp paired-end sequencing using the Illumina X-ten system (Illumina, San Diego, CA, USA).

Clean reads filtering
Fastq (16) was used to filter the raw reads. The low-quality parts of raw data were identified for removal as follows: reads with ≥10% unidentified nucleotides; reads with > 50% bases having phred quality scores of ≤20; and reads aligned to the barcode adapter. The quality of clean data after filtering is shown in Supplementary Figures S1A, B.

Somatic SNV and InDel identification
The Burrows-Wheeler Aligner (BWA) (17) completed the alignment of the clean reads against the human reference genome (GRCh38). Somatic SNV and InDel of multisample were called by MuTect (18).

Analysis of mutational signatures
The non-negative matrix factorization (NMF) method was used to obtain mutational signatures (19). There were 96 possible mutation types combined to obtain a matrix representing each feature. By cosine similarity, the identified signatures were compared to the previously determined consensus signatures (20) and confirmed.

Identification of mutated genes
Considering all variations, including somatic SNV and InDel, we identified the set of mutated genes. Genes with a p-value<0.05 were considered to be with high-frequency mutations. Significantly mutated genes (SMGs) refer to the genes whose mutation frequency was significantly higher in cancer than in controls. We used MuSiC (21) to search for SMGs. MuSiC is a method based on the mutational recurrence on all cancer samples. A convolution test was carried out on each mutation type. Only the mutated genes with FDR<0.01 were considered to be SMGs.

Gene functional enrichment analysis
The set of mutated genes with a p value<0.05 were used for gene functional enrichment analysis. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (22) and Gene Ontology (GO) (23) are commonly used databases. We used an in-house script that incorporated both databases for functional annotation analysis. The significance of gene group enrichment was defined by a modified Fisher's exact test and Q value< 0.05 was considered to indicate a statistically significant difference.

Calculation of TMB, MSI, and MATH
Nonsynonymous mutations in the coding region were selected to assess TMB. We used 38 Mb as the estimate of the exome size. All nonsynonymous mutations/exome lengths indicate the TMB of each sample (24).
Microsatellites are repetitive sequences with unit lengths varying from 1 to 6 bp. We evaluated the MSI following the MSIsensor method described by Niu and Ding (25).
Mutant-allele tumor heterogeneity (MATH) is a quantitative assessment of genetic intra-tumor heterogeneity, based on WES of tumors and matched normal tissues. The MATH score was calculated following the method described by Mroz and Rocco (26). Briefly, the MATH value of each tumor was calculated from the median absolute deviation (MAD) and the median of its mutant-allele fractions at tumor-specific mutated loci using: MATH = 100 * MAD=median:

Somatic CNV detection and analysis
VarScan 2 (27) was used to identify somatic CNVs with the following parameters: minimum coverage ≥20 and phred base quality ≥20. DNA copy packages in R/Bioconductor were used to normalize the log ratio of coverage between tumor and normal samples by circular binary segmentation. Recurrent somatic CNVs were identified by using GISTIC (28). GISTIC broad-level analysis was performed with a size threshold of 98% of a chromosome arm to differentiate between focal events and arm level. q ≤ 0.25 was considered significant for arm-level events and CNV regions. Then, residual q-values q ≤ 0.05 were used to determine the significance of focal CNV events. All oncogenes and tumor suppressor genes were identified according to the OncoKB ™ (https://www.oncokb. org/cancerGenes).

Detection and analysis of germline mutant genes
According to the study of Kuan-Lin Huang et al. (29), 152 genes were identified as cancer predisposition genes (Supplementary Table S1). We annotated and prioritized pathogenic variants in these 152 genes in our 64 patients.
We set the filter conditions as follows: firstly, only variants with a structure type of exonic or splitting were retained; we made requirements for allelic depth (AD) and read depth (RD), requiring AD ≥5 and AD/(RD + AD) ≥0.2; we then removed variants potentially contaminated by tumor-adjacent tissues. Variants were removed if they were found in both cancer and tumor-adjacent tissues, and the allele fraction in tumor-adjacent tissues was<0.3; finally, we collected all SNPs called by GATK and VarScan2, and indels that were called by at least two of the three callers (GATK, VarScan2, and Pindel).
We used several methods to score these variants, and finally obtained the sum of the scores of the variants. A higher score indicated higher pathogenicity. A final score of ≥1 was classified as pathogenic, while variants with scores of< 1 were discarded. The specific methods are as follows: Firstly, we required variants with a dominant mode of inheritance; variants would be removed if they did not meet this requirement.
Secondly, the pathogenicity of polypeptide changes was annotated in ClinVar; variants were given +7 points if consistent with previously identified pathogenic variants in ClinVar. Variants that resulted in a different amino-acid change at the same position were given +2 points.
Thirdly, we used HotSpot3D (35) to analyze the protein structure and pathogenicity; if a germline variant was found to be a recurrent somatic mutation among all cancer types in a HotSpot3D cluster, the variant was given +2 points.
Fourthly, we scored according to mutation frequency. Variants that did not exist or showed a very low frequency (MAF< 0.0005) in the ExAC dataset were given +2 points; common variants (MAF > 0.05) were given −8 points.
Finally, we used SIFT and PolyPhen for scoring. If a variant was identified as "damaging" or "deleterious" in SIFT (score<0.05) and "probably damaging" in PolyPhen (score >0.432), it was given +1 point. Conversely, if SIFT is greater than 0.05 and PolyPhen is less than 0.432, the 1 point would be taken away.

Statistical methods
Statistical analysis was carried out by using IBM SPSS Statistics 25.0 (IBM Corp., Armonk, NY, USA). Clinical data including age, gender, site of OC, smoking history, and excessive drinking history were compared between P-OC and MP-OC using the Mann-Whitney U test, Chi-square test, or Fisher's exact test. The difference of sequencing results between two groups such as TMB, MATH, and MSI was also analyzed by the Mann-Whitney U test, Chi-square test, or Fisher's exact test. The survival rate of MP-OC and P-OC was calculated by using the Kaplan-Meier method, and the survival curve was generated. We further performed univariate regression analyses and a multivariate Cox regression analysis. Candidates were subjected to a multivariate Cox proportional hazards regression with the enter selection of variables (a = 0.1), and the forest plot was drawn to show the hazard ratio of risk factors. The log-rank test was used to formally test the differences. Significance was determined based on a p< 0.05 threshold.

The proportion of MP-OC among oral cancers has increased over time
Data from SEER showed that the proportion of MP-OC among OCs increased over time. Among all OC patients between 1975 and 1985, MP-OC patients accounted for 12.73%, while the proportion increased to 25.32% between 2008 and 2018 ( Figure 1A).

MP-OC shows worse prognosis than P-OC
We first investigated the patients in the SEER database. The survival of patients with MP-OC was significantly worse, with a median survival time of 70.52 months (95% CI, 68.52-72.52), compared with 119.02 months (95% CI, 117.50-120.53) in the P-OC group (Kaplan-Meier survival analysis, p< 0.001) ( Figure 1B). To confirm this result, we enrolled 317 patients with OC from our center; their clinicopathologic features are listed in Table 1 (further details are provided in Supplementary  Table S2). There was no significant difference in the clinicopathologic features between the two groups (p > 0.05, Chi-square test). The median follow-up time of the 317 patients was 59 months. Our cohort paralleled our earlier findings; we found that the survival times were significantly shorter in MP-OC than in P-OC (p< 0.05, Kaplan-Meier survival analysis) ( Figure 1C). The median survival time of patients in MP-OC was 39.79 months (95% CI, 26.27-53.30) and 66.95 months in P-OC (95% CI, 63.33-70.57). To determine the independent risk factors for Kaplan-Meier results, univariate regression analyses and multivariate Cox regression analysis were used. The results showed that history of cancer and N stage were independent risk factor factors of OC patients ( Figure 1D).

General mutation characteristics of MP-OC and P-OC
To explore the mechanisms underlying the difference in prognosis, 9 MP-OC and 21 P-OC tumor samples were analyzed by WES. The clinicopathologic features of the cases are listed in Table 2 (further details are provided in Supplementary  Table S2).
An average of 91.17 million reads was acquired with an average base quality of 97.45% (Q20). The ratio of high-quality clean reads ranged from 99.24% to 99.84% (Supplementary Figure S1A). The average ratio of on-target coverage was 84.45%. The average target depth ranged from 66.16× to 452.29×, and 90% samples reached 80×. The quality control data for each case are provided in Supplementary Table S3.
Among the 30 collected patient tumors, we observed 6,761 non-synonymous SNPs and 2,272 synonymous SNPs (Supplementary Table S4). The type and distribution of SNPs and InDels are listed in Supplementary Table S4 Table S4). The most common base substitution in both groups was C > T. More C>G mutations were observed in MP-OC than in P-OC (p<0.05) (Figure 2A, middle; Supplementary Table S4).
CNV gain numbers of MP-OC trended higher than in P-OC, although the difference was not significant, and CNV loss was significantly less than in P-OC (p< 0.05, Mann-Whitney test).  Table S4).
Lastly, we conducted pathways analysis using GO and KEGG with the mutated genes in MP-OC and P-OC. The two groups shared some of the pathways. In the cellular component (CC), molecular function (MF), and biological process (BP) parts of GO pathways, there were zero, one, and four pathways unique to the MP-OC group ( Figure 4A). Among KEGG pathways, there were two pathways unique to the MP-OC group, including drug metabolism-cytochrome P450 and metabolism of xenobiotics by cytochrome P450 (Figure 4B; Supplementary Table S6).

Different CNV patterns in MP-OC and P-OC groups
We evaluated DNA fragment changes by CNV. We observed that mainly CNV gains rather than CNV losses occurred in MP-OC. Both CNV gain and CNV loss, and mixed patterns were observed in the P-OC group ( Figure 5A; Supplementary Table  S4). We highlighted some of the chromosome fractions that were different between the two groups. In these fragments, there were 99 oncogenes and 89 tumor suppressor genes, while 14 genes act as a double-edged sword (Supplementary Table S7). We labeled some genes that may impact the pathogenesis of MP-OC ( Figure 5B).

Germline mutations in MP-OC and P-OC group.
We also tried to analyze whether germline mutations were the possible cause of MP-OC. The clinicopathological features of 64 patients are shown in Supplementary Table S2. In our preliminary study, we found no significant differences between General mutation characteristics of 9 MP-OC and 21 P-OC patients. All somatic substitutions identified in this study were included, and the bar graph showed the subtype of SNP and InDel (left), proportion of base substitution types (middle), and number of copy number variation fragments (right) between MP-OC and P-OC groups, respectively. * indicates significant difference between the two groups (Mann-Whitney test, p< 0.05) (A). TMB and MSI of 9 MP-OC patients and 21 P-OC patients. The bar chart (left) shows the details of 30 patients, and the violin diagram (right) shows the comparison of MSI and TMB between the two groups, respectively. Mann-Whitney test, ns stands for no significance, which indicates p > 0.05. The cutoff between TMB-H and TMB-L is 20%, and the cutoff between MSI and MSS is 3.5% (B, C). the two groups (p > 0.05, nonparametric Mann-Whitney test) ( Table 3; Supplementary Table S8).

Discussion
Firstly, we analyzed data from the SEER database. We found that in recent years, the proportion of multiple primary cancer in the total OC pool has increased. The oral cavity is one body region with a high incidence of multiple primary cancers (36). Without a medical history, it is difficult to distinguish multiple primary cancers from primary cancers in the oral cavity. These distinctions are critical, as our study shows that the prognosis of MP-OC patients is worse than P-OC patients; this finding is consistent with previous studies (8, 37,38). However, the molecular mechanisms remain unexplored. With an increasing incidence of multiple primary OCs, there is growing importance in understanding the molecular mechanisms of poor prognosis. Therefore, we attempted to characterize the differences in the molecular mechanisms of MP-OC and P-OC by utilizing WES.
In general, mutations were more frequent in the MP-OC group. This group had a trend towards greater SNP, TMB, and MATH. TMB refers to the number of mutations per million bases in cancer. The MATH value was calculated from the MAD B A FIGURE 4 Pathway analysis of MP-OC and P-OC patients. Histogram showing the difference of GO between the MP-OC (upper) and P-OC (lower) group. Venn diagram showing the differences of molecular function, cellular component, and biological process in GO, respectively (A). Lists of the significant KEGG pathways in the two groups. Wayne diagram showing the differences between the two groups (B). B A FIGURE 5 CNV in MP-OC and P-OC groups. Histogram showing CNV gains (red) and CNV loss (blue) in 9 MP-OC patients and 21 P-OC patients. Patients with CNV gains and CNV loss in the two groups were counted in the pie chart (A). Composite of copy number profiles for the MP and P group, with gains in red and losses in blue. The regions that showed a difference in the frequency of copy number alterations between two groups were shaded in light gray rectangles; some genes in these regions are labeled (B). and the median of its mutant-allele fractions at tumor-specific mutated loci (39). MATH evaluates the genetic differences between tumor cells. Previous studies showed that TMB and MATH are not significantly different in multiple gastric cancer and multiple primary lung cancer (39,40), and our research supports these findings. However, the situation may be more complicated in MP-OC than in P-OC. A larger sample size is essential to confirm and refine these conclusions. In addition, we also compared the MSI of the two groups. No significant differences or tendencies were found. MSI refers to short and tandem repeat base sequences that have errors in the replication process, usually arising from DNA mismatch repair defects (dMMR). Previous studies found that colorectal cancer and head and neck cancer (HNC) with MSI-H show a higher risk of MPC (41,42). Multiple primary gastrointestinal cancer exhibits frequent MSI (43,44). MSI-H is not common in OC. A study across 39 cancer types showed an incidence of only 0.78% in head and neck squamous cell carcinoma (45). Our research shows similar results in OC.
Mutational signatures are different combinations of mutation types, often originating from different mutational processes (46). The mutation types of MP-OC and P-OC groups are shown in Figure S2. In our study, the mutational signatures of the MP-OC group were more abundant than in the P-OC group, which might result from more complicated pathogenic factors. Signature 1 was closely related to age. The age at diagnosis of patients with MPCs is older than that of patients with the same primary cancer (47, 48). Accumulated mutations accompanied by aging might play an important role in the occurrence of MPCs. Signatures 2 and 13 were related to the activity of the AID/APOBEC cytidine deaminase family. The APOBEC family has been shown to induce tumor mutations by aberrant DNA editing mechanisms (49). Infection with high-risk HPV, which is one of the pathogenic factors in the oral and maxillofacial regions, leads to increased APOBEC family activity (50)(51)(52).
We also found that the mutation spectrum between the two groups was significantly different. Among the top 30 frequently mutated genes, only 5 overlapped between groups. Similar results were found in our GO and KEGG pathways analysis. We carefully reviewed the specific KEGG pathways peculiar to MP-OC, in which the cytochrome P450 pathway is involved. Various studies have shown the significance of CYP 450 polymorphism in the susceptibility to cancer. In hormone-related cancers such as breast cancer (53) and prostate cancer (54), CYP 450, which is involved in steroid hormone metabolism, affects cancer susceptibility. Studies have shown that patients carrying a mutation of the  (58). OSCC patients who use tobacco also have some relevant characteristics, such as higher C>G transversion mutations (59). In our study, no smoker or betel nut chewer was enrolled into the MP-OC group for sequencing, which might have biased our results. Samples involving MPC smokers or nut-chewing patients will be collected in our future studies. CNV refers to a form of genomic structural change that results in abnormal gene copy numbers (46). It is an important way to regulate the expression and function of oncogenes and tumor suppressor genes, which lead to the evolution and progression of cancer (46, 60-63). We found different CNV patterns between the two groups. The main CNV pattern of MP-OC patients was gain. Cancer patients with CNV gain have a poorer prognosis (64), which might be one of the reasons why MP-OC patients have a poorer prognosis.
The causes of MPCs include accumulation of mutations accompanied by aging, field cancerization caused by physical and chemical factors, virus infection including HPV, iatrogenic effects of radiotherapy and chemotherapy, and germline mutations (41,(65)(66)(67). We preliminarily investigated germline mutations as a potential cause. Germline mutation refers to a congenital gene mutation in germ cells and leads to a higher risk of canceration of somatic cells. However, we found no significant differences in germline mutations between MP-OC and P-OC.

Conclusions
The biological behavior and molecular mechanism underlying MPC were different from primary cancers. MP-OC patients had a worse prognosis. By WES and comparative analysis of 9 MP-OC patients and 21 P-OC patients, our study described the differences in somatic point and fragment mutations between the two groups. TMB and MATH in the MP-OC group trended higher than in P-OC. The MP-OC group had more complex mutation signatures, suggesting that agerelated factors and the AID/APOBEC pathway play different roles in MP-OC versus P-OC. The mutated genes of MP-OC were enriched in the cytochrome P450 pathway of KEGG. Amplification is the main mode of CNV in MP-OC. Our study did not find germline mutations that play a key role in the pathogenesis of MP-OC. The mechanisms underlying MPC remain appealing but unresolved. There are several limitations in our preliminary study. Smoking and nut-chewing patients were not involved in the sequencing study. A larger cohort that includes smoking and nut-chewing patients will be essential for further studies. Samples should be carefully stratified and analyzed according to the different potential causes, including aging, field cancerization, radiation, and germline mutations. Moreover, WES alone is unable to resolve molecular mechanisms. Multi-omics studies incorporating analyses of the transcriptome and proteome will need to be incorporated to tease out the biological effects of specific mutations. In vivo and in vitro experiments need to be performed to confirm the specific potential driving genes and pathogenic pathways.

Data availability statement
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA-Human: HRA002206) and are publicly accessible at https://ngdc.cncb.ac.cn/gsa-human/s/ov8sC0i6.

Ethics statement
The studies involving human participants were reviewed and approved by Ethics Committee of Hospital of Stomatology, Sun Yat-sen University. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any data included in this article.

Author contributions
KL and JG drafted the article and reviewed the submitted version of the manuscript. QZ worked on the acquisition of data. KL, JG, QZ, LY, and XM worked on the analysis and interpretation of data. LY, XM, JC, GL, and YL critically revised the article and worked on the statistical analysis. GL and YL are responsible for administrative/technical/material support and study supervision. KL, JG, and YL worked on conception and design. YL and GL approved the final version of the manuscript on behalf of all authors. All authors contributed to the article and approved the submitted version.