The Plaque Microbiota Community of Giant Panda (Ailuropoda melanoleuca) Cubs With Dental Caries

Introduction

Dental caries, as irreversible tooth decay, may cause pulpitis, severe toothache, chronic infections and eventually lead to health problems such as anorexia, difficulty foraging and malnutrition, as well as gastrointestinal, cardiovascular, renal and other systemic disorders (Selwitz et al., 2007; Maruyama et al., 2014). After years of research, the etiology and related risk factors of human dental caries have been studied in detail and explained from the perspective of oral microbes (Fleming, 1965; Tanzer et al., 2001; Munson et al., 2004). Nevertheless, only a few studies have focused on dental disease in animals. Except those which have focused on pathogenic animal models (Loyola and Leyva, 1990), a small group have focused on dental caries of companion and domestic animals (Sturgeon et al., 2012; Sturgeon et al., 2014). There have been very few studies of caries in wild animal species (Cuozzo and Sauther, 2006; Li et al., 2013; Malmsten et al., 2015), especially in giant pandas (Wang, 1961; Jin et al., 2012).

With a unique digestive system and foraging patterns, giant pandas spend more than 10 hours per day eating hard bamboo, as a result complete and healthy teeth are especially important for this rare species of bear (Hull et al., 2011). Dental caries is one of the most harmful diseases affecting the quality and health of both wild and captive pandas (Jin et al., 2012; Jin et al., 2015). The prevalence of dental caries of ancient pandas was about 25% (Wang, 1961). The rate of caries of wild pandas living in Qinling mountain, Shaanxi Province was 22% (Shuangyun et al., 2009) and captive born individuals at Beijing Zoo reached 62.5% (Jin et al., 2012). This rate is much higher than the prevalence of caries reported for black bears (Ursus americanus) (11%) living in Northern Wisconsin (Manville, 1978), wild brown bears (Ursus arctos) (11.2%) and Alaskan grizzlies bears (Ursus arctos horribilis) (14.7%) living in Switzerland (Wenker et al., 1999). Dental caries has become an urgent problem for captive pandas, especially amongst cubs.

Breast milk is the only food source for wild giant pandas until they are weaned at 1.5 years of age, after which their diet will consist of more than 99% bamboo (Guo et al., 2018). On the contrary, formula milk will be suppled to ensure captive giant panda cubs have enough nutrition for growth before being weaned. Artificial formula is very different from breast milk in terms of composition and nutrition, such as the species and percentage of carbohydrates (Xuanzhen et al., 2005). Due to the high correlation between the development of dental plaque and food composition, the various species and percentage of compounds in the food will create the specific dental plaque. Compared with wild individuals, captive giant pandas are more likely to have caries due to higher proportions of concentrate as part of their daily food intake and fruit, which may increase the intake of carbohydrates (Stromquist et al., 2009; Jin et al., 2012). The diet of panda cubs is particularly monotonous, except for breast milk. According to the ecological plaque hypothesis, dental plaque biofilm becomes pathogenic when external challenges drive it towards a state with a high proportion of acid-producing bacteria (Xie et al., 2010; Head et al., 2014). With the rapid development of biotechnology and bioinformatics methods in recent years, high-throughput sequencing has been proven useful for a better understanding of both cultured and un-cultured bacteria and the microbiota-related diseases, which has opened the way toward a new understanding of the composition and structure of an animal’s microbiomes and the factors that affect them (Teng et al., 2015).

Currently, there is no study investigating the difference of the oral bacterial community between giant panda cubs with healthy teeth and those with caries. Therefore, based on 16S rRNA gene sequencing, the present research is aimed to increase the understanding of giant panda caries by comparing the dental plaque microbiota between healthy and caries-active giant panda cubs.

Materials and Methods

Animal Management and Welfare

All of the giant panda cubs in this study were socially housed in outdoor enclosures during the day and returned to an indoors den to stay with their mother at night. The procedures for sample collection and handling were approved by the Institutional Animal Care and Use Committee of the Chengdu Research Base of Giant Panda Breeding (IACUC No. 2018).

Samples

In this study, eight captive panda cubs born in 2017 at the Chengdu Research Base of Giant Panda Breeding were selected. Based on the veterinarian’s clinical diagnosis of dental health, the eight cubs were divided into two groups when they reached 1 year old: Caries-active group (CA, n = 4) and Caries-free group (CF, n = 4). In order to ensure that the sequencing results can truly reflect the composition of dental plaque microbiota, sampling was performed 4 h after eating in the afternoon between 17:00 – 18:00. All individuals had no record of antibiotic use in nearly three months before sampling. All sample collections were performed without anesthetic or restraint benefiting from the trust of the cubs with their keepers. Briefly, sterile swabs were used to collect dental plaque sample from all the tooth surfaces from giant panda avoiding contamination from gingival bleeding. Each swab with a plaque sample was immediately released into a sterile, labeled EP tube, and then immediately placed in a transfer box with dry ice. All samples were transported on dry ice to our laboratory for DNA extraction.

DNA Extraction and 16S rRNA Gene Sequencing

DNA was isolated from the sample swabs using the procedure described by Verdon et al. (Verdon et al., 2014) with the Tiangen TIANamp Swab DNA Kit (Tiangen, China) according to the manufacture’s instruction for isolation of DNA. The DNA quality of each sample was determined by 1% agarose gel electrophoresis, and sent to Shanghai Majorbio Bio-pharm Technology Co., Ltd (Shanghai, China) for further concentration assessment, amplified and sequenced according to their standard procedures (Wang et al., 2016; Yin et al., 2016). The V3–V4 regions of the 16S rRNA gene were amplified using the primers 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACGCGGGTATCTAAT-3’) that targeted conserved sequences found in bacteria. PCR reactions were performed in triplicate 20 μL mixture containing 4 μL of 5× FastPfu Bufer, 2μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase and 10 ng of template DNA. The amplicons were then extracted from 2% agarose gels and further purified by using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified by QuantiFluor-ST (Promega, USA) according to the manufacture’s protocols. Purified amplicons were pooled in equimolar and paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Illumina, San Diego, USA) according to the instruction. The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: SRP257211).

Statistical Analysis

Raw reads were demultiplex and quality-filtered using QIIME (version 1.9.1). The taxonomy of each 16S rRNA gene sequence was analyzed against the SILVA 128/16s bacteria database with a confidence threshold of 70%. All operational taxonomic units (OTUs) at 97% identity were obtained using UPARSE (version 7.0). The Venn diagram and bar graph of each group was displayed by R software (v3.1.1). The Alpha diversity index was calculated by Mother (v.1.31.2), and A ranked abundance curve was drawn by R software (v3.1.1) to explain both the richness and evenness of species, and the Wilcoxon rank sum test were performed on Alpha diversity index between groups to obtain the difference of species diversity between two groups. Unweighted pair-group method with arithmetic means (UPGMA) cluster analysis with Bary-curtis distance matrix was performed by R software (V3.5.1). Beta diversity was compared with the representative sequences of OTUs for each group using the Bray-Curtis similarity coefficient by Vegan package of R (V3.5.1) and the difference between two groups was performed using Adonis. After 999 iterations, the Principal Co-ordinates Analysis (PCoA) display diagram were obtained, which were used to study the similarity or dissimilarity of sample community composition between samples. R (v3.5.1) was used to calculate the significance of the difference test between the two groups for the bacteria in the top 10 relative abundance. LEfSe software1 was used to make a LEfSe cluster diagram and an LEfSe linear discriminant analysis (LDA) diagram. Partial Least Squares Discriminant Analysis (PLS - DA) model analysis was performed by R software (V3.5.1). Receiver Operating Characteristic Curve (ROC) model analysis was performed by Graphpad Prism 9. A Wilcoxon test was used to look for differences between groups by R (v3.4.1). Microbial functions were predicted using PICRUSt (version 1.0.0) and aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. In this study, differences were considered significant when P < 0.05 and extremely significant when P < 0.01.

Results

Sequencing Date

After the raw data was filtered and spliced, a total of 350,420 high-quality tags were produced in the CA and CF groups, with an average of 43,802 sequences per sample (ranging from 41,246 to 46,372, Table S1). The average sequence length was 426 bp, with the maximum length being 486 bp and the shortest length being 209 bp. The total numbers of OTUs obtained was 337, among which 218 OTUs were shared by both groups, 82 and 37 OTUs were uniquely identified in the CA and the CF groups, respectively (Figure 1A).

Figure 1 (A) OUT distribution in two groups. The red and blue circles represented samples of dental plaque in the CA and the CF group, respectively. The overlap denoted the OTUs shared by two groups. (B) Rank-Abundacnce curve. Curves of different colors show the relative percentage of the number of species (Y-axis) of each sample rank of the number of species (X-axis). The position of the horizontal coordinate at the end of the extension of the sample curve is the number of species in the sample. When the curve declines smoother indicates that the species diversity of the sample is higher, while the rapid and steep decline of the curve indicates that the proportion of dominant colonies in the sample is high and the diversity is lower.

Microbial Diversity Analysis

The Chao1 and ACE richness index was higher in CA group than in CF group, but there was no significant difference between groups by Wilcoxon test (Table 1). However, the Shannon diversity index was significantly different between groups CA and CF (3.03 vs. 2.30, p = 0.03), the Simpson diversity index was significantly different between groups CA and CF (0.09 vs. 0.18, p = 0.03), demonstrating the higher bacterial diversity of dental plaques with caries compared to caries-free (Table 1). Good’s coverage estimator for each group was over 95%, indicating that the current sequencing depth was sufficient to saturate the bacterial diversity of dental plaques. In addition, Simpson index indicated that the bacterial-community distribution in the plaque samples was very uneven, which was also observed in the rank-abundance curve which had a high slope and a long tail comprised of relative low abundance OTUs (Figure 1B). The 100 most abundant OTUs represented 99.14% of all sequences. Most of the remaining OTUs were present at relative low abundance.

Table 1 Alpha diversity indices for dental plaque bacteria in each group at 97% identity.

Community-Composition Analysis

The bacterial distribution was characterized in terms of the relative taxonomic abundances. A total of 13 phyla, 24 classes, 51 orders, 98 families, 189 genera and 268 species were detected in the plaque samples. When comparing at the phylum level, the CA and CF groups shared the same top five bacterial species with little difference in overall relative abundance. The most dominant phylum of samples in the CA group was Proteobacteria (70.97%), followed by Bacteroidetes (7.49%), Firmicutes (7.18%), Fusobacteria (7.62%), and Actinobacteria (6.48%; Figure 2A). Proteobacteria (72.11%) was also the most abundant phylum in the CF samples, followed by Bacteroidetes (16.10%), and Firmicutes (8.41%), Actinobacteria (2.72%), Fusobacteria (0.41%; Figure 2A). The estimated cumulative relative abundance of these five dominant phyla was above 99% of the identified OTUs. The dental plaque bacterial community of giant panda cubs in the CA and CF groups were similar in composition. However, there was the discrepancy in their relative abundance at phyla level. Further analysis of the relative abundance showed that the relative abundance of Bacteroidetes in the CF group was significantly higher than that in the CA group (p = 0.03); whereas, the relative abundance of Fusobacteria in the CA group was significantly higher than that in the CF group (p = 0.03; Figure 3A).

Figure 2 Microbial community bar plot. (A) Microbial community bar plot at the phylum level. Bar charts showing the relative abundance of all phyla detected in the dental plaque collected from the giant panda cubs in the CA (CA1-CA4) and the CF group (CF1-CF4). The identities of the microbiome were shown with color blocks on the right. (B) Microbial community bar plot at the genus level. Bar charts showing the relative abundance of all genera detected in the dental plaque collected from the giant panda cubs in the CA (CA1 – CA4) and the CF group (CF – CF4). The identities of the microbiome were shown with color blocks on the right.

Figure 3 (A) Comparison of bacterial differences among diffierent groups of samples at the phylum level in top 5. (B) Comparative analysis of different species of bacteria among groups at genus level in top 10 (0.01<p<= 0.05, marked with "*"; if p>0.05, not marked).

At the genus level, a total of 189 bacterial genera were identified in our study, 171 in the CA group and 143 in the CF group. Dental plaque samples in CA group had the highest unclassified Neisseriaceae (16.61%) content, followed by Lautropia (14.48%), Neisseria (13.61%), Actinobacillus (10.38%) and unclassified Pasteurellaceae (7.81%). While dental plaque samples in the CF group displayed the highest amount of unclassified Neisseriaceae, accounting for 32.28% of the total microbial relative abundance, the second most abundant was Actinobacillus (24.95%), followed by Porphyromonas (10.67%), Bergeyella (4.41%) and Lautropia (4.40%; Figure 2B). However, the plaque bacterial communities in the CF group of giant panda cubs were notably different compared to the CA groups and had fewer number of species. In accordance with analysis of the relative abundance at genus level, the top 10 species were selected to evaluate the significance of the difference test between the two groups. The results showed that there were differences in unclassified Neisseriaceae, Actinobacillus and unclassified Pasteurellaceae between these two groups. The relative abundance of unclassified Neisseriaceae and Actinobacillus in the CF group was significantly higher than that in the CA group (p = 0.03); unclassified Pasteurellaceae in the CA group was significantly higher than that in the CF group (p = 0.03; Figure 3B). Our results showed that the dental plaque bacterial community in caries-active group is different with those in caries-free group at the genus level.

Bacterial Community Structure

To gain insight into similarities in the bacterial community structures between the two study groups, PCoA of beta diversity analysis was performed based on the Bray-curtis distances, which demonstrated different community structures among CA and CF groups. No sample’s microbiota overlapped with others and, a clear segregation in community structures was exhibited between the two groups, with the first two principal components representing 57.53% and 15.59% of the total variations (Figure 4A). Adonis testing confirmed that a significant separation occurred between the two groups (R² = 0.484, p = 0.027). Unweighted pair-group method with arithmetic means (UPGMA) cluster analysis also revealed that the samples formed well-separated clusters corresponding to the two groups, indicating a difference in the bacterial community structures (Figure 4B).

Figure 4 (A) PCoA based on the Bray-curtis distance matrix to find principal coordinates. (Description) X-axis, PCoA axis1 and Y axis, PCoA axis2. The scale of the X-axis and theY-axis are the projection coordinates of the sample points in the two-dimensional plane, respectively. A dot represents each sample, and different colors represent different groups. (B) Samples Clustering result (Description, Bray-curtis). The same color represents the samples in the same group. Short distance between samples represents high similarity.

Core Microbiome

To examine the existence of an identifiable common core microbiome, we defined a core as the group of members shared among the microbial community and represented the core by overlapping areas in the circles in a Venn diagram, at 97% identity. We identified 218 OTUs which were shared between two groups, occupying 64.69% of all OTUs (337 OTUs) and 99.26% of all OUT relative abundance (Figure 1A). Furthermore, 11 shared phyla (84.62%), 20 shared classes (83.33%), 38 shared orders (74.51%), 72 shared families (73.47%), and 125 shared genera (66.14%). These shared taxonomic members can be regarded as the core microbiome of dental plaques.

Among the 218 core OTUs, the 10 most relative abundant OTUs were OTU 263 (19.36%), OTU 205 (11.20%), OTU 130 (7.36%), OTU 195 (4.50%), OTU 26 (3.45%), OTU 208 (3.11%), OTU 131 (2.56%), OTU 283 (2.52%), OTU 87 (2.43%) and OTU 209 (2.12%) (Figure 5). The top 10 relative abundance OTUs in core microbiome accounted for 58.61% of all sequences. It is worth mentioning that among the top 10 relative abundance OTUs in core microbiome, except OTU 209 belongs to Phylum Firmicutes and OTU 195 belongs to Phylum Bacteroidetes, other 8 OTUs all belong to Phylum Proteobacteria. In the Venn diagram, the unique OTUs in each group were also observed with 82 and 37 unique OTUs found in the CA and CF groups, respectively. These unique OTUs were in relative low abundance, only containing 1–378 sequences, and can be considered as a variable microbiome.

Figure 5 The barplot of top 10 relative abudance OTUs in core microbiota. Red color represent Caries-active group, blue color represent Caries-free group.

Potential biomarkers

To identify potential biomarkers, LEfSe analysis, PLS – DA analysis and ROC curve model was performed. The cladogram of the LEfSe analysis was showed in Figure 6B. The red and blue parts represented the groups of the CA and CF, respectively. The red and blue nodes in the cladogram denoted the bacteria playing a critical role in the CA and CF group, respectively. The yellow nodes corresponded to the bacteria that did not have an important role in each group. The LDA results showed that there was a significant difference between the genera in the dental plaque samples from the CA group and the CF group (Figure 6A). The genera found in the CA group belonged to unclassified Pasteurellaceae, Fusobacterium, Paenibacillus and Leptotrichia while the genera found in the CF group belonged to Actinobacillus, unclassified Neisseriaceae and Porphyromonas. Since LEfSe was a strict tool, we also generated a PLS-DA model to identify potential biomarkers on genus level more strictly. The key genera with VIP > 1 and p < 0.05 were considered to be important contributors to the model. As shown in Figure 6C, a clear segregation in community structures was exhibited by PLS – DA model between the two groups. a total of 5 genera with a VIP score > 1 and p < 0.05 were identified as key genera responsible for significant differences in the community composition. Among them, 3 genera were significantly enriched in CA group, including unclassified Pasteurellaceae (VIP = 3.644, p = 0.03), Neisseria (VIP = 4.815, p = 0.03) and Corynebacterium (VIP = 1.968, p = 0.04). The other 2 genera were significantly more abundant in CF group, including unclassified Neisseriaceae (VIP = 5.727, p = 0.02) and Actinobacillus (VIP = 5.483, p = 0.03). A further ROC analysis model showed that unclassified Pasteurellaceae, unclassified Neisseriaceae and Actinobacillus had optimized diagnostic and classification capability (AUC = 0.999, p = 0.021), followed by Neisseria and Corynebacterium (AUC = 0.875, p = 0.083, Figure 6D). Therefore, unclassified Pasteurellaceae, unclassified Neisseriaceae and Actinobacillus were considered to be the most valuable potential biomarkers.

Figure 6 Identification of potential biomarkers using the LEfSe test, PLS-DA and ROC curve. (A) Linear discriminant analysis (LDA) demonstrated distinct microorganism enriched in the CA group and the CF group. The graph shows the LDA scores obtained from linear regression analysis of the significant microorganism groups in the two groups. When the default LDA value is more than 4.0 and the p value is less than 0.05, the result corresponds to a differential species. (B) Linear discriminant effect size cluster tree for 16S rRNA gene sequencing analysis. Different colors indicate different groups. Colored notes represent a group, color shading over the notes indicate the significant microbe biomarker in the group causing a significant difference in abundance, and the biomarker name is listed in the upper right corner. The yellow notes represent the biomarkers that do not show significant differences in abundance in the groups. (red indicating CA, blue indicating CF and yellow indicating non-significant). (C) Partial Least Squares Discrimination Analysis (PLS-DA) result. (Description) X-axis, X - Variate 1 and Y axis, X- Variate 2. The scale of the X-axis and theY-axis are the projection coordinates of the sample points in the two-dimensional plane, respectively. A dot represents each sample, red circle represents CA group and blue triangle represents CF group. (D) Receiver Operating Characteristic Curve (ROC) for screening for potential biomarkers. (Description) X-axis, 1 - Specificity and Y axis, Sensitivity. The red dotted line represents the random selection probability, Area Under Curve (AUC) represents the distinction of different categories.

Predictive Function Gene Analysis

The functional genes of the samples were analyzed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. All the levels of energy related metabolic pathways on KEGG level 2, for example: amino acid metabolism, carbohydrate metabolism, energy metabolism, glycan biosynthesis and metabolism, lipid metabolism, nucleotide metabolism, xenobiotics biodegradation and metabolism in the CA group were significantly higher than those in the CF group (all the p value = 0.03, Figure 7). We also investigated bacterial functions associated with human diseases which were significantly enriched in the CA group, such as renal cell carcinoma, pathways in cancer, tuberculosis, African trypanosomiasis, chagas disease (American trypanosomiasis), Epithelial cell signaling in Helicobacter pylori infection and Amyotrophic lateral sclerosis (ALS) (Table S2).

Figure 7 Mann-Whitney U test result of metabolism pathway difference (KEGG level 2). The red box represent CA group and blue box represent CF group. The value of p obtained by Mann-Whitney U test, all the metabolism pathways are significantly different between the two groups (all the p = 0.03).

Discussion

A comprehensive and thorough investigation of the bacterial diversity of plaque microbiota is essential for understanding their etiologies and for developing effective prevention and treatment strategies of dental caries. The introduction of high-throughput 16S rRNA sequencing has provided new insights into the compositions and structures of microbial communities (Yang et al., 2012). At present, giant panda oral microbes still need to be explored, in order to better classify and identify microorganisms, we have selected the more reported V3 and V4 variable regions which yield the highest accuracies of phylogenetic assignment and Illumina mature system in our study (Xu et al., 2018; Lee et al., 2021).

Different from human plaque research, collection of plaque samples from animals (especially wild animals) is particularly difficult due to their food intake and behavior. In order to get accurate and effective plaque bacterial structure information without using anesthesia, there needs to be a mutual trust relationship between the keepers and giant panda cubs. The use of anesthesia was abandoned for two reasons: 1) sampling after anesthesia may have negative health effects for the animal and reduce trust between animals and keepers, making subsequent research sampling more difficult; 2) the composition of respiratory anesthesia may affect the bacterial structure of oral plaque (Sturgeon et al., 2013). For these reasons, without anesthesia, we obtained 350,420 high-quality sequences with an average of 43,802 sequences per sample, which was much higher than reported in previous human studies and animal studies (Keijser et al., 2008; Ling et al., 2010; Sturgeon et al., 2014). According to Hamady and Knight’s (2009) research, a depth of coverage of approximately 1000 sequences per sample was sufficient to enable detection of species at 1% relative abundance with reasonable accuracy. Thus, the sequencing depth in our study was reasonable and large enough to enable detection of the vast majority of bacterial species in plaque samples. However, there are still some bacteria in the results that cannot be accurately classified on genus level, speculating that aspects of the bacterial community of giant pandas and other wild animals is still unknown.

It is worth mentioning that in our study the bacterial diversity of plaque in giant panda cubs were significantly increased in the group with dental caries (Shannon: CA VS. CF, p = 0.03; Simpson: CA VS. CF, p = 0.03). In most of the previous studies of dental caries, the microbial diversity of individuals with caries is lower than that of healthy individuals, and the diversity index decreases with severity of dental caries, which supports the notion of the “ecological catastrophe” (Li et al., 2005; Marsh, 2006; Preza et al., 2008). In some reports, however, we found similar results to our study, and therefore we presume that due to the altered oral environment, the relative abundance of beneficial bacteria that support oral health is reduced and core microflora function is disturbed, leading to the colonization of many opportunistic bacteria (Luo et al., 2012). This implies that changes of oral bacterial community structure may be related to the shift from healthy to diseased.

Microbiota can reach an equilibrium by competing for niche and nutrition and mediate colonization resistance against pathogenic bacteria (Chassaing and Cascales, 2018). A suitable relative abundance of different plaque microorganisms contributes to the health of the oral environment, once established, the resident oral microbiota is fairly stable and resilient to exogenous encounters (Marsh and Zaura, 2017). This stability can, however, be lost when environmental stress factors overrule resilience which will cause a colonization of rare bacteria. Our study with PCoA analysis showed that the caries-active samples were more dispersed compare with the caries-free samples, which can be interpreted as an imbalance of the bacterial community.

The core microbiome is considered to be a key component of the basic functions of all organisms, and is enriched, selected and inherited through evolutionary processes (Lundberg et al., 2012; Lemanceau et al., 2017). The OTUs shared between the eight giant panda cubs were determined in order to establish the core microbiome, which was common among the eight samples. In our study, 337 OTUs corresponding to 189 unique genus-level identifications were present in all samples. Here we only reported the 10 bacteria on genus level whose proportion exceeds 1% of all sequences, representing the most core microbiome. The CA group has 45 more unique OTUs than the CI group, which is in agreement with the findings of Schloss et al. (2011), that some extremely low relative abundance of rare species could be detected and might be related to dental caries. In the future, we will increase the sample size and usage of samples from twin cubs and further strengthen the research on core microorganisms.

At the phylum level, 13 different phyla were identified in the plaque samples of both CA and CF groups. Five phylum had relative abundances higher than 1% (Figure 2A), of which the phylum Proteobacteria, Bacteroidetes, Firmicutes, Actinobacteria and Fusobacteria have been previously detected in oral microbiota studies of both humans and animals (Jiang et al., 2013; Willis and Gabaldon, 2020). From the phyla relative abundance, the structure of giant panda cubs’ oral microbiota is more similar to that of cats (Sturgeon et al., 2014) and cattle (Borsanelli et al., 2018), followed by kangaroo infants (Chhour et al., 2010) swine (Kernaghan et al., 2012), and then humans, chimpanzees and bonobos (Li et al., 2013), with the most differences with canines (Sturgeon et al., 2013) and marine mammals (Bik et al., 2016). Bacteroidetes was significantly enriched in the CF group and Fusobacteria was significantly enriched in the CA group. By LEfSe test, Fusobacteria and Fusobacterium were siginificantly enriched in caries-active cubs. Fusobacterium has been reported as one of the most effective co-aggregators of human dental plaque, and was enriched in the CA group (Kolenbrander et al., 1989). Fusobacterium could co-aggregated to early colonizers, and express very high autoinducer 2 activity. Supported by the increased sugar metabolism that resulted in increased Glycan biosynthesis in the CA group, may have formed a strong matrix of biofilm which helped bridging bacteria, such as Fusobacteria and Corneybacerium, link the early and late colonizers in the biofilm, ultimately causing caries and caries recurrence (Kalpana et al., 2020).

At the genus level, compared with the CF group, the CA group mainly changed from the replacement of the dominant bacteria. The relative abundance of unclassified Neisseriaceae and Actinobacillus significantly decreased in the CA group even though they still had a high percentage. The remaining spaces were occupied by Lautropia, Neisseria, unclassified Pasteurellaceae, Moraxella and Canocytophaga. In the present study, using an LEfSe test, unclassified Pasteurellaceae, Fusobacterium, Paenibacillus and Leptotrichia were siginificantly enriched in giant panda cubs with dental caries. Species in the family Pasteurellaceae exhibit urease activity, can survive in an acidic environment and influence the development of dental caries (Morou-Bermudez et al., 2015; Nomura et al., 2020). Paenibacillus is used for the prevention of dental caries because of its effect of hydrolyzing the biofilm formed by Streptococcus mutans (Kolahi and Abrishami, 2013). We also found cases of oral disease caused by Paenibacillus (Al Qarni et al., 2017). We speculate that the enrichment of Paenibacillus in the CA group is due to the increased concentration of sugar substrates resulting from the increased intake of carbohydrate components in the diet. Leptotrichia has been proved to have high saccharolytic potential, is able to ferment a large variety of mono- and disaccharides to lactic acid, and is well adapted to thrive under conditions that are conducive to the formation of caries (Thompson and Pikis, 2012). Certain Leptotrichia species were shown to be negatively associated with elevated urease activity in plaques, which is correlated with dental health, and therefore positively associated with caries (Morou-Bermudez et al., 2015; Richards et al., 2017). According the results of PLS-DA and ROC model analysis, unclassified Pasteurellaceae, unclassified Neisseriaceae and Actinobacillus were considered as a potential biomarker for the diagnosis of dental caries in giant panda cubs. Next step, we will increase the number of samples and further study the potential biomarkers. Besides, dental caries are not caused by a specific pathogenic bacterium, and the intrinsic association between the different dental plaque microorganisms and caries in giant panda cubs needs further in-depth study.

In this study, zoonotic pathogens were identified in each oral sample, particularly Streptococcus. These data suggest that those zoonotic pathogens may be a common commensal in giant pandas. These findings highlight the zoonotic potential of the oral flora of healthy giant pandas and the inherent risks of exposure to saliva and oral exposure.

In this study we compared the differences in metabolic pathways between the CA group and the CF group at KEGG level 3, in addition to deliberately comparing the metabolic pathways between the two groups at KEGG level 2. Not surprisingly, all the metabolism pathways at KEGG level 2 were enriched in the CA group including the top three with the highest relative abundance of amino acid metabolism, carbohydrate metabolism and energy metabolism. As with most previously reported findings, enhanced metabolic pathways are positively associated with the development of caries, especially carbohydrate metabolism (Moye et al., 2014; Strużycka, 2014; Wang et al., 2019).

As the first research of giant panda dental plaque microflora with caries, this study explored the difference of plaque microbiota composition and function of giant panda cubs with and without dental caries, in which the characteristics and core bacteria of those were analyzed. In addition, the differences of dental plaque microflora could be used as a potential biomarker for the diagnosis of dental caries in giant panda cubs. These findings can further deepen the understanding of the panda cubs’ oral cavity microenvironment and provide valuable information for disease treatment and captive management.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by Institutional Animal Care and Use Committee of the Chengdu Research Base of Giant Panda Breeding.

Author Contributions

RM, RH, J-LG, D-WQ, and Q-GY contributed to the conception of the study. RM, J-LG, and X-YZ performed the experiment. RM, J-LG, D-WQ, and Q-GY contributed significantly to analysis and manuscript preparation. RM, D-WQ, and Q-GY performed the data analyses and wrote the manuscript. S-JC, X-BH, RW, Y-PW, QZ, S-YD, J-CL, and YB helped perform the analysis with constructive discussions. All authors contributed to the article and approved the submitted version.

Funding

Funding was provided and supported by the National Natural Science Foundation of China (U21A20193, No.31772484); Project of Sichuan Department of Science and Technology (No. 2017JQ0026) and the Chengdu Giant Panda Breeding Research Foundation (No. CPF2017-20).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the staff at the Chengdu Giant Panda Research Base of Breeding for assisting in the collection of samples. We also thank Dr Ayala James for the manuscript grammar checking work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2022.866410/full#supplementary-material

References