This study was approved by the Ethics Committee of the Airforce Medical Center and the Ethics Committee of Peking University School and Hospital of Stomatology (No. PKUSSIRB-202164061). To estimate the sample size, we used the RnaSeqSampleSize method implemented in the R language by Zhao et al. which was based on the distributions of average gene read counts and dispersions estimated from real RNA-seq datasets of The Cancer Genome Atlas (TCGA) (Zhao et al., 2018). The estimated minimum fold change in prognostic gene expression between the two groups was set to 2.5 based on a former study investigating the salivary microbiome of OSA patients and normal controls (Jia et al., 2020), and the sample size of 50 would reach a false discovery rate (FDR) of 0.01 with a total of 20,000 genes for testing and 200 prognostic genes. Considering possible complications of the study, 60 participants were initially recruited.

From July to December 2019, patients aged 25–35 years who were diagnosed with obstructive sleep apnea were included from two sleep medicine centers in Beijing. Then, non-OSA controls matched with the OSA group by age, gender and BMI were voluntarily recruited from an enterprise’s male dormitory in Beijing. The inclusion criteria were as follows: 1) age >18 years; 2) male; 3) self-reported as having no alcohol abuse; and 4) self-reported as never smoking or having quitted smoking for more than 6 months. All participants signed informed consent after being told study details.

The exclusion criteria included a special diet preference, pet feeding, systemic diseases (hypertension, diabetes, etc.), salivary gland disease, oral mucosa infectious disease, untreated decays or pericoronitis, application of antibiotics or hormone drugs in the last 3 months, and periodontal and/or OSA treatment in the last 1 month. Among the 60 initially recruited participants, two voluntarily withdrew the periodontal examination, one failed to record sleep monitoring data due to connection failure, and three had saliva samples containing less than 0.5 ng/μl DNA. Finally, 54 participants, consisting of 27 OSA patients and 27 normal controls, completed the whole study with a dropout rate of 10%.

The periodontal examination was practiced for each participant by one experienced double-blinded dentist before polysomnography (PSG). A full mouth probing was conducted to record the probing depth (PD) and clinical attachment level (CAL) of each participant. All measurements were obtained using a Williams style periodontal probe (Hu-Friedy, Inc., Chicago, IL, USA) and lengths were rounded up or down to the nearest millimeter. CAL is determined as follows: 1) If the gingival margin is located on the anatomic crown, CAL is calculated by PD minus the distance from the gingival margin to the cementoenamel junction (CEJ). 2) If gingival margin coincides with the CEJ, CAL equals to PD. 3) If the gingival margin is located apical to the CEJ, CAL is PD plus the distance between the CEJ and the gingival margin.

In the diagnosis of periodontitis, we used the 2017 World Workshop on the Classification of Periodontal and Peri-implant Diseases and Conditions. A patient was a periodontitis case if: 1) Interdental CAL is detectable at ≥2 non-adjacent teeth, or; 2) Buccal or oral CAL ≥3 mm with pocketing >3 mm is detectable at ≥2 teeth (Tonetti et al., 2018).

Each participant underwent standard overnight PSG using a portable computerized device (SOMNOscreen™ plus, SOMNOmedics GmbH, Randersacker, Germany).

Before PSG, the blood pressure and heart rate were recorded. Accordingly, central electroencephalogram and electrocardiogram were taken. Chest and abdominal movement were traced by respiratory effort bands. Body position and oronasal airflow were recorded by pressure sensors, and oxygen saturation (SaO₂) was measured by a pulse oximeter. In the morning, blood pressure and heart rate were recorded again within 1 h after the rising.

Sleep stages and arousals, movements, and cardiopulmonary events were manually scored by one experienced sleep technician based on the scoring manual of the American Academy of Sleep Medicine (AASM) (Berry et al., 2020). The standard definition of an apneic event includes a minimum 10-second interval between breaths, with a neurologic arousal, a blood oxygen desaturation of 3% or greater, or both arousal and desaturation. Hypopnea is defined as an episode of shallow breathing (airflow reduced by ≥50%) during sleep, lasting for ≥10 s and usually associated with a fall in blood oxygen saturation attributable to partial obstruction of the upper airway. The apnea–hypopnea index (AHI) is calculated as the total number of apneas and hypopneas per hour of sleep. In adults, the presence of OSA was determined as an AHI score of ≥5, while mild OSA is defined as an AHI of at least 5 to 15 events per hour, moderate OSA as >15 to 30 events per hour, and severe OSA as >30 events per hour (Sateia, 2014).

Saliva samples were collected in the morning after the overnight PSG. All participants were required not to drink, eat, brush teeth, or floss for at least 2 h. Fifteen minutes before sampling, they were guided to gently rinse their mouth using purified water. Approximately 2 ml of whole saliva was collected by passive drool using Saliva Collection Aid (Salimetrics, State College, PA, USA). All samples were frozen at −20°C immediately after collection, and transferred to −80°C refrigerators in no more than 4 h until further processing.

Before DNA extraction, all samples were thawed at room temperature. Approximately 500 μl saliva of each sample was transferred into the centrifuge tube and centrifuged at 10,000 rpm for 15 min. Only the residues were collected and resuspended in 300 μl phosphate buffered saline (PBS). Then we used QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA) to extract DNAs according to the manufacturer’s protocols. The quality of the extracted DNA was checked by the spectrophotometer (Nanodrop 8000, Thermo Fisher Scientific, Waltham, MA, USA) to see whether the OD_260/280 ratio was 1.8–2.0. The integrity of DNA was verified by 1% agarose gel electrophoresis, where samples with less than 0.5 ng/μl DNA were excluded. The high-quality DNA was stored at −20°C for further sequencing.

Isolated DNA was used as a template for polymerase chain reaction (PCR) (Ko et al., 2019). The V3–V4 regions of the bacterial 16S rRNA gene were amplified with 341’F primer (CCTAHGGGRBGCAGCAG) and 805’R primer (GACTACHVGGGTATCTAATCC) (Hugerth et al., 2014) using a KAPA HotStart ReadyMix PCR Kit (KAPA Biosystems, Wilmington, MA, USA). The PCR program was set as follows: 98°C for 2 min, 26 cycles of 98°C for 20 s, 54°C for 20 s, 72°C for 2 min (Eriksson et al., 2019). The products were purified using an AxyPrep DNA Gel Extraction Kit (Axygen, Union City, CA, USA). The libraries were sequenced using the Illumina NovaSeq™ 6000 platform (Illumina Inc, San Diego, CA, USA) by Shanghai Personalbio Co., Ltd. (Shanghai, China) and all sequencing data were submitted to the NCBI Short Reads Archive Database under accession number SRP339526.

To process the raw sequencing data, the Vsearch (v2.13.4_linux_x86_64) software was used. First, amplicon reads with low quality and primers were trimmed and filtered, followed by chimeras removed and singletons filtered. Then selected high-quality sequences with a similarity threshold of 97% were clustered to the same operational taxonomic units (OTUs) using UPARSE (Edgar, 2013). In each OTU, the most abundant sequence was used as the representative sequence and was taxonomically classified into microbial taxa (phylum, class, order, family, and genus) according to the Ribosomal Database Project database and the Human Oral Microbiome Database (Dewhirst et al., 2010).

To compare the diversity or complexity of species in OTUs from OSA and non-OSA samples, alpha diversity analysis was applied through several indices, namely, observed species, Chao, Shannon, Simpson, Pielou’s evenness and goods coverage. The indices were calculated by Mothur (v1.31.2), and one-way ANOVA analysis was used for multigroup comparison. In beta diversity analysis, we used the unweighted UniFrac distances matrices to perform the principal coordinates analysis (PCoA) and analysis of similarities (ANOSIM), in order to ordinate the dissimilarity matrices of the samples.

To analyze the difference of species components in each group, the number of total OTUs, OTUs in the OSA and the non-OSA group was shown by a Venn diagram. Then, the distribution of OTUs in each sample were plotted in one histogram using R, version 3.0.3 (R Foundation for Statistical Computing, Vienna, Austria). Kruskal–Wallis tests were conducted to determine which taxonomic groups were significantly different between groups of samples.

Multivariate analyses using Multivariate Association with Linear Models (MaAsLin) was performed to test associations between OSA and salivary microbiome, taking periodontitis, age, BMI, and alcohol use into account as possible confounding factors. In the MaAsLin models, periodontitis (yes/no) and alcohol use (current/former/never) were defined as categorical metadata variables, along with age and BMI defined as continuous variables. A q-value <0.10 was used as significance cut-off for this taxa-phenotype associations (van der Meulen et al., 2019).

All statistical analysis, if not specified, was conducted using IBM SPSS Statistics for MacOS, version 26.0 (IBM Corp., Armonk, N.Y., USA), and a two-sided p-value of <0.05 was considered to be statistically significant. A normality test was practiced on all variances we collected.

We set OSA as the risk factor and periodontitis as the outcome, and the χ² test was applied to measure the difference in prevalence of periodontitis between variables. The differences between the mean values of PD were then evaluated between two groups using the independent t-test.

The comparison of sequencing data between two groups were firstly tested for the homogeneity of variance, and then followed by one-way ANOVA or Wilcoxon rank sum test.

In total, 54 male participants were enrolled in the study, with 27 participants in the OSA group and 27 participants in the non-OSA group according to the PSG diagnosis. The two groups were comparable with respect to other characteristics, as age, alcohol use, annual incomes, and level of education were basically matched. Additional participants demographics are presented in Table 1 .

With regard to the PSG diagnosis, almost all participants could be considered normal range (51.9% had BMI from 18.5 to 24.9 kg/m²) or pre-obese (46.3% had BMI from 25 to 29.9 kg/m²) according to the WHO criteria for obesity. A significantly higher BMI and a higher proportion of pre-obese participants were discovered in the OSA group.

Based on the periodontal examination, 30 participants (55.6%) had been diagnosed with periodontitis. The OSA group showed a significantly higher prevalence of periodontitis, but there was no statistically significant difference of mean PD or mean CAL between the two groups.

After data processing, a total of 4,133,794 high-quality reads were generated from 54 salivary samples (ranging from 57,599 to 112,104 reads), with an average of 76,552 ± 10,625 reads per sample. Then, all reads were clustered by using a 97% similarity cutoff. A total of 483 OTUs, with 341 and 434 OTUs in the OSA and the non-OSA group respectively were finally detected. As shown by the species accumulation curve ( Figure 1A ) and the rarefaction curve ( Figure 1B ), the curve gradually became flat with the increasing number of sequencing samples, indicating that the sample size is sufficient enough to represent the species composition as well as the diversity of microorganisms in the samples.

The alpha diversity of microbial community was calculated and then compared between the OSA group and the non-OSA group. The OSA group had a significantly lower microbial species richness than the controls, as demonstrated by the Chao index (p = 0.027) ( Figure 2A ) and observed species (p = 0.033) ( Figure 2B ). The OSA group also showed higher goods coverage (p = 0.023) ( Figure 2E ). Shannon ( Figure 2C ) and Simpson ( Figure 2D ) indexes indicated that there was no statistical difference of the microbial species diversity between the two groups. The microbial species evenness estimator (Pielou’s evenness) ( Figure 2F ) also showed no significant difference.

Beta diversity analysis was conducted to characterize the similarities or differences of microbial communities. PCoA to unweighted UniFrac distances were conducted in four categories: OSA/periodontitis, OSA/non-periodontitis, non-OSA/periodontitis, and non-OSA/non-periodontitis. A separation was observed among the four categories. Scatters of the OSA/non-periodontitis gathered to the left top of the plot, while scatters of the non-OSA/non-periodontitis category gathered to the right bottom ( Figure 3A ).

The distinction of microbial communities between the two groups were further tested by ANOSIM analysis ( Figure 3B ), showing a significant higher between-habitat diversity at the OSA group than that at the non-OSA group (R² = 0.065875, p = 0.005).

As shown by the Venn diagram ( Figure 4A ), the number of total OTUs in the OSA group was fewer than the non-OSA group (341 versus 434). There were 49 unique OTUs in the OSA group and 142 in the non-OSA group.

The taxonomic composition in samples of genus level was then generated into the heatmap ( Figure 4B ) to visualize the distribution of bacteria in each sample. While analyzing the detected bacteria at phylum, genus, and OTU levels, the microbial distribution of the top 10 taxa is provided in Figure 4C . The main detected bacteria phylum in the two groups included Firmicutes, Actinobacteria, and Proteobacteria, which made up more than 70% of the whole bacteria in each group. The top-most genus included Streptococcus, Neisseria, and Actinomyces. It was also observed that the microbial community structure differed from the OSA group to the non-OSA group. Among these top abundant microbes, Rothia (8.73% versus 5.01%) and Gemella (3.58% versus 2.56%) showed a significant abundance in the OSA group. In contrast, the abundance of Streptococcus (4.07% versus 2.17%) and Veillonella (4.07% versus 2.17%) was relatively higher in the non-OSA group. In OTU level, R. mucilaginosa (7.39% versus 3.92%) was significantly abundant in the OSA group, while the abundance of P. melaninogenica (4.36% versus 2.12%) was higher in the non-OSA group.

In the MaAsLin models, OSA was statistically associated with Peptostreptococcus and Granulicatella (q <0.1). No statistical associations were detected between OSA and periodontitis, age, BMI, and alcohol use.

At the genus level, the relative abundance of eight predominant genera which had been proved to initiate periodontal disease was visualized in the plot ( Figure 5A ). Further statistical analysis showed that only Prevotella were observed to have significant higher abundance in the OSA group than that in the non-OSA group. The relative frequency of other seven genera were found to have no significant difference between the OSA group and the non-OSA group ( Table 2 ).

At the species level, we included five specific species which have a clear proof of periodontal pathogenesis in this study, namely, Aggregatibacter actinomycetemcomitans (Aa), Prevella intermedia (Pi), Campylobacer rectus (Cr), Fusobacterium nucleatum (Fn) and Treponema denticola (Td). As shown by the box diagram ( Figure 5B ), the relative abundance of Aa, Pi, and Td was slightly higher in the OSA group. However, the differences of relative abundances of these five species related to periodontal diseases were not significant between the two groups ( Table 3 ).