The recruited patients were newly diagnosed with nasopharyngeal carcinoma by the Department of Otolaryngology, The First Affiliated Hospital of Sun Yat-sen University from April 2021 to September 2021. The patients and their families were aware of the study methods and obtained written consent. The present study protocol was approved by the ethics board of The First Affiliated Hospital of Sun Yat-sen University. Finally, 64 NPC patients were recruit in this study.The inclusion criteria were as follows:

1) Aged from 20 to 80 years old

2) Cancer clinical stage I - IV

3) WHO pathological type II-III

The exclusion criteria were as follows:

1) A history of antibiotic use in the past 3 month.

2) A history of radiotherapy or chemotherapy.

3) Severe periodontitis or dental caries.

4) Already have oral mucositis or recurrent attacks

All patients received VMAT radical radiotherapy at follows. GTV (gross target volume): 68.1Gy/30Fr; CTV1 (Clinical target volume 1/high-risk area): 60Gy/30Fr; CTV2 (Clinical target volume 2/low risk area): 54Gy/30Fr. All patients also received whole course concurrent chemotherapy (Docetaxel 60mg/m²+Cis-platinum 60mg/m²). Six month later after the chemotherapy and radiotherapy were down, all patients go to department of stomatology for sample collection. Samples were obtained directly from the mucosal surface of the soft palate by throat swab, especially at the surface of injury mucosa. Before sampling, patient rinsed mouse with water without disinfection mouthwash. The contact time with air was as short as possible, and the samples were flushed into the EP tube with sterile normal saline. The tubes were immediately transferred to the ultra-low temperature refrigerator at -80°C in the laboratory.

The evaluation of scoring scales of the severity of oral mucositis is based on WHO scale and OMAS scale, followed by reference (Hong et al., 2019). WHO scale evaluates in a categorical scale of 0 to 4 objective signs and patientreported symptoms. Patients of 0 score represent Normal healing (N) group. Patients of 1 to 2 score represent M group, and patients of 3 to 4 score represent S group. The OMAS scale is based solely on objective signs of erythema and ulceration. OMAS scores reported here could range from 0 to 45 and represent the aggregated scores from nine intra-oral sites evaluated. For the power calculation, the balanced one-way analysis of variance power calculation was achieved using SPSS and R package “pwr”. The statistical power value = 81.92% for SPSS, and 88.28% for R package. Hence, the sample numbers in each group were sufficient.

Microbial DNA was extracted from oral samples using the HiPure Stool DNA Kits (Magen, Guangzhou, China), and the following protocol was based on the manufacturer’s instructions. PCR amplification of the V3-V4 region of the 16S rRNA gene was performed under conditions of denaturation. Primers 341F: CCTACGGGNGGCWGCAG and 806R: GGACTACHVGGGTATCTAAT were used to amplify the 16S rRNA gene fragment. For amplicon sequencing, libraries were prepared with MiSeq library preparation Reagent Kit v3 (New England Biolabs, USA) and then sequenced on a NovaSeq 6000 Sequencer platform (Illumina, USA).

For data filtering, FASTP (version 0.18.0) was used to filter raw data from the Illumina platform. For sequence splicing, FLASH software (version 1.2.11) was used to combine clean reads into tags according to the threshold of minimum overlap of 10bp and maximum mismatch rate of 2%. Noisy sequences of raw tags were filtered under specific filtering conditions (Edgar, 2013) to obtain the highquality clean tags. Filled tags were clustered into Operational taxonomic units (OTUs) according to ≥97% similarity using UPARSE(version 9.2.64) software (Bolyen et al., 2019). All chimeric tags were removed using UCHIME algorithm (Quast et al., 2013) and finally obtained effective tags for further analysis. The tag sequence with highest abundance was selected as representative sequence within each cluster.

Alpha diversity analysis was calculated in QIIME software(version 1.9.1). The exponential dilution curve of diversity and rank-abundance curve were drawn using GGPlot 2 package in R language. Beta diversity analysis was performed using Muscle (version 3.8.31) for OTU representative sequence alignment. Multivariate statistical techniques including PCoA (principal coordinates analysis) of (Un) weighted unifrac, jaccard and bray-curtis distances were generated in R project Vegan package (Edgar et al., 2011) (version 2.5.3) and plotted in R project ggplot2 package (Nilsson et al., 2018) (version 2.2.1). To indicate species analysis, Venn Diagram package in R language was used for Venn analysis, and the R language UpSet R package was used to analyze shared unique species/OTUs between groups. Linear discriminant analysis (LDA) effect size (LEfSe) analysis was performed to identify taxa showing the most significant differences in microbial abundance between groups (Segata et al., 2011). Only taxa with LDA scores > 2.0 and p < 0.05 are shown. KEGG (Kyoto Encyclopedia of Genes and Genomes) metabolic pathway analysis of sample bacteria or archaea was performed using PICRUSt (Parks et al., 2014).

Associations between the clinical characteristics were performed by Pearson’s Chi-square test or Fisher’s exact test. The p value less than 0.05 were considered as statistically significant (*, <0.05; **, <0.01; ***, <0.001). The Turkey-HSD analysis was performed to analyze the alpha-diversity differences. In LEfSe analysis, phenotypic prediction Bugbase was revealed by contribution stacking diagram and analyzed by Turkey-HSD test, P <0.05, with * for significant correlation. Microbiome phenotype analysis was used by BugBase, and tested by Pearson’s Chi-square test (*, <0.05; **, <0.01; ***, <0.001). All data were analyzed using GraphPad Prism v. 6.01, and R package v. 3.3.2 (Xia et al., 2018).

A total of 64 cases were recruited in this study. All samples were collected at the six months after radiotherapy and chemotherapy. The patients were divided into three sub-groups based on the healing degree of oral mucositis. Patients with severe delayed healing were divided into S group, and patients with mild delayed healing were divided into M group, and patients with normal healing were divided into N group. Typical injury presentations of erythema and ulceration were mostly on the non-keratinized mucosa, such as soft palate and pharyngeal area ( Figure 1A ). The severity of oral mucosal injury was based on WHO score scales (Hong et al., 2019) ( Figure 1B ). Another commonly used scoring scale called OMAS, was also used to evaluate the severity of oral mucosal injury among two groups ( Figure 1C ). There was no significant differences in the clinical characteristics, including age, gender, BMI, alcohol drinking, smocking, diabetes status among the three groups ( Table 1 ).

The oral microbiota was assessed by 16S rRNA gene MiSeq sequencing. The relative bacterial evenness was evaluated by rank abundance curves, which exhibited similar patterns in all samples ( Figure 2A ). Alpha-diversity was used to calculate and evaluate the differences in bacterial diversity among the three groups, which was analyzed on OTU level and revealed by Sob, Shannon and Simpson index. The results showed that the Alpha-diversity of S and M groups were higher than N group based on Sob index, but was no significant difference between S and M groups. Based on Shannon and Simpson index, three was no significant difference among three groups ( Figures 2B-D ). In addition, the Venn diagram revealed the 375 of the total 806 OTUs were shared among the three groups. Notably, 85 OTUs were unique for the S group ( Figure 2E ). To reveal the microbiome space between three samples, beta-diversity was calculated by the unweighted UniFrac method of the principal coordinate nanlysis (PCoA) based on OTU level. The results showed that the alpha-diversity and beta-diversity of oral microbial communities in these three groups were not significantly different ( Figure 2F ).

The average composition of oral bacterial was shown at phylum, family, genus, and species levels respectively. At the phylum level, Bacteroidetes, Firmicutes, Fusobacteria, Proteobacteria, and Actinobacteria were the top five dominant bacterial phyla in the three groups. Actinobacteria was increased in the S group than M group, and was also added in S group than N group ( Figure 3A ). At the family level, Prevotellaceae, Veillonellaceae, Fusobacteriaceae, Leptotrichiaceae, and Flavobacteriaceae were the top five bacterial family in three groups, and Veillonellaceae was significantly enriched in S group than M and N groups. In addition, Neisseriaceae was decreased in S and M groups than N group. Most importantly, Actinomycetaceae was significantly added in S group than M and N groups ( Figure 3B ). At the genus level, Prevotella, Fusobacterium, Leptotrichia, Capnocytophaga, and Veillonella were the top five genus in three groups, and Veillonella was significantly increased in S group than N group ( Figure 3C ). At the species level, Leptotrichia_hongkongensis, Neisseria_mucosa, Capnocytophaga_granulosa, Streptococcus_oralis_subsp_dentisani_clade_058, Corynebacterium_matruchotii were the top five species in three groups ( Figure 3D ). The Circos circle diagram also showed the tendency, reflecting the composition proportion and distribution of dominant microbial communities of S group was different from another two groups at family and genus level ( Figures 3E, F ).

To study the function and metabolism changes in oral microbial communities, the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) was used to analyze all OTUs. The PICRUSt analysis identified top ten KEGG pathways in significant differences among all three groups ( Figures 4A-C ). The pathways of metabolism of degradation of organic acids and amino acids were significantly increased in S group compared with those of M group and N group, such as benzoate and valine degradation. In contrast, the pathways involved in glycosaminoglycan degradation, biofilm formation, flagellar assembly, and metabolism of xenobiotics by cytochrome P450 were inhibited in S groups related to thosed in M and N groups. In addition, pathways of D-glutamine metabolism, fatty acid biosynthesis, and aminoacyl-tRNA biosynthesis were down-regulated in M group compared with N group.

To study the difference of phenotype of oral microbiota among three groups, we used the BugBase analysis to reveal the characteristics of the microbial phenotype of all groups, which were divided into nine phenotypic classifications. We found that in the aerobic phenotype, the relative abundance of bacteria in S groups was significantly higher than those in M group and N group. Similarly, the relative abundance of bacteria of S group was also higher than those in M and N groups in the phenotype of biofilm formation. In addition, we also found that the relative abundance of both S and M groups were lower than N group in gram-negative phenotype, while that of S and M groups were higher than N group in gram-positive phenotype. We did not found significant difference between S and N groups in the phenotypes of anaerobic, contains mobile elements, facultatively anaerobic, pathogenicity, and oxidative stress tolerant. These results showed the severe delayed healing group had higher bacteria abundance level of aerobic, biofilm formation, and gram-positive, but lower abundance level of gram-negative phenotype ( Figure 5 ).

LEfSe analysis was used to determine and distinguish the composition of oral microbiome among all three groups. The phylogenetic tree showed the phylogenetic relationships among the three groups, and a cladogram map was drawn to show the distribution of the dominant groups ( Figure 6A ). The different taxa were further extracted and shown on the bar chart, and only when the LDA score was greater than the preset value of 4, there were significant species differences among the three groups ( Figure 6B ). The length of the bars indicated the degree of impact of species that differed significantly among three groups. The results revealed 8 taxa of the oral microbiome, including Actinobacteria and Veillonellaceae, which were extremely enriched in the S group. For the M group, 9 kinds of microbial biomarkers, including Firmicutes, Negativicutes, Bacteroidia, Prevotella, and Selenomonas, were the predominant flora. For the N group, there were 19 taxa of bacterial markers were accumulated, including Neisseria, Flavobacteriia, and Saccharibacteria_TM7.

There was significant dysbiosis of the oral bacterial microbiome in S groups, as shown by differences in bacterial composition, diversity, and function among the three groups. So we focused on the selected 36 bacterial biomarkers above and clustered them based on the abundance profiles ( Figure 6C ). We identified the one pathogen oral bacterial co-abundance groups (pCAGs) (e.g. Actinobacteria and Veillonellaceae), and two kinds of biofilm oral bacterial co-abundance groups (bCAGs) (e.g. Neisseria, Bacteroidia, Prevotella, and Flavobacteriia), (e.g. Bacteroidia, Prevotella, and Selenomonas). It indicated that the oral bacterial microbiome in S groups had greater difference and less connection between M and N groups. The connection of oral bacterial microbiomes between M and N groups were more closed.

In addition, a stringent network analysis was performed to obtain further insights with respect to correlative bacterial populations at the all levels. It reveals that the apparent dissimilarity of the microbial composition when comparing the three groups. The resulting network showed significant co-occurrence and anti-occurrence ( Figure 6D ). The oral pathogen of 5 kinds of Actinobacteria co-occurred with 3 kinds of Veillonellaceae, and anti-occurred with other biofilm oral bacteria. In the bCAGs groups, nearly all the bacteria were co-occurred with others. We also found the Saccharibacteria_TM7 were independent and had limited connection with other bacteria.

To clarify if the Veillonella and Actinomyces can be predictable markers of severe delayed healing of oral mucositis, we established the ROC curves to study the predictive ability of Veillonella and Actinomyces. The results indicates these two diagnostic biomarkers is effective and has diagnostic value ( Figures 7A, B ).