Abstract
Head and neck squamous cell carcinoma (HNSCC) is a common malignancy with poor prognosis and immune response, which plays an important role in tumor progression. Recently, immunotherapies have revolutionized the therapeutic means of malignancies including HNSCC. However, the relationship between immunophenotypes of HNSCC and its clinical response to immune-checkpoint inhibitors remains unclear. We aim to identify molecular subtyping related to distinct immunophenotypes in HNSCC. Consensus clustering algorithm was conducted for subtyping. Immunophenotypes between subtypes were compared according to infiltrating immunocytes, immune reactions, major histocompatibility complex (MHC) family, immunoinhibitory, immunostimulatory and immune scores. The relationship between immunophenotype and genotype was investigated from gene mutation and tumor mutation burden. The potential response of Immune-checkpoint blockade (ICB) therapy was estimated with TIDE and ImmuCellAI algorithms, and immune-checkpoint genes. The immune characteristics were also investigated. Biological functions were annotated by the gene-set enrichment analysis (GSEA) algorithm. Two distinct immune subtypes of HNSCC with different survival outcomes, biological characteristics, immunophenotype, and ICB response were identified. The subtype-1 was featured with better prognosis, more infiltrated immunocytes, stronger immune reaction, higher immune-related gene expression, higher immune-checkpoint gene expression (PD-1, PD-L1, and CTLA-4), and better ICB response. A higher immune response in subtype-1 was also revealed by GSEA. Subtype-1 possessed a higher immune response and more sensitivity to ICB therapy leading to a better prognosis. These findings may shed promising light on the immunotherapy strategy in HNSCC
Keywords
Introduction
Head and neck squamous cell carcinoma (HNSCC) is the most common head and neck malignancies.1,2 It is one of the top 10 most common malignancies in the world, with about 300,000 new cases and 140,000 deaths worldwide every year. 3 HNSCC has multiple adverse clinical characteristics such as strong local invasion, early lymph node metastasis, poor prognosis, etc. Despite the recent development in the diagnosis and treatment of HNSCC, the 5-year survival rate of patients remains only approximately 50%. 4 To date, the treatment strategies of HNSCC are mainly surgical resection, assisted with radiotherapy and chemotherapy. For patients with advanced and recurrent tumors, the survival rate is still very low despite the use of comprehensive treatment. It has been the unremitting pursuit for researchers to find new treatment methods in recent years.
More and more studies reported that the abnormal function of the immune system plays a crucial part in tumor development. HNSCC is an immunosuppressive malignancy and in its tumor microenvironment, there occurs lymphocytes dysregulation,5,6 incapacity of effector T cell,7,8 and weakened antigen presentation effect.9,10 Meanwhile, there were a large number of immunosuppressor cells, such as tumor-associated macrophages, myeloid-derived suppressor cells (MDSCs) and Tregs, which inhibited the tumor-killing effect of the immune system. 11 Immune suppression and immune escape caused by the abnormal function of immunocytes are closely associated with the occurrence and progression of HNSCC. 11 With the discovery of tumor immune checkpoints and their important role in abnormal immune function, immunotherapy targeting tumor immune checkpoint has become a heated topic in recent years. Targeting immune checkpoint PD-1 and CTLA-4 immunotherapy has made an effective clinical treatment effect and improved the prognosis of patients with tumors..10,12,13
However, major limitations can be found in immune-checkpoint blockade (ICB) therapy, that is its only effective in one-third of patients in most tumors. 14 Some patients have obvious immunotherapy responses and some do not, so there are significant individual differences in the clinical results of immunotherapy. 15 The tumor immune microenvironment is one of the factors that greatly influence the effect of immunotherapy, and the complexity of the tumor immune microenvironment makes it harder for immunotherapy to exert its effectiveness. 16
Hence, it is necessary to further investigate the patients’ overall immune status in order to identify the molecular subtypes of HNSCC thus improving the immunotherapy effect of patients with advanced HNSCC. We aimed to use a series of bioinformatic approaches to identify HNSCC immune subtypes and investigated their underlying immune characteristics that detects potentially targetable ICB subtypes
Materials and Methods
Data Acquisition
Since the primary solid tumor was the main study purpose for HNSCC researched in TCGA-HNSCC cohort, we excluded other cancerous tissue such as paracancerous tissue and metastatic carcinoma and retained only the primary solid tumor so as to avoid confounding factors, and in addition, samples missing prognosis information were also excluded. In this way, 494 HNSCC samples were involved in immune subtyping analysis. As for clinical feature adjusting analysis, only samples with complete clinical features were involved, which is 396 samples. The RNA-seq data and clinical information data were downloaded from the online database UCSC Xena (https://xenabrowser.net/datapages/). Since all of the data involved in our analysis are available to the public in the TCGA databases, ethical approval of the local institution was not needed for this study. All processes in the present study followed the relevant guidelines and policies of the TCGA databases. Immune-related gene lists were obtained from ImmPort (https://www.immport.org/shared/genelists) with 17 category immune genes. The somatic mutation (mutects) data were downloaded by R package “TCGAbiolinks.” All data were processed according to standard procedures including low expression gene removing, gene name mapping, expression quantile normalization, log2 conversion, missing values removing. 17 The tumor stemness indices were obtained from the previous study's supplementary file, in which whole TCGA tumors’ stemness indices were calculated. 18 The immune-checkpoint blockade response scores were obtained from the previous study, in which all TCGA tumors’ score has been calculated in the supplementary file. 14
Immune-Related Subtype Identification
The tumor microenvironment (TME), which includes blood vessels, lymph vessels, immune cells and mesenchymal cells, is a complex ecosystem of stromal cells and plays a critical role during tumorigenesis and progression. 19 The TME is a major contributing factor for patient outcomes and immune-checkpoint blockade (ICB) treatment. Patients harboring different immune subtypes results in distinct TME status, including major differences in the infiltrating immunocytes, and immune reactions surrounding the cancerous lesion, which no doubt could lead to different biological reactions and responses to drug regimen. Therefore, patients with specific immune characteristics and specific outcomes in the HNSCC are referred to as distinct immune subtypes, because their subgrouping is based on their immune characteristics. Immune genes were analyzed by univariate cox regression and genes with false discovery rates (FDR)<0.05 were regarded as survival-related immune genes. Prognosis-related immune genes were identified and a gene expression matrix that contained prognostically significant immune-related genes was formed. Unsupervised cluster technique was use for subtyping and unsupervised class discovery is a data mining technique for the detection of unknown possible groups of items based on intrinsic features and no external information. Consensus clustering a widely used class discovery tool with confidence assessments and item tracking by R package “ConsensusClusterPlus.” 20 The prognostically significant immune-related genes matrix was applied to consensus clustering and the parameter set as k-means cluster algorithm, euclidean distance, 1000 subsamples and 0.8 proportion of items to sample. The optimal number of subtypes was selected according to Delta area, consensus cumulative distribution function (CDF) and consensus matrix. We chose the ‘cleanest’ cluster partition according to the consensus matrix, where items nearly always either cluster together giving a high consensus (dark blue color) or do not cluster together giving a low consensus (white) (Figure 1D). We also chose the k according to CDF and Delta area, at which the distribution reaches an approximate maximum, which indicates a maximum stability and after which divisions are equivalent to random picks rather than true cluster structure (Figure 1B and C). Therefore, we determined 2 clusters/subtypes patients that distinct with each other in the expression pattern of prognostically significant immune-related genes. We named the 2 clusters/subtypes patients as immune subtype-1 and immune subtype-2, because they two are clustered by immune genes and also have distinct characteristics. Multivariate cox regression was performed to analyze immune subtypes and clinical features to explore the independence of immune subtypes.

Identification of HNSCC immune-related subtypes. (A) The top 20 immune-related genes closely associated with HNSCC overall survival were presented by forest plot. (B) Relative change in area under the CDF curve for k = 2 to 7. (C) Consensus clustering cumulative distribution function (CDF) for k = 2 to 7. (D) Heatmap of the matrix of co-occurrence proportions HNSCC samples. (E) Principal component analysis (PCA) of 266 survival-related immune genes between two HNSCC immune subtypes. The two first principal components (PC1, PC2) which explain most of the variables are plotted.
Immune Characteristics Analysis
The relative enrichment of infiltrating immunocytes and the activity of immune pathways were determined by single-sample gene-set enrichment analysis (ssGSEA).21,22 The gene sets used to determine the composition of immunocytes were obtained from the previous study, 23 which included 28 widely used tumor infiltrated immunocytes such as T cell family and B cell family, and we only presented the statistically significant immunocytes in the results figures. The gene sets used to evaluate the activity of immune-related pathways were obtained from the database ImmPort (http://www.immport.org), 24 which contains 17 immunologically relevant lists of genes curated with functions and Gene Ontology terms. The MHC genes, immunoinhibitory and immunostimulatory list were obtained from previous study supplementary file, 25 which are widely used three classes of molecules that are involved in tumor escape mechanisms. ImmuneScore and StromalScore were calculated by R package “estimate,” which is widely used to estimate the purity of carcinomas or immunocyte abundance. The ICB response and benefits prediction was performed by ImmuCellAI algorithm and TIDE algorithm14,26 to predict the response of ICB therapy.
Functional Enrichment Analysis
Functional enrichment analysis was conducted for annotating immune subtype biological characteristics by R package “ClusterProfiler” from 3 aspects including hallmark gene-sets, GO biological processes gene-sets and canonical pathways gene-sets. The marker genes for the two immune subtypes were identified by differentially expressed gene (DEG) analysis.
Mutation Landscape Analysis
The GDC DNA-Seq analysis pipeline identifies somatic variants within whole exome sequencing (WXS) and whole genome sequencing (WGS) data. Somatic variants are identified by comparing allele frequencies in normal and tumor sample alignments, annotating each mutation, and aggregating mutations from multiple cases into one project file. And what we used is the widely used MuTect2 Variant Aggregation and Masking pipeline. 27 Somatic mutation genes were analyzed by R package “maftools” and Frame Shift Del, Missense Mutation, Nonsense Mutation, Multi Hit, Frame Shift Ins, In Frame Ins, Splice Site, In Frame Del, Translation Start Site were considered as mutation. The top 20 mutation frequency of mutated genes in the TCGA cohort were presented. And the 20 genes mutation status between two immune subtypes were analyzed. Tumor mutation burden (TMB) is defined as the number of nonsynonymous alterations (SNVs or indels). 38 Mb was used as the estimation of the exome size for calculating the TMB. 28
Statistical Analysis and Visualization
R software 3.6.1 was used in this study for statistical analysis and visualization. R package “ggplot2,” “ggpuber” were used to make statistical plots including box-plots, violin-plots volcano-plots. Differences of continuous variables between two groups were analyzed by a two-sided Wilcoxon rank-sum test. OS is the length of time to death or any cause (lost patients are the last time of follow-up; Patients who were still alive at the end of the study, the end date of follow-up), and is evaluated by Kaplan–Meier analysis that is one of the best options to be used to measure the fraction of subjects living for a certain amount of time. The survival curve was plotted with the Kaplan–Meier method by using the R package “survminer,” and the two-sided log-rank test was employed for evaluating the difference between the two groups. The p-value was adjusted by the FDR method for multiple hypothesis testing. Dichotomous variables were compared using the chi-square test.
Results
Identification of HNSCC Immune Subtypes With Prognosis Significance
A list of 1674 immune genes were involved in univariate cox analysis and 266 genes were significantly related to overall survival (FDR<0.05; Supplemental Table S1). The top 20 significantly prognosis-related immune genes were presented (Figure 1A). Using the 266 immune genes, 2 distinct HNSCC immune subtypes were identified using the Consensus clustering algorithm (Figure 1B to D). To confirm the variation in immune genes might be an intrinsic feature that could characterize the individual differences, we performed principal component analysis (PCA) analysis and found that the expression of prognosis-related immune genes displayed distinct group-bias clustering and individual differences (Figure 1E). The two subtypes of HNSCC were significantly different in overall survival as were demonstrated by Kaplan–Meier (KM) analysis, in which subtype-1 had better prognosis and subtype-2 suffered from a worse one (Figure 2A; HR: 1.56, 95% CI: 1.19-2.06, Log-rank p < 0.05). After being adjusted by clinical features including age, gender, tumor location, margin status, grade and stage, immune subtype was still a factor significantly associated with prognosis, suggesting it is an independent factor for prognosis (Figure 2B and Table 1; HR: 1.54, 95%CI: 1.12−2.12. p < 0.05).

Prognosis significance of the two HNSCC immune subtypes. (A) The Kaplan–Meier estimates of the overall survival for subtype-1 and subtype-2. The overall survival differences between the two groups were determined by the two-sided log-rank test. It can be concluded that subtype-2 is significantly associated with worse overall survival (p <0.05). (B) Multivariate Cox regression analysis depicted the association between overall survival and HNSCC immune subtypes adjusted with other clinical factors when taking into account and represented by forest-plot, indicating immune subtype is an independent factor associated with prognosis.
The Clinical Information for the HNSCC Cohort Used in the Study.
The Relationship Between Immune Characteristics and Immune Subtypes
To explore the immune characteristics between the two subtypes, the ssGSEA algorithm was conducted to estimate the immunocyte and immune reaction scores. Compared with subtype-2, higher infiltration status of almost all of the immunocytes could be found in subtype-1, except for CD56 bright natural killer cells (Figure 3A). And similar results were observed in the immune pathways, except TGFb_Family_Member and Interferons immune pathway, all of the other immune pathways were more active in subtype-1 (Figure 3B). These results suggested subtype-1 possessed higher infiltration of immunocytes infiltrated and higher activities of immune reactions, which is an immune-active subtype. As for immune-related genes including the MHC family, immunoinhibitory and immunostimulatory gene expression, most of them have a higher expression level in subtype-1 (Figure 4A to D). This also validates that subtype-1 is an immune-active cluster.

The distinct immune characteristics between two immune subtypes. (A) The comparison of infiltrating immunocyte abundance between two immune subtypes. (B) The comparison of immune reaction activity between two immune subtypes.

The immunophenotypes are different between the two HNSCC immune subtypes. (A) The comparison of MHC gene expression between two subtypes. (B) The comparison of immunoinhibitory gene expression between two subtypes. (C and D) The comparison of immunostimulatory gene expression between two subtypes.
Mutation Pattern and Clinical Features Related to Immune Subtypes
To explore the relationship between clinical features and immune subtypes, the distribution of clinical features in two subtypes was investigated. We found only tumor location is associated with the immune subtypes while there were no significant correlations concerning other features (Figure 5A). The oral cavity site is associated with subtype-2. Then, the mutation pattern was also investigated. The top 20 mutated genes in HNSCC, TP53, TTN, CDKN2A, CSMD3, and PCLO were significantly different in two subtypes (Figure 5B). There was no significant difference in TMB between the two immune subtypes (Figure 5C). The commonly used immune-checkpoint molecules expression of PD-1, PD-L1, and CTLA-4 were also investigated between the two immune subtypes. Subtype-1 has a significantly higher expression level of PD-1, PD-L1, and CTLA-4 than subtype-2, suggesting it may have a better response of ICB (Figure 5D to F). Besides, the ImmuneScore and StromalScore of subtype-1 are also higher than subtype-2 (Figure 5G to H).

Associations of HNSCC immune subtypes with immune and carcinoma characteristics. (A) Comparing of age, gender, tumor location, marginal status, grade and stage between immune subtypes. The heatmap illustrates the association of different clinical characters with subtype-1 and subtype-2 patients. (B) Mutational landscape in TCGA-HNSCC cohort stratified by immune subtypes and presented by oncoplot. The distributions of the top 20 mutated genes were compared between subtype-1 and subtype-2. (C) Distribution and association of mutation burden in the subtype-1 and subtype-2 groups. (D to F) Associations between immune subtypes and known immune-checkpoint genes. PD1, CTLA-4 and CD274 relative mRNA expression were shown between two immune subtypes. (G and H) The distribution of immune and stromal scores between two immune subtypes and immune subtypes was significantly associated with the scores.
Association Between Immune-Checkpoint Blockade Response and Immune Subtypes
We found more sutype-1 patients are belonging to atypical and mesenchymal TCGA subtypes and more subtype-2 patients are basal and classical TCGA subtypes (Figure 6B). To estimate the ICB response of the two immune subtypes, two famous algorithms were conducted to predict the ICB response. ImmuCellAI indicated there are more ICB response patients in subtype-1 and TIDE revealed patients of subtype-1 can obtain more benefits from ICB therapy than subtype-2 (Figure 6C and D). Subtype-2 has higher TIDE scores than subtype-1, suggesting patients of subtype-2 have a higher chance of antitumor immune escape and exhibiting a lower response rate of ICB treatment (Figure 6E). To uncover the different ICB responses in two subtypes, the induction of T cell dysfunction in tumors with high infiltration of cytotoxic T lymphocytes (CTL) and the prevention of T cell infiltration in tumors with low CTL level were calculated. The T cell exclusion score in subtype-1 is lower than subtype-2 and the T cell dysfunction score of subtype-1 was significantly higher than subtype-2, suggesting subtype −1 may be more sensitive to ICB therapy when CTL is low or subtype-2 may be sensitive when CTL is high (Figure 6F and G). These results imply subtype-1 has a better ICB therapy response.

The relationships between the novel immune subtypes and ICB response. (A) The relationship between the novel immune subtype and Thorsson's immune subtyping. (B) The relationship between the novel immune subtype and TCGA based HNSCC subtype. (C) The difference of ICB response between two immune subtypes. (D) The difference of ICB benefits between two immune subtypes. (E to G) Tumor Immune Dysfunction and Exclusion algorithm evaluates the factors that underlie mechanisms of tumor immune escape including TIDE scores, dysfunction scores and exclusion scores.
Gene and Biological Function Characteristics of the Two Subtypes
DEG analysis revealed IGHA1, IGKV4 to 1, IGKC, IGKV1 to 5, IGLV3 to 21 were the up-regulated genes in subtype-1 most and DES, MMP10, AREG, KRT14, and KLK5 were the most highly expressed genes in subtype-2 (Figure 7A and B, Supplemental Table S2). They are the marker genes for the two subtypes and representing the subtype molecular features. GO analysis revealed the functions of subtype markers, which suggested subtype-1 mainly is characterized by more immunophenotypes such as leukocyte migration, cell surface receptors action of immune response and B cell-related immune process (Figure 7C and D). GSEA was conducted to address the biological function diversity of the two subtypes. Immunity and immunocyte regulation-related pathways were frequently enriched in subtype-1 suggesting subtype-1 has a more active immune response (Figure 7E to G).

Biological functional characteristics of immune subtypes. (A) Differentially expressed gene analysis revealed the marker genes for immune subtypes and the volcano plot demonstrated the genes specific to immune subtypes. (B) The expression status of immune subtype-specific marker genes. (C and D) The GO enrichment analysis for immune subtype-related genes for subtype-1 and subtype-2. (E to G) Gene-set enrichment analysis from three aspects (including hallmark gene sets, GO biological processes and canonical pathways) revealed the underlying biomolecular differences between the two immune subtypes, representing their biological features.
Discussion
Tumor immunotherapy has been regarded as a promising therapeutic strategy, and among the many therapies aimed at improving the body's antitumor immunity, immune-checkpoint therapy is considered to be the most promising. 13 In the tumor microenvironment, the immune-checkpoint pathway can suppress the antitumor effect of the immune system, to aid the immune escape of tumor. 29 The immune-checkpoint pathway works by binding ligands to receptors on T cells, such as CTLA-4 to CD84 and CD80, and PD-1 to PD-L1. By blocking this binding, monoclonal antibodies can improve antitumor immunity and kill tumor cells. Pilimumab, a monoclonal antibody against CTLA-4, is the first immunocheckpoint targeted drug approved by the US Food And Drug Administration (FDA) for clinical use in patients with metastatic malignant melanoma. 30 This was followed by the development of Nivolumab, an anti-PD-1 monoclonal antibody. These targeted drugs achieved amazing results in the treatment of tumors,31,32 also including HNSCC. 33 However, some patients have obvious immunotherapy responses and some do not, so there are significant individual differences in the clinical results of immunotherapy. 15 Hence, it is necessary to further investigate the patients’ overall immune status for identifying the molecular subtypes of HNSCC and improve the immunotherapy effect of patients with advanced HNSCC.
Recently, increasing studies have focused on investigating the molecular subtypes of HNSCC based on multi-omics or genome-wide profiles to promote the realization of personalized treatment and improve the survival rate of patients.34–36 However, the realization of molecular subtypes of immunity in HNSCC remains in the beginning phase, and the other types of malignant have been well studied. Yang et al. have identified an FL11 prostate cancer subtype using immune infiltration profiles as well as transcriptomic data, which is highly enriched in immune system processes, immune-related KEGG pathways and biological processes and showing the promising immune therapy targeting prostate cancer. 37 Chen et al. used single-cell transcriptome to reveal underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma, which is very enlightening for us and reminds us of our further deeper studies. 38 Sylvie Job et al. have investigated the heterogeneous tumor stroma composition and built a TME-based classification of intrahepatic cholangiocarcinoma (ICC) that detects potentially targetable ICC subtypes, which is a great reference to us for our HNSCC immune subtypes clinical practice. 39 Therefore, to well address the problem of HNSCC immune subtypes, we conducted several bioinformatics methods to systematically analyze the immunity in HNSCC to identify immune-related HNSCC subtypes. We did the following work and made relevant discoveries. First, identified two immune-related HNSCC subtypes which are closely associated with the prognostic outcomes of patients, and the immune subtype is an independent factor that affects overall survival. Second, the immunotherapy significance between the two subtypes was investigated from ICB response prediction and immune-checkpoint gene expression. We found subtype-1 has better ICB results. Third, the diversity of the tumor microenvironment was revealed from multiple immune characteristics, and results suggested subtype-1 is an immune-active subtype. Forth, the relationship of mutation pattern was also considered and we found some famous mutated TP53 and CDKN2A seem associated with subtype-2. Last, the biological characteristics of immune subtypes were analyzed and subtype-1 has more active immune actions than subtype-2. Besides, we explore the association between the novel immune subtypes and previously identified subtypes. There is no connection with Thorsson’s immune subtyping method, 40 indicating the novel immune subtypes are an independent or different subtyping mechanism from the previous study(Figure 6A). These results suggested subtype-1 may have a better immune response than subtype-2, resulting in its better prognosis and well ICB therapy. This can make promotion in carrying out personalized therapy and improve the patients’ survival rate and may inspire other researchers or clinicians for developing a better usage of immunotherapy, which may have clinical significance.
However, we must admit that our research has some shortcomings. First, our research is based on bioinformatics data analysis, which is theoretically feasible but not experimentally verified. 41 Experimental verification is necessary if a certain result is to be obtained. Secondly, the samples of this study are mainly from western populations, so it is unknown whether the conclusions of this study could apply to eastern populations, which may be biased. Therefore, an eastern population research cohort is needed to complete the study. Whereas, two different immune-based subtypes were identified in this study with distinct mechanisms of immune dysfunction associated with different patient outcomes. This immune classification might serve as a basis for the design of clinical trials aiming to test biology-guided immunotherapy options in HNSCC. 39
Conclusions
In conclusion, we identified two prognostically and clinically relevant immune subtypes of HNSCC. The two HNSCC immune subtypes differ in multiple aspects including prognosis, clinical outcomes, immune characteristics, mutation pattern and biological functions. Results also suggest that patients in the two subtypes may have different ICB therapy responses, which needs to be verified in future experiments. These findings of the immune microenvironment may shed new light on the strategy of immunotherapy in HNSCC.
Supplemental Material
sj-txt-1-tct-10.1177_15330338211045823 - Supplemental material for Identification of Immune Subtypes for Predicting the Prognosis of Patients in Head and Neck Squamous Cell Carcinoma
Supplemental material, sj-txt-1-tct-10.1177_15330338211045823 for Identification of Immune Subtypes for Predicting the Prognosis of Patients in Head and Neck Squamous Cell Carcinoma by Jing Sun, Guiqing Fang, Zhibin Zuo, Xijiao Yu, Lande Xue, Chong Li and Shu Li in Technology in Cancer Research & Treatment
Supplemental Material
sj-txt-2-tct-10.1177_15330338211045823 - Supplemental material for Identification of Immune Subtypes for Predicting the Prognosis of Patients in Head and Neck Squamous Cell Carcinoma
Supplemental material, sj-txt-2-tct-10.1177_15330338211045823 for Identification of Immune Subtypes for Predicting the Prognosis of Patients in Head and Neck Squamous Cell Carcinoma by Jing Sun, Guiqing Fang, Zhibin Zuo, Xijiao Yu, Lande Xue, Chong Li and Shu Li in Technology in Cancer Research & Treatment
Footnotes
Acknowledgments:
We would like to thank Zhang Xiaoqi from West China Stomatology Hospital for his help with our bioinformatic analysis.
Author Contributions:
Conceptualization, Jing Sun and Chong Li; Formal analysis, Guiqing Fang; Methodology, Jing Sun and Guiqing Fang; Project administration, Shu Li and Chong Li; Supervision, Shu Li and Chong Li; Visualization, Zhibin Zuo, Xijiao Yu and Lande Xue; Writing – original draft, Jing Sun; Writing – review & editing, Guiqing Fang and Chong Li.
Conflicts of Interest
The authors declare no conflict of interest.
Ethical Approval
Not applicable for this article.
Funding
The study was supported by grants from the Open Foundation of Shandong Province Key Laboratory of Oral Tissue Regeneration (SDDX202003)to Jing Sun and the second batch of Science and Technology Plan Projects of Jinan Municipal Health Commission (2020-3-49)to Jing Sun; The Dean's Reserch Fund of Jinan Stomatological Hospital(2021-01).
Supplemental material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
