In this study, we examined a 23-year-old elderly female giant panda with a tumor located on the right side of her oral cavity. The tumor originated from the right buccal mucosa and measured 9.2 cm × 5.0 cm, as visualized on a computed tomography scan. Immunohistochemical results of the tumor were, AE1/AE3 (-), CK (-), EMA (-), S100 (+), SMA (-), CD34 (-), CD31 (-), Desmin (-), Vimentin (+ +), Calponin (-), CD99 (-), Ki67 (positive rate about 1%), and the immunophenotype and biological behavior were indicative of a low-grade malignant fibroma. Microbial sampling from various niches (mucosa, dental plaque, tumor surface, and vestibular sulcus) within the oral cavity were collected using sterile cotton swabs. Oral hygiene and surgical removal of the tumor were performed after sampling. We obtained a total of 16 samples, whose names and abbreviations are provided in Table 1 . The samples were immediately stored in sterile phosphate-buffered saline on dry ice and kept at −80°C for further analysis. The samples were collected in June 2022 at the Dujiangyan base of the China Conservation and Research Center for the Giant Panda in Chengdu, Sichuan province. All animal experiments were approved by Animal Care and Use Committee of the China Conservation and Research Center (Letter no: CCRCGP2020003).

Total genomic DNA samples were extracted using the OMEGA Soil DNA Kit (M5635-02; Omega Bio-Tek, Norcross, GA, USA), following the manufacturer’s instructions, and stored at −20°C prior to further analysis. The quantity and quality of the extracted DNA samples were measured using the NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively.

Polymerase chain reaction (PCR) of the full-length bacterial 16S rRNA genes was performed using the forward primer 27F (5’-AGAGTTTGATCMTGGCTCAG-3’) and the reverse primer 1492R (5’-ACCTTGTTACGACTT-3’). Thereafter, sample-specific 16-bp barcodes were incorporated into the primers for multiplex sequencing in the second step of PCR. Each PCR sample contained 5 μL of buffer (5×), 5 μL of GC buffer (5×), 0.25 μL of Q5 DNA polymerase (5 U/μL), 2 μL (2.5 mM) of dNTPs, 1 μL (10 μM) each of forward and reverse primers, 2 μL of DNA template, and 8.75 μL of ddH₂O. PCR was conducted under the following conditions: initial denaturation at 98°C for 2 min; followed by 25/10 cycles (for first and second PCR amplification steps, respectively) of denaturation at 98°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 90 s; and final extension at 72°C for 5 min. The PCR amplicons were purified using Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN) and quantified using the PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). Subsequently, the amplicons were pooled in equal amounts and subjected to single-molecule real-time sequencing using the PacBio Sequel platform at Shanghai Personal Biotechnology Co., Ltd (Shanghai, China). PacBio circular consensus sequencing (CCS) reads were generated from multiple alignments of sub-reads to reduce sequencing errors. In CCS, DNA polymerase reads a ligated circular DNA template multiple times to generate a consensus sequence from multiple reads of a single molecule. The raw sequences were filtered for a minimum of 3 passes through the PacBio SMRT Link portal (version 5.0.1.9585) until a minimum predicted accuracy of 99% was achieved. The predicted accuracy threshold was defined as the level below which CCS acted as noise. Files generated by the PacBio platform were trimmed to < 2,000 bp amplicon size.

Bioinformatics was conducted using QIIME2 2022.11 (Bolyen et al., 2019) with slight modification. Briefly, raw sequence data were subjected to demultiplexing using the demux plugin followed by primer cutting using the cutadapt plugin. Sequences were then merged, quality filtered, and dereplicated using the fastq-mergepairs, fastq-filter, and derep-fullength functions, respectively, in the Vsearch plugin. Non-singleton amplicon sequence variants (ASVs) were aligned with mafft (Katoh et al., 2002) and used to construct a phylogeny with fasttree2 (Price et al., 2009). Alpha- and beta-diversity metrics were estimated using the diversity plugin with samples rarefied to 989 sequences per sample. Taxonomy was assigned to ASVs using the classify-sklearn naïve Bayes taxonomy classifier in the feature-classifier plugin (Bokulich et al., 2018) against the SILVA Release 132 database (Koljalg et al., 2013). The raw reads were submitted to the NCBI Sequence Read Archive (SRA) database with the assigned Accession Number: PRJNA1053871.

Sequence data analyses were performed using the QIIME2 and R packages (v3.2.0). ASV-level alpha diversity indices, such as the Chao1 richness estimator (Chao, 1984), Observed species, Shannon diversity index (Shannon, 1948a; Shannon, 1948b), Simpson index (Simpson, 1949), Pielou’s evenness (Pielou, 1966), and Good’s coverage were calculated using the ASV table in QIIME2 and visualized as box plots. ASV-level ranked abundance curves were generated to compare the richness and evenness of the ASVs among samples. Beta diversity analysis was performed to investigate the structural variation of microbial communities across samples using Bray-Curtis metrics (Bray and Curtis, 1957) and visualized via principal coordinate analysis (PCoA) and unweighted pair-group method with arithmetic means (UPGMA) hierarchical clustering (Ramette, 2007). Significant variations in the microbiota structure among the groups were assessed by permutational multivariate analysis of variance (PERMANOVA) (McArdle and Anderson, 2001), analysis of similarities (Clarke, 1993; Warton et al., 2012), and permutation test of multivariate homogeneity of groups dispersions (Anderson et al., 2006) using QIIME2. The microbiota composition and abundance were visualized using MEGAN (Huson et al., 2011) and GraPhlAn (Asnicar et al., 2015). Venn diagram was generated using the R package “VennDiagram” to visualize the shared and unique ASVs among groups, based on their occurrence, regardless of their relative abundance (Zaura et al., 2009). Linear discriminant analysis effect size (LEfSe) was performed to detect differentially abundant taxa across groups using the default parameters (Segata et al., 2011). Microbial functions were predicted by phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt2) (Douglas et al., 2020) using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) database.

To conduct a comparative analysis of the oral microbiota, the samples obtained from the right and left regions were classified as tumor side and non-tumor side samples, respectively. The samples included in this study are enumerated in Table 1 . A total of 269,344 sequences were acquired (16,834/sample), which were merged, quality filtered, and dereplicated using the Vsearch plugin, to obtain 200,817 ASVs (12,551/sample). To achieve consistent sequencing depth across samples, the sequence data for each sample was extracted using the rarefaction method, resulting in 989 ASVs per sample. Rarefaction analysis was conducted to confirm the validity of the sequencing data. The rarefaction curves for all the samples demonstrated a plateau, indicating that the sequencing depth was adequate to capture the majority of gene diversity ( Supplementary Figure S1 ). The distribution of sequenced lengths ranged between 572–1899 bp, with the majority of sequences (n = 23,230) being 1463 bp ( Supplementary Figure S2 ).

Taxonomic richness varied at the different levels across all samples, resulting in the identification of 31 phyla, 55 classes, 89 orders, 136 families, 219 genera, and 402 species ( Figure 1A ). At the species level, the upper left molar and lower left molar samples exhibited the highest abundance, with 248 and 209 species, respectively. In contrast, the species content of the lower left buccal mucosa was the lowest, amounting to only 30 species. According to the abundance ranking curve ( Figure 1B ), the microbial composition of the left dental plaque appears to be relatively uniform, while those of the other regions appear to be uneven, as observed by the steep slope of the curves.

The top 20 abundant species identified in our study were Ottowia sp. canine oral taxon 014, Spodiobacter cordis, Porphyromonas sp. feline oral taxon 110, Moraxella sp. canine oral taxon 396, Neisseria shayeganii, Capnocytophaga sp. H4358, Fusobacterium simiae, Neisseria sp. feline oral taxon 145, Xanthomonadaceae bacterium feline oral taxon 091, Aggregatibacter aphrophilus, Lautropia sp. canine oral taxon 060, Globicatella sp. feline oral taxon 122, Glaesserella parasuis, Cardiobacterium sp. canine oral taxon 177, Neisseria sp. feline oral taxon 078, Leptotrichia sp. canine oral taxon 345, Neisseria sp. VA252/2008, Campylobacter sp. feline oral taxon 100, Streptococcus parasuis, and Streptococcus minor ( Supplementary Table S1 ).

Owing to the complex and diverse nature of the oral microbiota composition of the giant panda, we conducted a comprehensive detection of its microbial species. The results revealed the presence of Ottowia sp. canine oral taxon 014, S. cordis, Moraxella sp. canine oral taxon 396, N. shayeganii, Capnocytophaga sp. H4358, Xanthomonadaceae bacterium feline oral taxon 091, A. aphrophilus, Lautropia sp. canine oral taxon 060, and Cardiobacterium sp. canine oral taxon 177 in various niches of the oral cavity. However, the other species were confined to specific regions within the oral cavity, suggesting a complex and diverse microbial composition in these areas.

The alpha diversity indexes Chao1 and observed species represent community richness, Shannon and Simpson indices represent community diversity, Pielou’s evenness index represents community evenness, and Good’s coverage represents the coverage of species in the community. In the tumor side and non-tumor side samples, the mean Chao1 indices were 109.70 and154.91, observed species indices were 107.93 and 153.39, Shannon indices were 5.05 and 5.47, Simpson indices were 0.93 and 0.95, and Pielou’s evenness indices were 0.76 and 0.79, respectively. However, the Good’s coverage was 99% for both samples, suggesting that the sequencing depth was sufficient to indicate the diversity of the samples ( Figure 2 ). These results indicate no significant difference in the alpha diversity of the oral microbiota in the tumor side and non-tumor side of the giant panda.

The beta diversity analysis was conducted to assess the similarity of microbial structure among different samples using the Bray-Curtis distance algorithm and PCoA analysis. The results revealed a difference between the community structure of the tumor side and non-tumor side samples, with the first two principal coordinate components accounting for 32.9% and 20.7% of the total variation, respectively ( Figure 3A ). The PERMANOVA analysis confirmed a significant difference between the tumor side and non-tumor side samples (P < 0.05). Similarly, the UPGMA analysis also demonstrated a notable difference between the microbial communities of the tumor side and non-tumor side samples ( Figure 3B ). Furthermore, the microorganisms on the teeth (including canine teeth and molars) as well as oral mucosa exhibited prominent clustering. Altogether, these results indicate a significant difference in the microbial composition of the tumor side and non-tumor side samples.

The number of ASVs between the tumor side and non-tumor side samples was significantly different, with the tumor side containing 331 ASVs and non-tumor side containing 643 ASVs, among which 169 were common to both sides ( Figure 4A ). These results indicate that the microbial abundance is relatively low in the tumor side compared to the non-tumor side.

The top 20 most abundant species exhibited a distinct distribution pattern within the oral cavity of the giant panda. Specifically, the abundance of F. simiae (99.6%), Porphyromonas sp. feline oral taxon 110 (99.3%), Campylobacter sp. feline oral taxon 100 (98.9%), Neisseria sp. feline oral taxon 078 (98.3%), Xanthomonadaceae bacterium feline oral taxon 091 (86.4%), Cardiobacterium sp. canine oral taxon 177 (83.7%), and Leptotrichia sp. canine oral taxon 345 (82.0%) was significantly higher in the tumor side, with F. simiae, Porphyromonas sp. feline oral taxon 110, Campylobacter sp. feline oral taxon 100, and Neisseria sp. feline oral taxon 078 present exclusively in the tumor side of the oral cavity. In contrast, the abundance of G. parasuis (86.9%), Neisseria sp. VA252/2008 (86.1%), and S. parasuis (84.5%) was significantly higher in the non-tumor side of the oral cavity. The outer periphery of the tumor exhibited the highest abundance of Porphyromonas sp. feline oral taxon 110 and F. simiae ( Figure 4B ).

The cluster heatmap of the groups revealed a significant difference in microbial clustering between the tumor side and non-tumor side of the oral cavity ( Figure 4C ). The tumor side exhibited a more uniform microbiome compared to the non-tumor side, with F. simiae, Porphyromonas sp. feline oral taxon 110, Campylobacter sp. feline oral taxon 100, Neisseria sp. feline oral taxon 078, Xanthomonadaceae bacterium feline oral taxon 091, Cardiobacterium sp. canine oral taxon 177, and Leptotrichia sp. canine oral taxon 345 as the primary strains. Notably, F. simiae, Porphyromonas sp. feline oral taxon 110, Campylobacter sp. feline oral taxon 100, and Neisseria sp. feline oral taxon 078 were present exclusively in the tumor side. Furthermore, N. shayeganii was found to be present in low abundance in the mucosa surrounding the tumor, whereas it was significantly more abundant in the corresponding regions in the non-tumor side of the oral cavity (68.3%); however, the difference in the abundance of N. shayeganii in the dental plaques of tumor side (56.1%) and non-tumor side (43.9%) regions was not significant.

After investigating the variations in microbial community composition, we determined the specific species in the tumor side and non-tumor side that contribute to these differences, using LEfSe analysis. We identified robustly differential species (marker species; P < 0.05) between the tumor side and non-tumor side samples, based on the linear discriminant analysis scores (> 4). The predominant marker species in the tumor side were F. simiae, Porphyromonas sp. feline oral taxon 110, Xanthomonadaceae bacterium feline oral taxon 091, Moraxella sp. canine oral taxon 396, Campylobacter sp. feline oral taxon 100, and Neisseria sp. feline oral taxon 078, while the predominant species in the non-tumor side were Capnocytophaga sp. H4358, G. parasuis, and Globicatella sp. Feline oral taxon 122, with Streptococcus being the primary genus (marker taxon; Figure 5 ). These results suggest that the marker species around the oral tumor were F. simiae, Porphyromonas sp. feline oral taxon 110, Xanthomonadaceae bacterium feline oral taxon 091, Moraxella sp. canine oral taxon 396, Campylobacter sp. feline oral taxon 100, and Neisseria sp. feline oral taxon 078. Among them, F. simiae, Porphyromonas sp. feline oral taxon 110, Campylobacter sp. feline oral taxon 100, and Neisseria sp. feline oral taxon 078 are consistent with the results of the taxonomic analysis conducted on the tumor side.

KEGG metabolic pathway analysis revealed that the oral microbiota of the giant panda was predominantly associated with metabolism, especially cofactor and vitamin metabolism, amino acid metabolism, and carbohydrate metabolism, as well as genetic information processing, including replication and repair activities. The abundance of metabolic pathways was found to be higher on the tumor side compared to the non-tumor side; however, the overall metabolic pattern remained consistent on both sides. ( Figures 6A, B ). According to the KEGG orthologous group description of each metabolic pathway, we predicted that the metabolic disparity between the tumor side and non-tumor side may be the bacterial invasion of the epithelial cells (P < 0.001; Figure 6C ).