key: cord-0014056-s87h4yva authors: Huang, Ling-Yu; Hsieh, Yi-Ping; Wang, Yen-Yun; Hwang, Daw-Yang; Jiang, Shih Sheng; Huang, Wen-Tsung; Chiang, Wei-Fan; Liu, Ko-Jiunn; Huang, Tze-Ta title: Single-Cell Analysis of Different Stages of Oral Cancer Carcinogenesis in a Mouse Model date: 2020-10-31 journal: Int J Mol Sci DOI: 10.3390/ijms21218171 sha: 79d29a2c443d0e3ec40a814eedcbffc635622936 doc_id: 14056 cord_uid: s87h4yva Oral carcinogenesis involves the progression of the normal mucosa into potentially malignant disorders and finally into cancer. Tumors are heterogeneous, with different clusters of cells expressing different genes and exhibiting different behaviors. 4-nitroquinoline 1-oxide (4-NQO) and arecoline were used to induce oral cancer in mice, and the main factors for gene expression influencing carcinogenesis were identified through single-cell RNA sequencing analysis. Male C57BL/6J mice were divided into two groups: a control group (receiving normal drinking water) and treatment group (receiving drinking water containing 4-NQO (200 mg/L) and arecoline (500 mg/L)) to induce the malignant development of oral cancer. Mice were sacrificed at 8, 16, 20, and 29 weeks. Except for mice sacrificed at 8 weeks, all mice were treated for 16 weeks and then either sacrificed or given normal drinking water for the remaining weeks. Tongue lesions were excised, and all cells obtained from mice in the 29- and 16-week treatment groups were clustered into 17 groups by using the Louvain algorithm. Cells in subtypes 7 (stem cells) and 9 (keratinocytes) were analyzed through gene set enrichment analysis. Results indicated that their genes were associated with the MYC_targets_v1 pathway, and this finding was confirmed by the presence of cisplatin-resistant nasopharyngeal carcinoma cell lines. These cell subtype biomarkers can be applied for the detection of patients with precancerous lesions, the identification of high-risk populations, and as a treatment target. as differences between isocitrate dehydrogenase mutants astrocytoma and oligodendroglioma in distinct tumor microenvironments and their signature genetic events. Both tumor types share similar developmental hierarchies and lineages of glial differentiation [12] . Moreover, for head and neck squamous cell carcinoma, subtypes are characterized by their malignant and stromal composition, and p-EMT has been established as an independent predictor of nodal metastasis, grade, and adverse pathologic features [13] . Oral carcinogenesis is the sequential development of the normal oral mucosa into potentially malignant disorders and finally into oral cancer. In clinical practice, the carcinogenesis of cell subtypes cannot be analyzed in the same patient because of ethical concerns that we cannot let the patient's precancerous lesion develop into cancer without any intervention. The variations that exist among individual patients also make it difficult to compare the carcinogenesis sequence of cell subtypes between individuals. These concerns limit the analytic power of scRNA-seq. A mouse model with arecoline and 4-nitroquinoline 1-oxide (4-NQO) cotreatment, which mimics the etiology of oral cancer in humans, was utilized to better investigate oral cancer carcinogenesis [14] . Mice were exposed to drinking water with 4-NQO (200 µg/mL) and arecoline (500 µg/mL) and were sacrificed at 8, 16, 20, or 29 weeks. Except for the mice treated for 8 weeks and sacrificed at 8 weeks, all were treated for 16 weeks and then given water for the remaining weeks until sacrificed. Lesions on the tongue were excised for analysis. The animal model applied in this research minimized individual heterogeneity and thus, it can be used during the analysis of different carcinogenesis cell subtypes in different patients. Suitable biomarkers must be identified for detecting oral mucosa cells that have gradually formed precancerous lesions and preventing them from evolving into malignant cancer cells. In clinical treatment, patients are targeted on the basis of the average expression of all genes in a whole tissue. However, an increasing number of studies have confirmed that the entire tumor ecosystem is highly complex and consists of a variety of cell groups with varying abilities to worsen disease progression. Through scRNA-seq gene detection and analysis, a detailed and holistic understanding of the entire cell population and different cell subpopulations can be obtained, and the primary factor influencing the survival of patients with tumor progression and metastasis can be effectively identified. In this study, 4-NQO (200 µg/mL) and arecoline (500 µg/mL) were given to mice in their drinking water every day to simulate the carcinogenic effects of cigarettes and betel nut in transforming healthy oral mucosa cells into precancerous lesions and eventually malignant lesions ( Figure 1A ). Figure 1B displays changes in the tongues of mice. In terms of their gross appearance, the tongues exhibited white lesions after 16 weeks of treatment. Hematoxylin and eosin (H&E) staining results revealed squamous hyperplasia of tongue cells after 16 weeks of treatment ( Figure 1C ), tongue tumors ( Figure 1B ) and tongue squamous cell carcinoma ( Figure 1C ) after 16 weeks of treatment, and 4 weeks of follow-up and invasive tumor formation at the end of the experimental period ( Figure 1B ). H&E staining results at 29 weeks revealed advanced tongue squamous cell carcinoma ( Figure 1C ). Tongue lesions at 16 and 29 weeks were excised for scRNA-seq, RNA, and protein expression analyses. The results indicated the poor differentiation of squamous cell carcinoma at 16 weeks and advanced squamous cell carcinoma at 29 weeks. Figure 1 . The carcinogen 4-nitroquinoline 1-oxide (4-NQO) (similar to carcinogens found in cigarettes) and arecoline (the carcinogen found in betel nut) were used to induce head and neck cell carcinogenesis in mice. (A) Four experimental and control groups were examined. Mice in the experimental group were fed drinking water containing 4-NQO and arecoline and then monitored until being sacrificed. (B) The progression of precancerous lesions on the tongue to cancer was monitored when mice were alive, and they were sacrificed at different time points. With the increasing duration of treatment in the experimental group, white lesions began to appear on the tongues of mice. Tumors were observed on the tongues of mice at week 16 and had become invasive carcinoma by week 29. (C) The tongues of mice in the control and experimental groups underwent hematoxylin and eosin staining for histopathological examination. The carcinogen 4-nitroquinoline 1-oxide (4-NQO) (similar to carcinogens found in cigarettes) and arecoline (the carcinogen found in betel nut) were used to induce head and neck cell carcinogenesis in mice. (A) Four experimental and control groups were examined. Mice in the experimental group were fed drinking water containing 4-NQO and arecoline and then monitored until being sacrificed. (B) The progression of precancerous lesions on the tongue to cancer was monitored when mice were alive, and they were sacrificed at different time points. With the increasing duration of treatment in the experimental group, white lesions began to appear on the tongues of mice. Tumors were observed on the tongues of mice at week 16 and had become invasive carcinoma by week 29. (C) The tongues of mice in the control and experimental groups underwent hematoxylin and eosin staining for histopathological examination. In this study, mouse tongues were excised for scRNA-seq analysis at 16 and 29 weeks. Each cell was wrapped in oil droplets to generate GEMs (Supplementary Materials Figure S1 ). Inside the GEMs, intracellular mRNA was reverse transcribed into cDNA, and Qubit fluorometric quantitation (Thermo Fisher Scientific) was used for the concentration measurement. The data were as follows: 0.162, 0.501, and 0.716 ng/µL for the 16-week control group; 0.126, 0.446, and 0.584 ng/µL for the 16-week experimental group; 0.132 and 3.100 ng/µL for the 29-week control group, and 0.288, 5.340, and 5.280 ng/µL for the 29-week experimental group. The required number of cycles in the library construction was according to its concentration. After library construction, the concentrations were 33.8, 19 .0, and 25.0 ng/µL for the 16-week control group; 32.0, 14.0, and 18.0 ng/µL for the 16-week experimental group; 21.0 and 14.0 ng/µL for the 29-week control group, and 37.6, 23.4, and 26.8 ng/µL for the 29-week experimental group. The 11 groups of samples were sequenced for next-generation sequencing analysis, and the results were mapped with the mouse genome mm10 and divided into 16-and 29-week control and experimental groups. All cell reads are presented in Table 1 . The total number of sequencing reads varied according to the number of cells packed each time, and the mouse genome (mm10) accounted for more than 92% of the map. The Cellranger (v3.1.0, 10x Genomics) was used to analyze the sequencing reads for each sample, and the results of mapping with the mouse reference genome (mm10). The results of all samples were aggregated, and corrections were made on the basis of the library size of each sample. Each cell had read counts of 31,053 genes, which were analyzed using the Seurat R package and normalized according to the total count for each cell. Genes with the greatest differences between cells were identified; 1000 genes were selected from 31,053 genes for subsequent analysis (Supplementary Materials Figure S2a ), and then, principal component analysis was performed to linearly reduce the dimensionality of 1000 gene data (Supplementary Materials Figure S2b ,c). The first 10 principal components were used to calculate clustering and dimensionality reduction and provide data visualization of t-distributed stochastic neighbor embedding (tSNE) (Figure 2A,B) . The cells were clustered into 17 subtypes. The cells in the four groups (16-and 29- week control and experimental groups) were plotted separately with tSNE ( Figure 2C ). The top 10 genes with the highest expression compared with other cell populations were identified and represented by a clustering heat map ( Figure 2D ). Most cell subtypes had a high expression of specific genes. The genes corresponding to each cell subtype in the cluster heat map are listed in Supplementary Materials Table S5 . After cells were separated into 17 subtypes, the 10 genes with the highest expression in each subtype of cells were identified and compared between subtypes. Yellow indicates that the gene expression is higher than in other cell subtypes. The tSNE map offered a comparison of the control and experimental groups with respect to the proportion of a subtype of cells within and between cell groups ( Table 2 ). The proportions for each cell subtype are displayed in a pie chart ( Figure 3A ). The proportions of the seventh and ninth cell subtypes in each group are displayed in Figure 3B . The proportion of the seventh cell subtype was lower in the 29-week control group than in the 16-week control group, and that of the ninth cell subtype was similar for both of these groups. The proportions of the seventh and ninth cell subtypes were higher in the 29-week experimental group than in the 16-week experimental group ( Figure 3C ). The seventh cell subtype was catalogued as stem cell and the ninth subtype was catalogued as keratinocyte by PanglaoDB [15] . The results revealed increased oncogenic activity in these two cell subtypes. Gene set enrichment analysis was performed to determine the regulatory pathways involved in carcinogenesis. The tSNE map revealed that the seventh and ninth cell subtypes exhibited an increasing trend in the proportion of cells in the experimental group; their differences in expression genes at 29 and 16 weeks in the experimental group are listed in Supplementary Materials Table S6 . Analysis of the molecular signature database (MSigDB) hallmark gene set collection was performed to identify three relevant regulatory pathways for the expression level of different genes in the seventh cell subtype in the 16-and 29-week experimental groups. The increasing expression pathways were MYC_targets_v1, Allograft_rejection, and Oxidative_phosphorylation, as represented by an enrichment plot ( Figure 3D ); the decreasing expression pathways were Heme_metabolism, Protein_secretion, and Mitotic_spindle ( Figure 3E ). The calculation method is based on the fold change of all genes in the 16-and 29-week experimental groups. A line graph of enrichment scores obtained by the cells of different genes was plotted for a comparison of the selected gene group data within the range between the maximum and minimum peak position; the gene groups at the left and right end points of the graph (leading-edge subset) are important genes involved in regulating the pathway. The normalized enrichment scores (NESs) of the regulatory pathways ranged between −2 and 3, and the calculated p value was significant ( Figure 3F ) (Supplementary Materials Table S7 ). In the ninth cell subtype, an enrichment plot was generated to display the top three related regulatory pathways with the greatest increase and the top three with the greatest decrease in gene expression levels between the 16-and 29-week experimental groups. Those with the greatest increase were MYC_targets_v1, Oxidative_ phosphorylation, and Unfolded_protein_response ( Figure 3G) ; those with the greatest decrease were KRAS_signaling_up, IL2_STAT5_signaling, and TNFα_signaling_via_NFkB ( Figure 3H ). The NESs of the first three regulatory pathways ranged between −3 and 3, and the p value was significant ( Figure 3I ) (Supplementary Materials Table S8 ). The gene expression clusters of the most significant regulatory pathways of the seventh and ninth cell subtypes in the 16-and 29-week experimental groups were MYC_targets_v1, as represented by incremental points in Figure 3J ,K. For the seventh and ninth cell subtypes, the average expression of the most highly expressed genes in the MYC_targets_v1 pathway were increased, and the proportion of cells that expressed these genes was also increased in the 29-week experimental group compared with the 16-week experimental group. The ratio of the average expression among cell expression of these genes exhibited a downward trend in the 29-week experimental group compared with the 16-week experimental group. The involvement of the genes RANBP1, MCM5, EIF3B, PSMA6, NPM1, and HSP90AB1 in the most significant regulatory pathways of the seventh and ninth cell subtypes was further confirmed in the nasopharyngeal carcinoma cell line (HONE-1), and cisplatin was used to screen for the expression of these genes in the drug-resistant HONE1-CIS6 cell line. RT-PCR results revealed that these genes were more highly expressed in the drug-resistant HONE1-CIS6 cell line than in the HONE-1 cell line ( Figure 4A) . Additionally, these genes had higher protein expression levels in the drug-resistant HONE1-CIS6 cell line than in the HONE-1 cell lines ( Figure 4B ) (Supplementary Materials Figure S4 ). Oral cancer is the sixth most common cancer among men worldwide [16] . In Taiwan, it ranks fifth in mortality among all cancers in the general population and fourth in the male population [17] . Surgery is currently the primary approach for treating oral cancer. For advanced stages, chemotherapy, radiation, and immune therapies are the primary treatment options. However, the efficacy of these treatments is limited when the diagnosis is made at a later stage [18] . Therefore, identifying new biomarkers for detecting precancerous lesions and preventing their progression to oral cancer are pivotal. The composition of cancer cells is diverse and complex in all human cancers. In the same tumor, cells may vary in terms of their genetic, epigenetic, and transcription characteristics and proteomic dysregulation, all of which affect the overall disease progression. If evaluations are based on the characteristics of total cells, then the result may be failure to apply a treatment modality that can effectively target the most oncogenic tumor cell subtype, thus considerably reducing the therapeutic effects and patients' likelihood of survival [12, 19] . For example, some cancer cell subtypes have a high degree of proliferation, and some subtypes are more invasive and metastatic. The entire tumor microenvironment is highly heterogeneous, and higher heterogeneity among cells corresponds to weaker therapeutic effects, resistance to disease treatments, and an increased risk of recurrence [12] . Traditional RNA expression analysis is often employed to study the characteristics of a group of cells by averaging the characteristics displayed by specific cell subtypes. The use of single-cell technology Oral cancer is the sixth most common cancer among men worldwide [16] . In Taiwan, it ranks fifth in mortality among all cancers in the general population and fourth in the male population [17] . Surgery is currently the primary approach for treating oral cancer. For advanced stages, chemotherapy, radiation, and immune therapies are the primary treatment options. However, the efficacy of these treatments is limited when the diagnosis is made at a later stage [18] . Therefore, identifying new biomarkers for detecting precancerous lesions and preventing their progression to oral cancer are pivotal. The composition of cancer cells is diverse and complex in all human cancers. In the same tumor, cells may vary in terms of their genetic, epigenetic, and transcription characteristics and proteomic dysregulation, all of which affect the overall disease progression. If evaluations are based on the characteristics of total cells, then the result may be failure to apply a treatment modality that can effectively target the most oncogenic tumor cell subtype, thus considerably reducing the therapeutic effects and patients' likelihood of survival [12, 19] . For example, some cancer cell subtypes have a high degree of proliferation, and some subtypes are more invasive and metastatic. The entire tumor microenvironment is highly heterogeneous, and higher heterogeneity among cells corresponds to weaker therapeutic effects, resistance to disease treatments, and an increased risk of recurrence [12] . Traditional RNA expression analysis is often employed to study the characteristics of a group of cells by averaging the characteristics displayed by specific cell subtypes. The use of single-cell technology provides a new approach for exploring the entire tumor microenvironment and clarifying the functions of individual cells. For triple negative breast cancer, which is the most aggressive and chemotherapy-resistant type of breast cancer, scRNA-seq can be used to identify specific cells that induce cancer formation and affect disease progression [20] . Single-cell technology can also be used to elucidate the mechanisms of cell resistance and genomic and transcription dysregulation [21] . In ovarian cancer, malignant cells and cancer-related fibroblasts express genes related to the JAK-STAT signaling pathway, and inhibiting this pathway was reported to promote effective anti-tumor ability in the primary cell cultures and xenograft models of patients [22] . Tumor-infiltrating lymphocytes are a potential benefit of immunotherapy and enhance clinical response to cancer, and analysis of the T-cell population in hepatoma tissues revealed that the LAYN gene in infiltrating T cells can inhibit the function of regulatory T cells and CD8 T cells [10] . The use of scRNA-seq technology can enable the further exploration of factors affecting the course of a disease and reveal the characteristics and control pathways of specific cell subtypes, thus facilitating disease prevention and detection and improving treatment strategies to prevent disease progression or recurrence. In this research, scRNA-seq was used to investigate the different stages of carcinogenesis in an oral cancer mouse model and revealed associations of oncogenic clusters 7 (stem cells) and 9 (keratinocytes) with the MYC_targets_v1 pathway. Similar to previous research [23, 24] , an experiment was conducted in the present study wherein 4-NQO, which is a soluble carcinogen in cigarettes, was added to the drinking water of mice in the experimental group. Carcinogens in water can disrupt DNA and induce the production of reactive oxygen species. The other carcinogen used in this study, arecoline, is the main alkaloid in betel nut. It is known to have cytotoxicity and genotoxicity and induce mutations, which can cause histological or other biological changes [25] . In mice, the combination of the two drugs can mimic the occurrence of oral cancer, ranging from epithelial cell hyperplasia and cell dysplasia to invasive squamous cell carcinoma [14] . H&E staining and images captured every week confirmed the carcinogenesis of the mice tongues. The findings regarding the carcinogenesis observed in the male C57BL/6J mouse model in this study can also serve as a reference for minimizing the effects of variance among individual patients during analysis of carcinogenesis in different cell subtypes. In this study, 10x Chromium single-cell technology was used to extract healthy-tongue and tumor cells in the 16-and 29-week control and experimental groups for single-cell sorting. The factors of carcinogenesis that caused the cells to become specific malignant cell subtypes were analyzed in four groups (16-and 29- week control and experimental groups). The cells were divided into 17 subtypes, and two cell subtypes with little difference in the proportion or expression level between the 16and 29-week control groups were selected; the cell proportion and expression exhibited significant upregulation in the 16-and 29-week experimental groups. These two cell subtypes were inferred to be influential oncogenic factors in carcinogenesis. Gene function enrichment analysis [26, 27] revealed that the most significant regulatory pathway involved in the two cell subtypes is MYC_targets_v1. MYC is not only a regulatory gene that affects gene expression but also a proto-oncogene. Its three related genes, c-myc (MYC), l-myc (MYCL), and n-myc (MYCN), have been found to be involved in many human cancers, with the regulation of some signal pathways leading to carcinogenesis [28] . Additionally, c-myc is involved in some physiological functions, including cell cycle regulation, ribosome biogenesis, protein synthesis, stem cell function regulation, and cell apoptosis [29] . MYC mutations have been related to Burkitt's lymphoma and leukemia [30, 31] ; they are activated through gene rearrangement and gene amplification, which affect other genes. In addition, the overexpression of the MYC gene is also found in breast, cervical, colorectal, lung, and gastric cancers and leads to cell proliferation [32] [33] [34] [35] [36] . The overexpression of the MYC gene has also been noted in most patients with oral cancer. For patients with T3 or T4 oral cancer, MYC gene expression is related to the likelihood of survival and can serve as a key biomarker for treatment [37] . In the seventh cell subtype (stem cells) involved in the MYC_targets_v1 pathway, increased RANBP1, MCM5, and EIF3B gene expression was noted in the 29-week experimental group; the average expression level and cell proportion noticeably increased. The regulation of c-Myc mRNA is accompanied by specific RNA-binding proteins and also occurs at the translation level. For example, CELF1 inhibits Myc translation by reducing the binding of c-Myc mRNA to HuR, and it provides new insights into the molecular function of RNA-binding proteins [38] . The interaction between RNA-binding proteins and c-Myc mRNA may become a new target for cancer treatment [39] . MCM5 is a DNA replication licensing factor. The protein is upregulated during the transition from G0 to G1/S of the cell cycle, and it is involved in numerous MYC regulatory pathways related to DNA replication [40] . EIF3B is a member of the eukaryotic translation initiation factor 3 subunit B and is related to the transcription initiation of cancer-related genes. Gene function enrichment analysis revealed that EIF3B is involved in the regulatory pathway and overexpressed in patients with hepatoma [41] . The ninth cell subtype (keratinocytes) had a larger proportion of cells and increased PSMA6, NPM1, and HSP90A expression in the 29-week experimental group. PSMA6 is a gene encoding the proteasome α6 subunit protein, which can assemble the 20S proteasome complex, and the proteasome is a key component of the ubiquitin-proteasome system [42] . This system has an important function in malignant transformation and regulates gene expression by degrading transcription factors, such as c-Myc. The ubiquitin-proteasome system participates in the development of many malignant cancers [43] . Nucleophosmin is a highly expressed protein located in the nucleolus that moves between the nucleus and cytoplasm. Its main functions are ribosome maturation and export, participation in centrosome replication, cell cycle progression, histone assembly, and response to various stimuli. In acute myeloid leukemia and some solid tumors, NPM1 gene overexpression is related to mitotic index and metastasis [44] . The Cancer Genome Atlas database was used to analyze the relationship between NPM1 mRNA expression and the survival rate of patients with head and neck cancer. High expression of NPM1 is related to a lower likelihood of survival (Supplementary Materials Figure S3 ). Nucleophosmin is crucial to the oncogenic activity of c-Myc. Overexpression of nucleophosmin stimulates c-Myc to induce hyperproliferation and transformation [45] . Subtypes of heat shock protein 90 include HSP90AA (encoded by HSP90AA1 and HSP90AA2) and HSP90AB (encoded by HSP90AB1). Studies have demonstrated that c-Myc directly activates the transcription of HSP90A, and this transcriptional induction occurs in different tissues and is not related to cell proliferation. Since the c-Myc binding site is near the HSP90A gene promoter, c-Myc induces HSP90A expression and can control many message transmission pathways, including cell transformation [46] . The drug-resistant human nasopharyngeal carcinoma cell line HONE1-CIS6 was originally detected in HONE-1 through the use of 10 µM cisplatin in culture medium. HONE1-CIS6 could represent a cisplatin-resistant cell subtype of HONE-1; such subtypes are more aggressive and are associated with higher rates of disease recurrence. RANBP1, MCM5, EIF3B, PSMA6, NPM1, and HSP90AB1 were identified through RT-PCR analysis of HONE-1 and HONE1-CIS6 as the genes most significantly involved in the MYC_targets_v1 regulatory pathway in the seventh and ninth cell subtypes; moreover, MCM5, RANBP1, and HSP90AB1 were significantly overexpressed in HONE1-CIS6. Male C57BL/6J mice raised on the Tainan campus of the National Institutes of Health, Taiwan, were used in experiments in this study in compliance with the guidelines of the Laboratory Animal Care and Use Committee of National Cheng Kung University. Mice were randomly divided into a control group and experimental group: the control group received normal drinking water, and the experimental group received drinking water with 4-NQO (200 µg/mL) (SIGMA-ALDRICH, N8141-5G) and arecoline (200 µg/mL) (TCI, A0523). All mice in the experimental group, except for those sacrificed at week 8, were treated for 16 weeks; then, they were either sacrificed or given normal drinking water thereafter until being sacrificed at weeks 20 or 29. Observation and analysis confirmed that the treatment stimulated healthy oral mucosa cells in mice to evolve into precancerous lesions and then into cancerous and malignant tissues. After mice were sacrificed through carbon dioxide asphyxiation, incisional biopsy was performed on lesions on their tongues, and the specimens obtained were fixed in 10% buffered formalin. The specimens were stained with hematoxylin and eosin (H&E) for histopathological examination, and whole lesions were excised for subsequent scRNA-seq analysis. Tongue lesions were rinsed with 1× phosphate-buffered saline (PBS) (VWR, 97062-730), cut with a scalpel and soaked in gentleMACS™ C tubes (MACS, 5190703020) containing cell culture medium (DMEM/F12; Simply, CC131-1000) and tissue dissociation enzyme (neutral protease Dispase II (100 U/MG; SIGMA, D4693-1G) and Collagenase IV (0.1 mg/mL; SIGMA, C5138-100 MG)). Then, they were placed in an automatic tissue homogenization/single-cell separator (gentleMACSTM Octo Dissociator with Heaters) and maintained at 37 • C for 10 min to separate tissues into single cells. Subsequently, the lesions were placed in an incubator and maintained at 5% CO 2 and a temperature of 37 • C for 3 min and then filtered through a 40-µm cell strainer (Falcon, 352340) . The cell-containing suspension was collected, and 1× red blood cell lysis buffer (10× red blood cell lysis buffer; BD, 555899) was added. After lysis of red blood cells, the suspension was centrifuged at 1500 rpm for 5 min to remove the supernatant, washed with PBS containing 0.04% bovine serum albumin (Cyrusbioscience, 101-9048-46-8), and then centrifuged in PBS. Next, 10 µL of the final cell suspension was mixed with 10 µL of trypan blue (VWR, 72-57-1), and an automatic cell counter was used to determine the number of cells. Finally, the concentration was adjusted to 1000 cells/µL for microfluidic cell sorting. 10× Genomics Platform and 3 Gene Expression Library Construction [47] were performed follow the manufactural recommendation protocols and list in Supplementary Materials Table S1. The 3 gene expression library, constructed as described in the previous section, was subjected to next-generation sequencing with the Illumina NextSeq 500/550 High-Output v2.5 Kit (150 cycles) (20024907) and sample indices. The parameters for paired-end, single indexing sequencing were as follows: Read 1, 28 cycles; i7 Index, 8 cycles; Read 2, 91 cycles; and Library Loading NextSeq 500/550, 1.5 pM [48] . After next-generation sequencing, Cellranger (v3.1.0, 10x Genomics) was used to analyze sequencing reads [49] and perform comparisons with original data for the mouse reference genome (mm10) by using STAR (RRID:SCR_005622) [50]-where "Count" is used for reads alignment, filtering, barcode counting, MMI counting, cell clustering, and gene expression analysis with the cell barcode and "aggr" is used to aggregate sample results [51] -and output the cloupe file to the Loupe Cell Browser (10x Genomics, v.3.1.0) [52] . The cell lines of the human nasopharyngeal carcinoma cell line (HONE-1) and drug-resistant human nasopharyngeal carcinoma cell line (HONE1-CIS6) screened using 10 µM cisplatin (Hospira Australia Pty Ltd.) were employed in the present study [53] . These two cell lines were provided by the Laboratory of Distinguished Professor Dar-Bin Shieh, Institute of Oral Medicine, National Cheng Kung University, Tainan, Taiwan. The cell lines were cultured directly from liquid nitrogen storage cells that were validated on 12/10/2013 (Supplementary Materials Figures S5 and S6) . The cell lines were cultured with an RPMI 1640 cell culture medium (Gibco TM , 31800-014) containing 2 mM L-glutamine (Simply, CC515-0100) with the addition of 2 g/L sodium bicarbonate (Millipore, 144-55-8), 10% fetal bovine serum (Gibco, A38401), 1% L-glutamine (Simply, CC515-0100), and 1% antibiotic/antimycotic (Simply, CC501-0100). For the drug-resistant human nasopharyngeal carcinoma cell line (HONE1-CIS6), 10 µM cisplatin was added to the cell culture medium to preserve the cell characteristics. RNA Expression Assay details were listed in Supplementary Materials Table S3 . Protein Expression Assay details were listed in Supplementary Materials Figure S4 . In conclusion, a mouse model of oral tongue cell carcinogenesis induced by 4NQO and areocline administration enabled the transition of cells from healthy to precancerous and finally to cancerous to be studied with minimal individual heterogeneity. Tongue lesions obtained from the 16-and 29-week experimental groups were analyzed through scRNA-seq and compared with those obtained from the control group; 17 cell subtypes involved in carcinogenesis were identified, and the seventh (stem cell) and ninth (keratinocyte) oncogenic subtypes had the most significant involvement in the pathway related to MYC. The results were further validated by the presence of the cisplatin-resistant cell line HONE1-CIS6 with significant overexpression of MCM5, RANBP1, and HSP90AB1. These cell subtype biomarkers in oral cancer carcinogenesis can be applied in the detection of patients with precancerous lesions in the identification of high-risk populations and as a treatment target. Liquid Biopsy in Oral Cancer Oral cancer: Etiology and risk factors: A review Identification of Potential Biomarkers and Survival Analysis for Head and Neck Squamous Cell Carcinoma Using Bioinformatics Strategy: A Study Based on TCGA and GEO Datasets Gene Expression Subtype Predicts Nodal Metastasis and Survival in Human Papillomavirus-Negative Head and Neck Cancer Single-Cell RNA Sequencing in Cancer: Lessons Learned and Emerging Challenges Platforms for Single-Cell Collection and Analysis Droplet-based single cell RNAseq tools: A practical guide Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq Landscape of Infiltrating T Cells in Liver Cancer Revealed by Single-Cell Sequencing Single cell sequencing reveals heterogeneity within ovarian cancer epithelium and cancer associated stromal cells Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer Co-treating with arecoline and 4-nitroquinoline 1-oxide to establish a mouse model mimicking oral tumorigenesis A web server for exploration of mouse and human single-cell RNA sequencing data Global cancer statistics Relationship of Oxidative Stress, Inflammation, and the Risk of Metabolic Syndrome in Patients with Oral Cancer Controlled Drug Delivery Systems for Oral Cancer Treatment-Current Status and Future Perspectives Evolution of the cancer stem cell model Single-cell RNA-sequencing: The future of genome biology is now Chemoresistance Evolution in Triple-Negative Breast Cancer Delineated by Single-Cell Sequencing A single-cell landscape of high-grade serous ovarian cancer The increase of oncogenic miRNA expression in tongue carcinogenesis of a mouse model K14-EGFP-miR-31 transgenic mice have high susceptibility to chemical-induced squamous cell tumorigenesis that is associating with Ku80 repression Current Concepts about Areca Nut Chewing Molecular signatures database (MSigDB) 3.0 The Molecular Signatures Database (MSigDB) hallmark gene set collection Tristetraprolin impairs myc-induced lymphoma and abolishes the malignant state Small molecules targeting c-Myc oncogene: Promising anti-cancer therapeutics Myc suppression of Nfkb2 accelerates lymphomagenesis Deep sequencing of MYC DNA-binding sites in Burkitt lymphoma MYC-Driven Pathways in Breast Cancer Subtypes GDF15 promotes the proliferation of cervical cancer cells by phosphorylating AKT1 and Erk1/2 through the receptor ErbB2 Modeling oncogenic signaling in colon tumors by multidirectional analyses of microarray data directed for maximization of analytical reliability MYC Drives Progression of Small Cell Lung Cancer to a Variant Neuroendocrine Subtype with Vulnerability to Aurora Kinase Inhibition YAP/TAZ Initiates Gastric Tumorigenesis via Upregulation of MYC c-MYC expression in T (III/IV) stage oral squamous cell carcinoma (OSCC) patients Competition between RNA-binding proteins CELF1 and HuR modulates MYC translation and intestinal epithelium renewal Targeting MYC in cancer therapy: RNA processing offers new opportunities MYC: Connecting selective transcriptional control to global RNA production Expression of eukaryotic translation initiation factor 3 subunit B in liver cancer and its prognostic significance Perilous journey: A tour of the ubiquitin-proteasome system Potential for proteasome inhibition in the treatment of cancer Molecules that target nucleophosmin for cancer treatment: An update Nucleophosmin interacts directly with c-Myc and controls c-Myc-induced hyperproliferation and transformation Direct activation of HSP90A transcription by c-Myc contributes to c-Myc-induced transformation Chromium 10x Single-Cell 3' mRNA Sequencing of Tumor-Infiltrating Lymphocytes RNA Sequencing in B-Cell Lymphomas Immune cell profiling of COVID-19 patients in the recovery stage by single-cell sequencing Ultrafast universal RNA-seq aligner SCelVis: Exploratory single cell data analysis on the desktop and in the cloud Leveraging Single-Cell RNA Sequencing Experiments to Model Intratumor Heterogeneity Gelsolin regulates cisplatin sensitivity in human head-and-neck cancer The authors declare no conflict of interest.