This was an observation and comparative study. The participants were prospectively enrolled in an investigation of tonsil and gut microbiomes at the Department of Otolaryngology, Chang Gung Memorial Hospital (Linkou Main Branch, Taoyuan, Taiwan) between March 2017 and January 2019 (33). The Institutional Review Board of Chang Gung Medical Foundation approved the study (Approval No.: 201507279A3). Both parents and participants aged 6 years or older provided written informed consent. Our research adhered to the revised Declaration of Helsinki (34) and complied with the STROBE guidelines (35).

Children aged 5-12 years, exhibiting an AHI of ≥ 5.0 events/hour or an AHI of ≥ 2.0 events/hour with at least one associated morbidity (e.g., elevated blood pressure, daytime sleepiness, growth retardation), were considered eligible (36). We excluded patients with craniofacial, neuromuscular, or chronic inflammatory disorders (37). Those with acute inflammation or conditions requiring antibiotic treatment were only eligible for stool sample collection after a minimum of two weeks following remission.

From a database of 66 children (33), we meticulously selected six children with OWO and another six with NW based on age, sex, and AHI. The participants were subsequently categorized based on their BMI z-scores into the ‘OWO’ group (BMI z-score ≥ 1.0) and the ‘NW’ group (BMI z-score > -2.0 and < 1.0) (38). Figure 1 displays the study flow diagram.

Our study assessed various factors linked to OSA severity, including AHI, apnea index, arousal index, mean blood oxygen saturation (SaO₂), minimum SaO₂, total sleep time, and sleep stages. Following the 2012 American Academy of Sleep Medicine Manual guidelines (39), we employed standard full-night in-lab polysomnography. The detailed protocol for this polysomnography has been previously documented (37, 40, 41).

Stool samples were collected, snap-frozen in liquid nitrogen, and stored at -80°C. DNA extraction utilized a fecal DNA isolation kit (MoBio Laboratories, USA). DNA concentrations were ascertained with a NanoPhotometer P360 system (Implen, USA) and standardized to 1 ng/ul for subsequent analyses of full-length or 5 ng/μl for subsequent analyses of partial length 16S rRNA amplicon sequencing.

The full-length 16S genes encompassing V1–V9 hypervariable regions were amplified using barcoded 16S gene-specific primers as detailed in PacBio^®’s guide (42). Each primer is designed to contain a 5’ buffer sequence (GCATC) with a 5’ phosphate modification, a 16-base barcode and the degenerate 16S gene-specific forward or reverse primer sequences (Forward:5’Phos/GCATC- 16-base barcode - AGRGTTYGATYMTGGCTCAG -3’, Reverse: 5’Phos/GCATC- 16-base barcode – RGYTACCTTGTTACGACTT -3’) (43). The degenerate base identities are defined as follows: R = A or G; Y = C or T; M = A or C. A total of 2 ng of gDNA was used for the polymerase chain reaction (PCR), employing the KAPA HiFi HotStart ReadyMix (Roche, USA) under specified PCR conditions. Post-reaction, the PCR products were examined on a 1% agarose gel. Samples exhibiting a prominent band around 1,500 bp were selected and purified using AMPure PB Beads for subsequent library preparation.

The SMRTbell library preparation followed the PacBio^®’s guidelines as above. Briefly, an equal molar amount of each barcoded PCR product was pooled, and between 500-1,000 ng of the pooled amplicon sample underwent DNA damage repair. This was followed by end-repair, A-tailing, and ligation steps to attach the universal hairpin adapters to double-stranded DNA region. After purification with AMPure PB beads to discard adapter dimers, the SMRTbell library was prepped with Sequel II primer 3.1 and Sequel II Binding Kit 3.1 for primer annealing and polymerase binding. Finally, sequencing was executed in the circular consensus sequence (CCS) mode on a PacBio^® Sequel IIe instrument to yield HiFi reads with a predicted accuracy (Phred Scale) of 30.

CCS reads were determined based on a minimum predicted accuracy of 0.9, with the least number of passes set at three, as per PacBio^®’s official workflow using the SMRT Link software. Only CCS reads exceeding a quality score of Q30, termed Q30 HiFi reads, advanced to the succeeding phase. Post-demultiplexing, the HiFi reads underwent further processing using DADA2 (version 1.20; https://qiime2.org/) to extract amplicons with single-nucleotide resolution (44). Trimming and filtering were set to a maximum of two expected errors per read. The DADA2 algorithm can discern exact amplicon sequence variants (ASVs) from the full-length 16S rRNA gene with near-zero error rates. For each distinct sequence, tools such as the feature-classifier (45) and classify-consensus-vsearch (46) algorithm in QIIME 2 (v2022.11) were employed to annotate taxonomy classifications based on data from the National Center for Biotechnology Information (NCBI) 16S ribosomal RNA database (47). To evaluate sequence similarities among diverse ASVs, a multiple sequence alignment was executed using QIIME 2’s MAFFT (48) tool against the NCBI 16S ribosomal RNA database. A phylogenetic tree, illustrating the relationship of representative ASV sequences, was constructed using QIIME 2’s phylogeny FastTree (49).

We targeted the V3–V4 hypervariable region of the 16S rRNA gene for sequencing, amplified with specific primers: 341F (5’-CCTACGGGNGGCWGCAG-3’) and 806R (5’-GACTACHVGGGTAT CTAATCC -3’). Amplification followed the 16S Metagenomic Sequencing Library Preparation protocol (Illumina^®, USA). PCR was carried out with 12.5 ng of genomic DNA, using KAPA HiFi HotStart ReadyMix (Roche, USA). The PCR products were verified via 1.5% agarose gel electrophoresis. DNA from samples with a prominent ~500 bp band was purified using AMPure XP beads (Beckman Coulter, USA). Quality assessment of the indexed PCR products was conducted using the Qubit 4.0 Fluorometer (Thermo Scientific, USA) and the Qsep100TM system (BiOptic, Taiwan). The sequencing library was constructed as per the aforementioned Illumina^® protocol, and sequencing was executed on an Illumina^® MiSeq platform, yielding paired 300-bp reads.

Post-sequencing, raw reads were demultiplexed based on barcodes. The paired-end reads underwent primer and adapter sequence removal via QIIME 2 cutadapt plugin (50). ASV construction involved the QIIME 2 DADA2 plugin (v2021.4), which facilitated quality filtering, dereplication, denoising, and more (44). Taxonomy classification and phylogenetic tree were performed as above.

To account for sequence depth variations, ASVs were rarefied to the minimal sequence depth using QIIME’s script (single_rarefaction.py). Both alpha and beta diversity analyses utilized this normalized data. The community’s relative abundance and evenness accounting for alpha diversity were assessed using the Shannon, Simpson, and Peilou indices (51). For beta diversity, weighted and unweighted UniFrac were calculated (52). Statistical evaluations employed principal coordinate analysis (53) and non-metric multidimensional scaling (NMDS) (54). Furthermore, Welch’s t-test with the Benjamini-Hochberg procedure was performed to control false discovery rate using the STAMP software (v2.1.3) (55). The presence of statistically significant biomarkers was ascertained using the linear discriminant analysis effect size (LEfSe) analysis (56). Community structure differences were determined using analysis of similarity (57) and permutational multivariate analysis of variance (58). Functional abundances were predicted from full-length 16S rRNA data using Tax4Fun2 and the V3–V4 region for the Tax4Fun2 (59) and Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology database (60). For highlighting pivotal functional profiles, the ggplot2 and microeco R packages were coupled with metastat statistics (61).

The normality of the continuous variables was assessed using the Kolmogorov-Smirnov test. For variables that followed a normal distribution, results were presented as the mean with standard deviation (SD). In cases where the distribution was non-normal, we reported the median alongside the interquartile range (IQR). To evaluate the differences between the two weight statuses, the unpaired Student t-test or Mann-Whitney U test were employed for continuous variables, while the Fisher exact test was used for categorical variables. For within-group comparisons, we applied either the paired Student t-test or the Wilcoxon signed-rank test, depending on the nature of the continuous variables. All statistical procedures were performed using multiple software tools: R (versions 4.3.1, R Foundation for Statistical Computing, Vienna, Austria), SPSS (version 27.0, IBM Corp., Armonk, NY, USA), and GraphPad Prism 10.0 for Windows (Graph Pad Software Inc., San Diego, CA, USA).

Of the 66 children initially diagnosed with OSA, 36 were excluded because their stool samples were collected after December 31, 2017. Of the 30 children that remained, 17 had both OSA and OWO, and 13 had OSA with NW. Further exclusion based on age, sex, and weight status criteria left 12 participants for the final analysis: 6 children with concurrent OSA and OWO, and another 6 with OSA and NW, as illustrated in Figure 1 . The final cohort consisted of 10 boys and 2 girls, with a median age of 6 years (IQR: 5–10 years). This is detailed in Table 1 . The clinical variables were similar between the OWO and NW groups, except for the BMI z-score, which showed a significant difference (p = 0.002).

Table 2 outlines the ASV counts, and corresponding sequence counts for each sample. With a Q score threshold of 30, the full-length 16S rRNA sequence via PacBio^® yielded significantly fewer ASVs than the V3–V4 region sequenced with Illumina^® (16044 (15554–18552) vs. 147952 (141946–152740), p < 0.001). Notably, almost all ASVs derived from full-length 16S rRNA sequencing were classified, resulting in a negligible rate of unassigned sequences (0–0.2% range) across phylum to species levels. In contrast, the V3–V4 region had a progressive increase in unassigned ASVs, especially from family (0.8%–14.5%) to species tiers (97.9%–99.9%). Moreover, the coverage rates from full-length V1–V9 region surpassed those of V3–V4 across all taxonomic levels (all p < 0.01).

Upon normalization by relative abundance, there were no significant disparities in the taxonomic profiles between full-length 16S rRNA sequencing and V3–V4 region from the phylum to the genus levels ( Figures 2A–E ; all p > 0.05). However, distinctions emerged at the species level ( Figure 2F ; p < 0.001). In terms of alpha diversity, the Simpson and Pielou indexes were both appreciably higher for full-length 16S rRNA sequencing as compared to the V3–V4 region ( Figure 3A ; p = 0.021 and p = 0.001, respectively). Delving further, it was observed that both unweighted ( Figure 3B ) and weighted ( Figure 3C ) UniFrac distances, as computed for individual samples from patients with NW and OWO, exhibited pronounced differences between the two NGS methodologies (p = 0.001 and p = 0.002, respectively). Regarding beta diversity ( Figure 3D ), the unweighted UniFrac distances associated with the full-length 16S rRNA sequencing were considerably lower than those of the V3–V4 region (p = 0.001). Conversely, the weighted UniFrac distances recorded for the full-length 16S rRNA sequencing surpassed those for the V3–V4 region (p = 0.001).

There were notable differences in ASVs at both the phylum and genus levels between the OWO and NW groups, as evidenced by both the full-length 16S rRNA sequencing and the V3–V4 region methodologies ( Figures 4A, B ; all p < 0.05). However, the full-length method provided deeper insights compared to the V3–V4 region. This distinction was apparent when evaluating the F/B ratio ( Figure 4C ), the alpha diversity ( Figure 4D ), beta diversity ( Figure 4E ), and the number of differentiated ASVs ( Figures 5A, B ), helping to better distinguish children with OWO from those with NW.

Using the Tax4Fun2-based NMDS plot, a marked differentiation in the habitat-specific functional profile of gut microbiota emerged between children with OWO and those of NW status. This was clearly highlighted when analyzing data via the full-length 16S rRNA sequencing, paired with the level 3 pathway of the KEGG orthology database (stress value = 0.065) ( Figure 6A ). In parallel, the Tax4Fun NMDS plot, derived from V3–V4 segment data and aligned with the same functional prediction pathway, particularly pronounced divergences in functional profiles between the two weight groups (stress value = 0.021) ( Figure 6B ).

By harnessing the Welch’s t-test in tandem with the Benjamini-Hochberg procedure to manage the false discovery rate, we pinpointed 73 distinct level 3 pathways using full-length 16S rRNA genes and another 61 pathways using the V3–V4 region ( Figure 6C ). All pathways were identified with both p- and q-values falling below 0.05. Intriguingly, both NGS methods shared 18 predictive functions. Difference in the proportion of overlapping level 3 functional predictions were not significant (25% vs 30%, p = 0.562). Figure 7 delves deeper into these shared 18 level 3 pathways. While the majority (n = 16) exhibited consistent distribution patterns across both methods, two pathways — specifically, fructose and mannose metabolism, and glycosaminoglycan biosynthesis-chondroitin sulfate — displayed distinct patterns between the two methods.

Additionally, of the 73 identified through full-length methodology, 35 (48%) were affiliated with the level 1 metabolism pathway, while from the 61 pathways discerned through the V3–V4 region, 28 (46%) resonated with the same category. Interestingly, five metabolism pathways—comprising fructose and mannose metabolism, galactose metabolism, tetracycline biosynthesis, steroid hormone biosynthesis, and pentose and glucuronate interconversions—exhibited consistency between the two methodologies. The proportion of overlapping level 1 metabolism pathways were similar (86% vs 82%, p = 0.737).

Utilizing the Spearman correlation test, a set of bacteria—Bacteroides ovatus, Bifidobacterium pseudocatenulatum, and Streptococcus parasanguinis ATCC 15912—were identified in association with BMI z-score. These were discerned through the combined use of Welch’s t-test and the Benjamini-Hochberg procedure. Additionally, Bacteroides uniformis was singled out via the LEfSe method as having relevance to BMI z-score. Concurrently, an analysis using Tax4Fun2 highlighted 17 metabolic pathways that bore associations with the BMI z-score ( Figure 8 ).

The intricate interconnections amongst BMI z-score, the four gut bacterial species, and the metabolic pathways are vividly depicted in the network graph. This representation underscores the multifaceted associations between BMI z-score and these critical components of the gut microbiome. Strikingly, the trio of bacteria—Bacteroides ovatus, Bifidobacterium pseudocatenulatum, and Streptococcus parasanguinis ATCC 15912—showcased active interactions with the metabolic pathways ( Figure 9 ).

Comprehensive 16S rRNA gene sequencing revealed that alpha diversity indices were not correlated with polysomnographic measures. However, distinct genera demonstrated significant associations: Bifidobacterium levels correlated with AHI and Escherichia with AI. Additionally, Phocaeicola abundance was inversely related to AI, and Bifidobacterium abundance was inversely associated with minimal SaO₂. Extended analysis using the V3–V4 regions of the 16S rRNA gene confirmed these relationships and further identified a significant correlation between Bifidobacterium and AI ( Table 3 ).

Metabolic function predictions based on full-length 16S rRNA gene sequencing indicated significant associations between microbial metabolism in diverse environments, propanoate metabolism, and AI, as well as between glycine, serine, and threonine metabolism and mean SaO₂. Conversely, sequencing of the V3–V4 regions elucidated a notable association of peptidoglycan biosynthesis with minimal SaO₂ ( Table 4 ).

In this research, we focused on comparing two sequencing methods: full-length sequencing and the V3–V4 region sequencing of the 16S rRNA gene. One of the most striking observations was that full-length sequencing demonstrated superior taxonomic coverage compared to the V3–V4 region as previous studies (27, 62). This suggests that the full-length method captures a broader spectrum of microbial diversity, potentially providing a more comprehensive picture of the microbial community. This comprehensive taxonomic resolution could be particularly beneficial in studies where fine-scale differences in microbial communities are of interest (63).

While the full-length sequencing yielded fewer ASVs after denoising, it’s crucial to consider the quality and relevance of these ASVs. Fewer ASVs doesn’t necessarily imply reduced microbial diversity but may indicate reduced noise or spurious sequences (64). The reduced rate of unassigned sequences in the full-length method further underscores its precision. Unassigned sequences, often a source of ambiguity in microbial studies, can stem from sequencing errors, chimeric sequences, or genuine novel sequences that aren’t yet represented in reference databases (65).

Despite similar relative abundances of both methods, the significant variations in both alpha (within-sample diversity) and beta diversities (between-sample diversity) between the two methods have substantial implications. The higher alpha diversity in the full-length sequencing might suggest that it can detect a wider range of taxa within a single sample, potentially unearthing rare or low-abundance species that might be missed by the V3–V4 method (66). Differences in beta diversity, on the other hand, could imply that community compositions derived from the two methods might not be directly comparable (22). This is pivotal for studies looking at differences between groups, such as those based on health status, geography, or other variables.

While our results highlight the potential advantages of full-length 16S rRNA gene sequencing, researchers should carefully consider their objectives before selecting a method. For instance, studies that require rapid results or are constrained by budget might still benefit from the V3–V4 region due to its established protocols, faster turnaround, and lower per-case costs (67). Conversely, projects demanding in-depth taxonomic resolution might favor the full-length approach.

The observed superiority of full-length sequencing in discerning the differences in gut microbiota composition and functional prediction between OWO and NW groups versus the V3–V4 approach is a significant finding. The enhanced depth offered by full-length sequencing might be analogous to viewing an image in higher resolution (66). Where the V3–V4 approach might provide a general view of the microbiota landscape, full-length sequencing offers a detailed map. In this study, the performance of the full-length sequencing to differentiate children with OWO from children with NW was better than that of the V3–V4 approach, regarding F/B ratio, alpha diversity, and beta diversity. Therefore, the full-length sequencing is particularly cost-effective when the sample size was relatively small. This can be especially significant when identifying minor bacterial taxa that may have substantial impacts on host physiology (67, 68). For example, the full-length method revealed six bacterial species of 29 representative ASVs, whereas the V3–V4 region methodology identified 0 bacterial species of 13 representative ASVs. Furthermore, Bacteroides ovatus, Bacteroides uniformis, Bifidobacterium pseudocatenulatum, and Streptococcus parasanguinis ATCC 15912 were closely related to BMI z-score by using the full-length method.

The differences in functional predictions identified by each sequencing method underscore the complex and multifaceted nature of the gut microbiome’s role in host metabolism (69). In this study, we pinpointed 73 distinct level 3 pathways using full-length 16S rRNA genes and another 61 pathways using the V3–V4 region. Differences in gut microbiota composition can impact nutrient absorption, gut barrier function, inflammation, and interactions with other organ systems (70). Furthermore, the metabolic pathways, whether involved in short-chain fatty acid production, amino acid metabolism, or other processes, might hold the key to understanding how gut microbiota influence weight and overall health (71–73). Using the full-length method, 35 differentiated metabolic pathways were identified, whereas the V3–V4 region method revealed 28 differentiated metabolic pathways. Furthermore, Bacteroides ovatus, Bifidobacterium pseudocatenulatum, and Streptococcus parasanguinis ATCC 15912 showcased active interactions with the metabolic pathways and closely related to BMI z-score by using the full-length method. Therefore, beyond just identifying which bacterial species are present, the full-length sequencing might offer insights into the potential functions of these bacteria. Moreover, understanding how various strains and species are related can give insights into how the gut microbiome might have evolved in response to dietary, environmental, or other external factors in children with OWO versus children with NW.

The ability of full-length sequencing to associate specific bacterial species with the BMI z-score is not just a technical achievement but also a clinically relevant finding. Such associations might provide the foundation for personalized therapeutic interventions. For instance, seaweed Undaria pinnatifida administration significantly increased lean-related Bacteroides ovatus and short-chain fatty acids and tricarboxylic acid cycle intermediates and decrease high-fat diet-induced body weight gain in high-fat diet-fed mice (74). As mentioned above, carbohydrate, amino acids, and lipid metabolisms are weight-related metabolic pathways. Interestingly, the role of Bifidobacterium pseudocatenulatum on these pathways is strain dependent. For example, Bifidobacterium pseudocatenulatum CECT 7765 can reduce both metabolic and immunological dysfunctions related to obesity in HFD-fed mice (75), whereas Bifidobacterium pseudocatenulatum JCLA3 involves the carbohydrate and amino acids metabolism (76). Streptococcus parasanguinis ATCC 15912 is an atypical viridans streptococcus (77) and its clinical significance is not well studied. In this study, Streptococcus parasanguinis ATCC 15912 is positively associated with BMI z-score and involves amino acid and lipid metabolic pathways. This is a novel finding; however, its exact mechanism of weight gaining should be further investigated. Nevertheless, understanding the metabolic pathways prevalent in different weight categories could guide dietary recommendations.

While the study focuses on weight differences, it indirectly sheds light on the broader interactions between the host and its microbiota. The gut microbiota’s influence isn’t limited to weight but extends to the immune system, mental health, and more (78, 79). The identified differences in gut microbiota composition and function could have ripple effects across various physiological systems. These findings set the stage for further studies. Longitudinal analyses could explore whether these microbiota differences precede weight changes or result from them. Additionally, mechanistic studies could elucidate how the identified bacterial species influence the BMI z-score and whether modulating their populations can have tangible health benefits.

Our study offers preliminary evidence that full-length and V3–V4 16S rRNA sequencing methods identify significant correlations between certain gut microbiota and polysomnographic parameters in pediatric OSA. Both techniques demonstrated the capability to discern relationships between the composition of gut microbiota and OSA-related parameters. Notably, taxa such as Bifidobacterium, Escherichia, and Phocaeicola, which are either anaerobic or facultative anaerobic and typically reside in the human gastrointestinal tract, were examined. Previous research suggests that anaerobic conditions may promote the proliferation of Bifidobacterium and Escherichia (80, 81), corroborating our findings of their positive correlations with OSA severity indicators. Conversely, Phocaeicola, known for its distinctive oxygen tolerance (82), demonstrates an inverse correlation with the AI, supporting our observations of its relationship with OSA severity. Hence, NGS methods are valuable in elucidating the contribution of specific gut microbes to OSA pathogenesis.

Distinctly, the full-length sequencing provided a more comprehensive insight into potential OSA-associated metabolic pathways compared to the V3–V4 technique. The full-length approach has been instrumental in identifying correlations between microbial metabolism in diverse environments and AI, as well as in discovering and characterizing new microbial lineages that significantly enrich our understanding of microbial diversity (83). The associations of propanoate metabolism and the metabolism of glycine, serine, and threonine with OSA were also observed (84, 85). The impact of peptidoglycan biosynthesis, vital for bacterial structural integrity and stress tolerance (86), in relation to OSA and hypoxic conditions is novel and merits further exploration. These observed differences likely reflect the distinct genomic regions targeted by each sequencing method, indicating that each technique may provide a unique perspective on microbiome complexity. Furthermore, full-length 16S rRNA sequencing could offer a more nuanced understanding of the pathophysiological interplay between OSA and the gut microbiota.

The comparison underscores the importance of method selection in microbiome studies. There are some study limitations. First, the choice of sequencing method and analytic pipeline might lead researchers to different conclusions or emphasize certain findings over others (67, 87). It’s a reminder that while partial-length (such as V3–V4) sequencing is more common and might be more cost-effective, there might be trade-offs in terms of the depth and breadth of data acquired (88). Additionally, expanding the scope to include larger cohorts and other disease-associated microbiomes could enrich the preliminary results presented herein. The observed differences in microbial compositions between the two groups, as revealed by the full-length sequencing, can pave the way for longitudinal studies. Second, although certain bacterial species were more predominant in children with OWO than in children with NW and might play roles in metabolism, appetite regulation, or other physiological processes linked to weight, the cause-effect relationships between gut bacteria and weight categories need further investigations. These studies can assess if the microbial differences are a cause or consequence of being overweight or obese and if modulating the microbiota can have therapeutic benefits. Targeted interventions, such as probiotics or dietary changes, could potentially be designed to modulate the gut microbiota in favor of weight regulation.

The full-length 16S rRNA sequencing offers a promising avenue for microbial research, particularly when high taxonomic resolution is required. However, the choice between it and the V3–V4 region should be context-driven, factoring in the study’s goals, available resources, and the potential trade-offs of each method. Furthermore, the observed advantages of full-length sequencing in differentiating gut microbiota and predictive functions between children with OWO and children with NW, coupled with its ability to correlate these functions with OSA-related indicators, highlight its potential in advancing our understanding of the complex relationship between the gut microbiome and weight. It emphasizes the need for methodological rigor in microbiome meta-analysis and suggests a rich area for future investigations. The associations drawn between gut microbiota composition, metabolic pathways, and BMI z-score offer a promising avenue to understand weight management and overall health better.