Tongue coating microbiome composition reflects disease severity in patients with COVID-19 in Nanjing, China

ABSTRACT

Our purpose is to investigate the relationship between the microbiota of patients' tongue coating microbiota and the severity of COVID-19, and to identify the severity of COVID-19 patients' condition as early as possible. The participants were categorized into three groups: healthy controls (Con group) consisting of 37 individuals, patients with mild to moderate symptoms (M group) comprising 49 individuals, and patients with severe and critical symptoms (S-C group) consisting of 44 individuals. We collected oral swabs from all participants and performed 16S rRNA gene sequencing to analyze the microbiome. The α and β diversity differences were assessed respectively. Additionally, we employed the Linear Discriminant Analysis Effect Size (LEfSe) analysis to evaluate taxonomic differences among the three groups. Our findings revealed a significantly higher richness of tongue coating microbiota in both the S-C group and M group compared to the Con group. When compared with Con group, decreased Prevotella, Neisseria, Fusobacterium and Alloprevotella, and over-expressed Streptococcus and Rothia in M and S-C group were identified. LEfSe analysis indicated a greater abundance of Pseudomonas, Acinetbacter, Lactobacillus, Corynebacterium, Rothia in S-C group. Our study suggests a potential association between tongue coating microbiome and the severity of COVID-19 patients.

KEYWORDS:

• COVID-19
• oral microbiota
• tongue coating
• characterization
• severity

Introduction

COVID-19 is an acute respiratory disease caused by a novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Its clinical manifestations are mostly non-specific symptoms [1]. The most common symptoms are fever and respiratory symptoms. While most COVID-19 patients exhibit asymptomatic or mild symptoms, approximately 14% of patients demonstrate either normal or severe manifestations, characterized by evident signs of infiltration or respiratory damage. Additionally, around 5% of patients progress to a critically ill state [2]. Critical illness in COVID-19 can lead to severe pneumonia, pulmonary edema, acute respiratory distress syndrome, multi-organ failure, and ultimately, death. Importantly, individuals of any age without pre-existing health conditions can develop critical illness following COVID-19 infection. Therefore, it is imperative to ascertain specific biomarkers that can identify individuals at a high risk of progressing to severe illness.

Numerous studies have indicated the significant involvement of the microbial community in COVID-19 [3,4]. Alterations in the richness and diversity of intestinal and airway microbial populations have been observed in COVID-19 patients. Notably, the proportion of beneficial bacteria in the saliva of COVID-19 patients tends to decrease, while the prevalence of pathogenic bacteria shows an upward trend [5–7]. Cao et al. demonstrated a depletion of butyrate-producing bacterial communities and an increase in opportunistic pathogens in severe cases of COVID-19, in contrast to mild to moderate cases. The reduction in fecal butyrate levels was associated with elevated plasma levels of the pro-inflammatory cytokine IL-10 and the chemokine CXCL-10 [8]. To date, there have been no investigations examining the relationship between the tongue coating microbiota and the severity of COVID-19 patients.

In this study, our objective was to evaluate and compare the diversity and composition of the oral microbiota in the tongue coating of COVID-19 patients with mild symptoms, severe symptoms, and a control group. This analysis aims to shed further light on the potential role of the microbiota in the pathogenesis of COVID-19.

Materials and methods

Study participants

A total of 93 COVID-19 patients and 37 healthy control volunteers (referred to as the Control group) were recruited from Nanjing First Hospital and Nanjing Integrated Traditional Chinese and Western Medical Hospital between December 2022 and January 2023. Among the COVID-19 patients, 49 individuals presented with mild to moderate symptoms (referred to as the M group), while 44 patients exhibited severe and critical symptoms (referred to as the S-C group). These patients were confirmed positive for SARS-CoV-2 through quantitative reverse transcription PCR (RT-qPCR) analysis performed on nasopharyngeal swabs collected by hospital staff. The specific classifications were as follows: 1) mild, indicating the absence of radiographic evidence of pneumonia; 2) moderate, characterized by pneumonia accompanied by fever and respiratory tract symptoms; 3) severe, characterized by respiratory rates ≥ 30 breaths per minute, oxygen saturation ≤ 93% when breathing ambient air, or PaO2/FiO2 ≤ 300 mmHg; 4) critical, representing cases with respiratory failure requiring mechanical ventilation or organ failure necessitating intensive care. The study protocol was approved by the institutional review board of Nanjing Integrated Traditional Chinese and Western Medical Hospital, and all the experiments were performed in accordance with approved guidelines and regulations.

Samples collection and DNA extraction

To ensure the precision of tongue coating collection and minimize the influence of food residues, several measures were implemented. Participants were instructed to observe a fasting period, refrain from smoking and consuming alcohol, and avoid performing any oral hygiene procedures for a duration of 4 hours prior to sample collection. A sterile throat swab was employed, and it was gently rolled along the base of the tongue, ensuring approximately 30 brush strokes to collect the tongue coating. The swab was then immersed in a centrifuge tube containing phosphate buffer solution, while continuous stirring facilitated the transfer of microorganisms from the swab to the buffer. This process was repeated twice, and the collected samples were combined to ensure an adequate quantity of tongue coating samples. Following that, the collected tongue coating samples were subjected to centrifugation at 4°C and 5000 g for 5 minutes, with subsequent removal of the supernatant. The remaining material was stored at −80°C in a refrigerator. All items used were sterile. The DNA extraction process adhered to a previously established protocol [9].

PCR amplification

The bacterial 16S rRNA gene V3–4 hypervariable region was amplified using the primers 338F (ACTCCTACGGGGAGGCAGCAG) and 806 R (GGACTACHVGGGTWTCTAAT). To enable sample identification, a 10-digit barcode sequence was added to the 5’ end of both the forward and reverse primers. The PCR amplification was conducted using the Mastercycler Gradient (Eppendorf, Germany) with a reaction volume of 20 μl. The reaction mixture contained 4 μl of 5×FastPfu Buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl of the forward primer (5 μM), 0.8 μl of the reverse primer (5 µM), 0.4 μl of FastPfu Polymerase, 0.2 μl of BSA, 10ng of template DNA, and the final volume was adjusted with ddH2O to 20 μl. The cycling parameters for PCR amplification were as follows: an initial denaturation step at 95°C for 3 minutes, followed by 27 cycles of denaturation at 95°C for 30 seconds, annealing at 55°C for 30 seconds, extension at 72°C for 45 seconds, and a final extension step at 72°C for 10 minutes. To minimize PCR biases, three PCR products from each sample were pooled together. The PCR products were purified using the QIAquick Gel Extraction Kit (QIAGEN, Germany), quantified using RT-PCR, and subsequently subjected to sequencing.

Sequencing processing and analysis

The purified amplicons were combined in equimolar ratios and subjected to paired-end sequencing using an Illumina Novaseq 6000 PE250 platform (Illumina, San Diego, USA). The obtained raw data underwent filtering using QIIME (v 1.8.0), which involved the removal of dereplicated reads and those shorter than 150bp. The filtered reads were then processed using Unoise3 to generate amplicon sequence variants (ASVs) through denoising. To assign species classification information to each ASV, a comparison was made against the SILVA database (version 138). All raw reads have been deposited in the NCBI Sequence Read Archive (SRA) database under the accession number PRJNA970271.

Statistical analysis

To analyze continuous variables, we employed independent t-tests, White’s nonparametric t-test, and Mann-Whitney U tests. For categorical variables between groups, we used either Pearson’s chi-square test or Fisher’s exact test, depending on the validity of assumptions. The QIIME software was utilized to assess α diversity by calculating the Shannon index and Simpson index. To compare diversity differences among groups, we employed Partial Least Squares Discriminant Analysis (PLS-DA) for β diversity testing. To identify key microbes associated with different groups, we performed Linear Discriminant Analysis Effect Size (LEfSe) analysis, using an LDA threshold of 3. Additionally, we used Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) analysis to predict Kyoto Encyclopedia of Genes and Genomes (KEGG) biochemical pathways. Statistical analysis was conducted using SPSS V19.0 (SPSS Inc., Chicago, IL) and STAMP V2.1.3. Graphs were prepared using GraphPad Prism version 6.0 (San Diego, CA). A p-value <0.05 was considered statistically significant.

Results

Baseline characteristics of participants

A total of 37 healthy volunteers, 49 patients with mild to moderate disease, and 44 patients with severe and critical disease were included in this study. The demographic characteristics of all participants are presented in Table 1. It was observed that the S-C group had a higher average age compared to the Con group. Furthermore, significant differences were found in terms of age, hypertension, and coronary heart disease.

Table 1. Clinical characteristics of enrolled patients and healthy controls.

Characteristics	Control group(n=37)	M group (n=49)	S-C group (n=44)	p value
Age	57.73±12.56	74.45±13.26	81.40±13.67	.027*
Sex (male)	16	30	31	.063
Hypertension (Yes)	14	30	35	.045*
Diabetes (Yes)	8	10	13	.132
Coronary heart disease (Yes)	8	15	18	.036*

Microbial diversity and richness between three groups

After sequencing and quality filtering, more than 11 million tags and a total of 6209 ASVs were obtained with the dominant length of tags locating among 400-440bp (Figure 1a). As shown in the rarefaction curves in Figure 1b, the numbers of ASVs in the plot increased with the numbers of sequences with tends to reach the platform stage. This means that sequencing quantity for each sample was sufficient and reliable (Figure 1b).

Figure 1. Quality Control and Basic Analysis. (a) the abscissa is the sequence length of tags, and the ordinate is the number of tags.(b) Rarefaction curves for ASVs.

Microbial diversity and richness among the three groups

The microbial α diversity and β diversity were applied to analyze the microbiota biodiversity and composition among the three groups. Our study compared tongue coating microbiota richness and diversity among Con, M and S-C groups. In the α-diversity analysis, the Chao 1 showed that a significant higher richness of tongue coating microbiota in both S-C group and M group compared with Con group (S-C VS Con, p = 0.0000, M VS Con, p = 0.0062, Figure 2a). Compared with Con group, the Shannon index of the M group showed a decreasing trend (p = 0.87, Figure 2b). Whereas, the S-C group owned significantly lower Shannon index, in comparison with M group (p = 0.0014) and Con group (p = 0.0007).

Figure 2. The microbial α diversity and β diversity analysis in different groups. (a) Chao 1 index was higher in M and S-C group than in Con group. (b) The S-C group had significantly lower Shannon index, in comparison with M group and Con group. (c) A Venn diagram displayed the overlaps among groups. (d) PLS-DA revealed different microbial community structures in the three groups.

In addition, the Venn diagram indicated that 1270 of the total 6209 ASVs were shared among the three groups, with 175, 589 and 773 OUTs were unique for control, M and S-C group respectively (Figure 2c). About β diversity, principal coordinates analysis (PCoA) of the Bray – Curtis distance was performed to visualize the bacterial composition dissimilarity among the three groups. The result demonstrated that the microbiota from tongue coating samples differed among Con, M and S-C groups (Figure 2D).

Changes in tongue coating microbiota composition among the three groups

As depicted in Figure 3a-b, distinct bacterial compositions were observed at the phylum and genus levels among the different groups. We examined the distribution of taxa at these levels to uncover the unique characteristics of each group. The five dominant bacterial phyla in all groups were Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria, and Fusobacteria.

Figure 3. Comparison of relative abundance among each group. Barplots of the relative abundance of the main bacterial taxa at (a, b) phylum and (c, d) genus level for Con, M and S-C groups. (*<0.05, **<0.01).

In the healthy control group, the predominant phyla consisted of Firmicutes (19.3%), Bacteroidetes (36.2%), Proteobacteria (24.6%), Actinobacteria (3.7%), and Fusobacteria (13.1%), with an additional ~ 3.1% of other bacteria (Figure 3a). At the genus level, Prevotella (27.0%) was the main contributor to the microbiota profile, followed by Streptococcus (6.2%), Neisseria (17.9%), Veillonella (8.2%), Rothia (1.9%), Leptotrichia (4.9%), Fusobacterium (8.2%), and Alloprevotella (5.5%) (Figure 3c).

The relative abundance of Firmicutes and Actinobacteria was higher in the mild to moderate (M) and severe and critical (S-C) groups compared to the healthy control (Con) group (p = 0.000). Conversely, Bacteroidetes and Fusobacteria were more abundant in the Con group (p = 0.000) (Figure 3b).

In comparison to the Con group, the M and S-C groups exhibited decreased levels of Prevotella (p = 0.0010), Neisseria (p = 0.0010), Fusobacterium (p = 0.0010), and Alloprevotella (p = 0.0040), along with an overexpression of Streptococcus (p = 0.0070) and Rothia (Figure 3d).

Characterized microbial taxa associated with COVID-19 patients

We accomplished the LEfSe analysis to distinguish the significant tongue coating bacteriome among the three groups. The bar graph and cladogram of different taxonomic levels with LDA scores > 3 were presented in Figure 4, respectively. At the genus level, increased Pseudomonas, Acinetobacter, Lactobacillus, Corynebacterium, Rothia were detected as the top five biomarkers in S-C group. We observed a high level of Streptococcus, Veillonella, Leptotrichia and Actinomyces at the genus level in M group. The enriched species which were found in Con group included Prevotella, Neisseria, Fusobacterium, Alloprevotella, Aggregatibacter, and Porphyromonas.

Figure 4. LEfSe analysis showed the most abundant taxa from the phylum to the genus level among Con, M and S-C groups.

Functional analysis of tongue coating microbiota among the three groups

Finally, Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) was conducted to predict the metagenomes and identify the KEGG pathways involved in each group.

Compared to Con group, M group showed outstanding upregulations of microbial genes involved in the protein families of signaling and cellular processes, carbohydrate metabolism, membrane transport, signal transduction, sorting and degradation, folding and cellular community-prokaryotes, while down-regulated genes involved in the protein families of genetic information processing, glycan biosynthesis and metabolism, metabolism of cofactors and vitamins, metabolism, energy metabolism, replication and repair, poorly characterized, translation and lipid metabolism (Figure 5a).

Figure 5. The function prediction of the three groups. Differential KEGG pathways were analyzed using PICRUSt for three groups. Significant differences between Con and M group (a), Con and S-C group (b), and M and S-C group (c) were presented respectively.

Microbiota associated to S-C group were characterized by a higher potential for Protein families: signaling and cellular processes, signal transduction, membrane transport, xenobiotics biodegradation and metabolism, amino acid metabolism, cellular community-prokaryotes, carbohydrate metabolism and lipid metabolism, while Con group showed robustly reduced Protein families: genetic information processing, glycan biosynthesis and metabolism, translation, Protein families: metabolism, metabolism of cofactors and vitamins, replication and repair, energy metabolism, nucleotide metabolism, sorting and degradation and folding (Figure 5b).

Moreover, when comparing M group and S-C group, S-C group-associated microbiota showed a significantly increased Protein families: signaling and cellular processes, signal transduction, xenobiotics biodegradation and metabolism, amino acid metabolism, cellular community-prokarayotes, membrane transport, poorly characterized and lipid metabolism. Conversely, it displayed a consistently decreased Protein families: genetic information processing, translation, Protein families: metabolism, glycan biosynthesis and metabolism, replication and repair, metabolism of cofactors and vitamins, nucleotide metabolism, energy metabolism and folding (Figure 5c).

Discussion

COVID-19 is an acute respiratory infectious disease caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The ongoing and unpredictable duration of the epidemic, coupled with the high risk of infection, poses a significant threat to the lives of patients [10]. Some individuals with COVID-19 May rapidly develop hypoxemia and/or dyspnea, while those with severe cases may experience sepsis, acute respiratory distress syndrome, multiple organ failure, and challenging-to-correct metabolic acidosis [11]. Currently, there is a limited repertoire of drugs available for the treatment of COVID-19 in clinical practice, with treatment approaches often focused on symptomatic relief, anti-infection measures, antiviral interventions, among others [12]. Consequently, the early differentiation of disease severity among COVID-19 patients holds paramount importance in determining patient prognosis. Oral bacteria are known pathogens in various periodontal diseases, contributing to systemic inflammation and exhibiting associations with systemic conditions such as diabetes, lung disease, and rheumatoid arthritis [13–15]. Growing evidence indicates a correlation between the oral microbiome and the development of COVID-19 [16,17]. Recent research has identified a decline in oral microbial diversity and an increase in dysbiotic species as potential predictors of COVID-19 [6]. This highlights the potential utility of the oral microbiome as a diagnostic tool for SARS-CoV-2 infection. However, studies examining the link between dysbiotic oral microbiota and COVID-19 patients with varying disease severity have been limited.

To investigate whether COVID-19 patients with varying disease severity exhibit distinct microbial community profiles compared to healthy controls, we employed 16S rRNA gene sequencing to analyze the tongue coating microorganisms. Our results revealed a significantly higher richness of tongue coating microbiota in both the mild to moderate (M) and severe and critical (S-C) groups compared to the control (Con) group. However, no statistically significant differences were observed between the M and S-C groups. Additionally, we observed lower microbial diversity in the S-C group compared to the Con group, which is consistent with previous studies [6].

At the phylum level, both the mild to moderate (M) and severe and critical (S-C) groups exhibited a higher relative abundance of Firmicutes and Actinobacteria, while showing a lower abundance of Bacteroidetes and Fusobacteria compared to the control (Con) group. Notably, the lower Bacteroidetes/Firmicutes ratio observed in COVID-19 patients was not as pronounced. Previous studies have reported a lower Bacteroidetes/Firmicutes ratio in the gut microbiota of obese individuals, suggesting it as a potential biomarker [18–20]. However, its association with inflammation in COVID-19 patients remains unstudied. This discrepancy may reflect the unique pathogenesis of COVID-19.

In the genus-level analysis, both the mild to moderate (M) and severe and critical (S-C) groups exhibited significantly increased proportional abundance of Streptococcus and Rothia, while showing decreased abundance of Prevotella and Neisseria compared to the control (Con) group. Streptococcus encompasses symbiotic, pathogenic, and opportunistic aerobic Gram-positive organisms in humans and animals [21], primarily inhabiting body cavities and the skin [22]. Previous studies have also identified Streptococcus as an indicator of ventilator-associated pneumonia or COVID-19 severity [23,24]. Elevated levels of Streptococcus have been associated with increased expression of inflammatory cytokines such as IL-18, TNF-α, IFN-γ, leading to worsened clinical outcomes [25–27]. Additionally, the increased abundance of Streptococcus has been linked to up-regulated adhesion receptors that facilitate viral entry during respiratory virus infections [28].

Rothia, a Gram-positive aerobic bacterium, is commonly found in the oral and respiratory tracts [29]. While typically considered commensal, Rothia can act as an opportunistic pathogen causing various infections, particularly in individuals with compromised immune function [30]. A study has shown an association between Rothia and susceptibility to secondary bacterial pulmonary infections in patients infected with the avian H7N9 virus [31]. Notably, the role of Rothia in the pathogenesis and disease severity of SARS-CoV-2 infection, as well as its potential impact on the susceptibility of COVID-19 patients to secondary bacterial lung infections, requires further investigation.

In the study of Jitvaropas et al. [32], their LefSe analysis found that the phylum Proteobacteria, class Gammaproteobacteria, order Burkholderiales, family Burkholderiaceae, genus Lautropia, and species Lautropia mirabilis were significantly dominant in the asymptomatic COVID-19-positive group. Our LEfSe analysis pointed out that at the genus level, Pseudomonas, Acinetobacter, Lactobacillus, Corynebacterium, Rothia were detected as the top five biomarkers in S-C group, Streptococcus, Veillonella, Leptotrichia, Actinomyces and Chloroplast were detected as the top five biomarkers in M group while Prevotella, Neisseria, Fusobacterium, Alloprevotella, Aggregatibacter, and Porphyromonas were biomarkers for Con group.

Apart from compositional changes in bacterial taxa, we also predicted alterations in function among the groups. In our study, we observed a rise of Protein families: signaling and cellular processes among the three groups, which suggested a potential underlying mechanism of COVID-19.

In conclusion, our data provide evidence that some tongue coating microbes may be related to the severity of COVID-19 patients, and these associations could have guiding significance for the eradication or vaccination prevention strategies of high-risk groups. Prospective research is still needed, followed by animal experiments to clarify the causal relationship between them.

However, this study has two important limitations that need to be acknowledged. Firstly, the relatively small sample size in each group restricts the generalizability of our findings, and larger studies with more participants are needed to obtain more robust and reliable results. Secondly, as a cross-sectional study, it emphasizes the need for future prospective trials to establish a causal relationship and fully elucidate the role of the microbiota in determining the severity of COVID-19.

Authors’ contributions

Zongdan Jiang, Zhenyu Zhang and Wanli Liu conceived, organized and supervised the project, and proofread the manuscript. Zongdan Jiang, Lu Yang and Qu Yi collected and analysed the data. Zongdan Jiang and Xuetian Qian drafted the manuscript. Kunhan Su and Yuzhen Huang supervised statistical analysis. All authors approved the final version of the manuscript. All the authors critically revised and approved the final version.

Availability of data and material

All raw reads were stored in NCBI Sequence Read Archive (SRA) database, and the accession number is PRJNA970271.

Ethics approval

The study protocol was approved by the institutional review board of Nanjing Integrated Traditional Chinese and Western Medical Hospital, and all the experiments were performed in accordance with approved guidelines and regulations.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was supported by Nanjing soft science Research Program (202206029). This work also was supported by Nanjing Science and Technology Information Research Institute Project (NJKJZX2022A002).