Ethical approval was obtained from the Institutional Review Board (IRB) of Hallym University Dongtan Sacred Heart Hospital (IRB no. 2019-05-007-002). All steps of the study were conducted in accordance with the principle outlined in the Declaration of Helsinki concerning biomedical studies involving human participants. Written informed consent was obtained from the parents of each child.

The study was carried out as a prospective study between May 2020 and February 2021. The OME group consisted of children who solely required ventilation tube insertion. The inclusion criteria for ventilation tube insertion were the presence of COME persisting for > 3 months from the onset or diagnosis date. The control group consisted of children receiving tonsillectomy and/or adenoidectomy due to obstructive symptoms such as snoring and mouth breathing. Patients who had ingested oral antibiotics within 4 weeks before sample collection were excluded from the study.

Adenoid samples were collected during surgery under general anesthesia at Hallym University Dongtan Sacred Heart Hospital. To prevent contamination of the adenoid samples, we observed the adenoid with direct vision using a 70-degree rigid telescope (KARL STORZ, Tuttlingen, Germany) when the ventilation tube insertion was performed. The adenoid size was evaluated using an endoscopic grading scale as per the study by Cassano et al. that delineated adenoid enlargement into four distinct categories: grade 1, representing adenoid occupying less than 25% of the choanal area; grade 2, where the adenoids occupied 25% to 50% of the choanal area; grade 3, indicating adenoids occupying 50% to 75% of the choanal area; and grade 4, denoting adenoids occupying 75% to 100% of the choanal area (Cassano et al., 2003). Swabs were collected through an intraoral approach using sterile cotton (482CE; COPAN, Brescia, Italy). Collected swabs were placed in a collection tube and the remainder of the stem was discarded. The tubes were then capped and immediately stored at -80 °C until DNA extraction. Fecal samples were collected using an OMNIgene gut tube (DNA Genotek, Ontario, Canada) according to the manufacturer’s instruction. Collected fecal samples were immediately stored at -80 °C until DNA extraction. We examined clinical data related to COME, including a history of previous ventilation tube surgery, allergy testing (multiple allergens simultaneous test; MAST), the phenotype of middle ear effusion (serous or mucoid), and the interval between the last antibiotic usage and sampling. If multiple class 1 adenoids tested positive for at least one antigen in the MAST, it was considered positive. The criteria for defining mucoid effusion include the presence of mucus strands upon gross examination and the absence of gravitational flow of the effusion (Val et al., 2018; Lin et al., 2003).

Metagenomic DNA was extracted from adenoid tissue and fecal samples using the RNeasy PowerMicrobiome Kit (Qiagen, Inc., Valencia, CA, USA). The extracted DNA was purified using the DNeasy PowerClean Pro Cleanup Kit (Qiagen) according to the manufacturer’s instructions. The DNA concentration was measured using the BioPhotometer D30 with a μCuvette G1.0 (Eppendorf, Hamburg, Germany). To eliminate potential contaminations during the experimental process of sequencing-based studies (Eisenhofer et al., 2019), negative controls were included at every step. A total of 15 negative samples were sequenced along with the samples. These comprised sampling tubes, DNA-free water added to the DNA extraction kit, DNA-free water added to the purification kit, and DNA-free water added to the sequencing library preparation kit.

The DNA extracted from adenoidal (n=33) and fecal samples (n=32) were used to prepare the WMS library using a Swift 2S™ Turbo DNA Library Kit with Adapter (Swift Biosciences, Ann Arbor, MI, USA) following the manufacturer’s instruction. Index polymerase chain reaction (PCR) was performed using a 2S Combinatorial Dual Indexing Kit (Swift Biosciences). The library size was verified using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). The concentration of the library was measured using Qubit™ dsDNA HS Assay Kits (Thermo Fisher Scientific, Waltham, MA, USA). The equimolar concentration of each library was calculated using a TaKaRa PCR Thermal Cycler Dice Real Time System III (TaKaRa Bio, Inc., Shiga, Japan) with the GenNext NGS Library Quantification Kit (Toyobo, Osaka, Japan). Libraries were pooled and sequenced using the Illumina NovaSeq system (250-bp paired ends).

Sequence reads obtained from the NovaSeq system were analyzed as described previously (Lee et al., 2018, 2022). Adapter removal and quality filtering were performed using Trimmomatic with default options (Bolger et al., 2014). Paired-end sequences were merged using PEAR v.0.9.11 (Zhang et al., 2014). Human gene features were removed using BBMap with a reference human genome (HG-19). Taxonomic profiles were obtained using MetaPhlAn v.3.0 (Segata et al., 2012), and functional profiles were obtained using HuMAnN v.3.0 (Beghini et al., 2021). The resultant UniRef90 IDs were converted to KO. A total of 610,003,360 reads were obtained from the WMS analysis.

Total bacterial amounts in collected samples were estimated using quantitative real-time PCR based on the 16S rRNA gene, as described previously (Lee et al., 2022). The DNA extracted from samples was amplified in a 25-μL reaction mixture containing 12.5 μL of TB Green Premix Ex Taq (Tli RNaseH Plus) (TaKaRa Bio, Inc.), 2 μM of each primer, and 1 μL of DNA template (a 10-fold dilution series of sample DNA) with the primer 340F (5′-TCC TAC GGG AGG CAG CAG-3′) and 518R (5′-ATT ACC GCG GCT GCT GG-3′) with a Thermal Cycle Dice Real-Time System III (TaKaRa Bio, Inc.). The amplification conditions were 95 °C for 30 s, followed by 40 cycles of denaturation at 95 °C for 5 s and annealing at 60 °C for 30 s. Each sample was measured in triplicates. We quantified the bacterial amount by comparing threshold cycle values to a standard curve generated from parallel reactions of serial dilutions (1 × 10¹- 1 × 10⁷) of the 16S rRNA gene from the Escherichia coli K12 w3110 strain. Regression coefficients (r²) for all standard curves were ≥ 0.99.

To mitigate potential contamination during the experimental process, decontamination was performed using the R package ‘decontam’ based on the sequences of the negative controls (Davis et al., 2018). Decontamination was performed with a threshold of 0.5 prevalence in the ‘isContaminant’ function. After decontamination, the analyzed features were selected based on a > 10% sample prevalence to increase the relevance of the results.

The Shannon diversity index was calculated using the R package ‘microbiome’. The principal coordinates analysis (PCoA) plot was obtained based on Bray-Curtis dissimilarity.

The influence of covariates, including age, sex, date of last antibiotic administration, duration of COME, adenoid size, allergies, and tube surgery, on the differences in the microbiome was analyzed using the ‘envfit’ function in the R package ‘vegan’ (v.2.5-7). The significance of covariates was calculated using permutational multivariate analysis of variance (PERMANOVA), and p-value < 0.05 were considered statistically significant.

Multivariate regression analysis was performed using multivariate association with linear models (MaAsLin2) (Mallick et al., 2021). The taxonomic and functional features were normalized by relative abundance and cumulative sum scaling (CSS) methods, respectively.

The important species for determining the microbiota age were selected through random forest regression with the ‘randomForest’ R package. The chronological age of the microbiota in children without COME was used as a training set to estimate the microbiota age (Galazzo et al., 2020). The N top discriminatory species were identified using the ‘rfcv’ function in the randomForest package. The increased node impurity index (IncNodePurity) to determine the importance of discriminatory species was calculated using the ‘importance’ function in the randomForest package. The relative abundance of important species was compared among groups and visualized in a heatmap using the ‘ggplot2’ R package. Spearman’s correlation was used to analyze the correlation between the estimated microbiota age and chronological age.

The indicator species for each group (control, serous-, and mucoid-COME) in children aged 6–12 years were determined using the ‘multipatt’ function of the ‘indicspecies’ R package with 10,000 permutations. Results with p < 0.05 and an indicator value > 0.5 were considered indicators for each group. The indicator species are presumed to be the main contributors to the observed microbiota differences across the groups.

The interspecies correlations in adenoid microbiota were analyzed using FastSpar with 1,000 bootstraps (Watts et al., 2019). The corresponding correlation network was visualized using the ‘qgraph’ R package. Significant correlations with p < 0.05 and a correlation value < |0.2| were demonstrated in the network.

Functional features in the adenoid microbiome were analyzed based on the KO category. The relative abundance values, represented as Reads Per Kilobase (RPK), for the KO categories and gene families, were normalized using the CSS method. The KO categories and gene families were compared between groups using MaAsLin2. Functional features between groups were depicted visually on volcano plots to highlight significant differences. The x-axis displays the fold-change of gene families (log₂ fold-change), and the y-axis demonstrates significance with -log₁₀ (p-value) calculated using MaAsLin2. Pearson’s correlation analysis was conducted to explore the relationship between indicators for mucoid COME and significantly different gene families (p < 0.001).

The variables of clinical data between groups were compared using a two-sided Fisher’s exact test. The differences in taxonomic and functional features between groups were determined using the Wilcoxon-rank sum test in the R software. Dunn’s multiple comparison tests were used to identify the differences among groups using the ‘dunn.test’ R package. The Benjamini-Hochberg false discovery rate multiple testing correction was utilized to adjust the p-value where applicable. The beta-diversity of taxonomic and functional features was visualized using the PCoA plots based on Bray-Curtis dissimilarity, and significance was determined using PERMANOVA (Adonis from the package vegan with 999 permutations). Results with p < 0.05 and q (adjusted p-value) < 0.05 were considered statistically significant.

The number of participants, age, and sex were not significantly different between the control and COME patient groups ( Table 1 ). We analyzed 80 WMS data (adenoidal sample = 33, fecal sample = 32, and negative control sample = 15) in this study. The Decontam pipeline was utilized to trim potential contaminant sequences, identified based on the sequences detected in 15 negative controls, encompassing the sampling to library preparation processes ( Supplementary Figure 1 ). Microbiota detected in negative samples significantly differed from those in adenoidal and fecal samples (p < 0.001). Further microbiome analyses were performed after removing potential contaminants.

The covariates significantly associated with the variation of microbiota were determined using the EnvFit model ( Supplementary Table 1 ). The variation in adenoid microbiota was observed to be associated with COME disease (r² = 0.097, p < 0.05). Age was a significant factor influencing the composition of the adenoid microbiota in the control group (r² = 0.415, p < 0.05). Conversely, in the COME group, the type of middle ear fluid (MEF), whether serous or mucoid, was identified to be a significant determinant of adenoid microbiota variation (r² = 0.278, p < 0.05).

The diversity of adenoid microbiota in COME with mucoid MEF was lower than those in the control and COME with serous MEF groups (p < 0.05; Figure 1A ). On the other hand, the bacterial amounts in the adenoid microbiota of COME with mucoid MEF were higher than those in the other groups. The bacterial diversity and amounts in the gut microbiota did not differ significantly among the groups ( Figure 1B ). The adenoid microbiota demonstrated significant differences among the groups based on Bray-Curtis dissimilarity in the PCoA plot (p < 0.01; Figure 1C ). This result was consistent with the influence of covariates as determined by the EnvFit model ( Supplementary Table 1 ). Notably, Actinobacteria and Proteobacteria were identified as significantly different phyla between the groups. However, no significant differences were observed in gut microbiota among the groups in the PCoA plot or in terms of phylum composition ( Figure 1D ).

Previous studies have reported the peak prevalence of OME between the ages of 2 and 5 years, coinciding with the period of rapid adenoid growth observed between the ages of 3 and 5 years (Dhooge et al., 2005; Papaioannou et al., 2013). Therefore, we analyzed the variation in microbiota according to age and compared the differences between the 2–5 years and 6–12 years age groups. In agreement with these previous studies, age was identified as a significant covariate in the adenoid microbiota of the control group. However, this age-dependent variation was not observed in the gut microbiota, but only in the adenoid microbiota of the control group ( Figures 1E, F ).

We analyzed the adenoid microbiota across different age groups. Random forest modeling after cross-fold validation identified eighteen species, including Porphyromonas somerae and Gemella sanguinis, as influential age-dependent species ( Figure 2A ). The relative abundances of these species according to age were perturbed in children with COME compared to those in controls. The predicted age of the adenoid microbiota was positively correlated with the chronological age in the control group (rho = 0.47, p < 0.05). In contrast, the adenoid microbiota in the COME group exhibited age-independent patterns ( Figure 2B ).

Our analysis revealed age-associated perturbations in the adenoid microbiota compositions among children with COME. The alteration in adenoid microbiota was analyzed according to the type of MEF within each age group. Although the diversity was significantly different between control and COME groups (p < 0.05), bacterial amounts and microbiota compositions did not differ among the groups according to MEF types in children aged 2–5 years (p > 0.05; Figure 2C ). In contrast, alterations in microbiota corresponding to MEF types of COME were evident in children aged 6–12 years ( Figure 2D ). The diversity was lowest in the mucoid COME group compared to that in the other groups (p < 0.01), whereas the bacterial amounts were highest in the mucoid COME group (p < 0.05). The composition of microbiota exhibited significant differences according to MEF types in the PCoA plot (p < 0.01). The relative abundances of Bacteroidetes and Firmicutes were higher in the control group than in the mucoid COME group, whereas the relative abundance of Proteobacteria was higher in the mucoid COME group than in the control group (p < 0.05). Actinobacteria was higher in the serous COME group than in the mucoid COME group (p < 0.05). These results suggest that the adenoid microbiota in the COME group is differentially affected by the age of the children. In children aged 6–12 years, the composition of the adenoid microbiota differed significantly among COME groups depending on MEF type. The gut microbiota did not significantly differ among the groups in both age categories of 2–5- and 6–12-year-old children ( Supplementary Figure 2 ).

We analyzed the indicator species for each group in children aged 6–12 years to study the shift in adenoid microbiota according to MEF types. Fusobacterium nucleatum, G. sanguinis, and Eubacterium sulci were identified as indicators in the control group ( Figure 3A ). The relative abundances of F. nucleatum and E. sulci decreased in the COME group compared to those in the control group, whereas the relative abundance of G. sanguinius was higher in the serous COME group than in the other groups (p < 0.05). Thirteen species, including Streptococcus salivarius, Veillonella parvula, and Neisseria mucosa, were identified as indicators for the serous COME group, and their abundances were higher in the serous COME group than in the other groups (p < 0.05). Streptococcus pneumoniae and H. influenzae were indicators for the mucoid COME group. Although indicator species were commonly detected in all groups, their abundances varied across the groups.

These compositional differences in the relative abundances of indicator species across the groups could be partly driven by interactions among the indicators within each group. Significant positive correlations were observed among species within each group, except for those in the mucoid COME group ( Figure 3B ). Positive correlations were observed among the indicators of the control group, with higher correlation values than those noted in the other groups. Complex interactions between indicators were detected within the serous COME group, whereas indicators for the mucoid COME groups were not correlated with each other. Additionally, indicators for the serous COME groups were correlated with indicators for the control and mucoid COME groups. However, no direct correlation was identified between indicators for the control and mucoid COME groups. In the serous COME group, N. mucosa and Cardiobacterium hominis interacted negatively with E. sulci and F. nucleatum in the control group, respectively. This interaction could be associated with the reduced proportion of E. sulci and F. nucleatum in the serous COME group. H. influenzae in the mucoid COME group interacted positively with Rothia dentocariosa in the serous group, however, the organism interacted negatively with Capnocytophaga leadbetteri. S. pneumoniae in the mucoid COME group interacted positively with Streptococcus sp. F0442 in the serous COME group.

The adenoid microbiome may play different roles in children with COME depending on the age group, as the microbiome demonstrates distinct alterations. To analyze the potential influence of the altered microbiome, functional features based on the KO were compared among groups within each age category. The functional features of the adenoid microbiome were not significantly different among the groups according to MEF types in the children aged 2–5 years (p > 0.05; Figure 4A ). Although seven gene families differed among the groups (p < 0.01), they did not significantly contribute to the functional differences of the adenoid microbiome according to MEF types. However, the functional features of the adenoid microbiome significantly differed among groups within 6–12-year-old age group (p < 0.01; Figure 4B ). A total of 120 gene families exhibited significant differences across the groups (fold changes > 2 and p < 0.01). Among these 120 significant features, 18 gene families had a p-value < 0.001. The number of features that differed between the control and mucoid COME groups (11 gene families with fold changes > 2 and p < 0.001) was higher than in the other comparisons (two gene families between control and serous COME groups; five gene families between serous and mucoid COME groups). These differences in functional features among the groups were consistent with previous microbiota results ( Figure 2 ).

We analyzed the correlated gene families with indicators for the mucoid COME groups to identify potential functions of the altered adenoid microbiome. We identified nine gene families that were significantly correlated with S. pneumoniae and H. influenzae (p < 0.05; Figure 5A ). Four gene families (K21602, K08092, K22025, and K02037) demonstrated a positive correlation with two indicators, and their abundances increased in the mucoid COME group (p < 0.05; Figure 5B ). In contrast, five gene families (K08986, K09787, K05601, K00797, and K03665) exhibited a negative correlation with two indicators, and their abundances decreased in the mucoid COME group (p < 0.01). The D-erythronate catabolism pathway may be activated in the mucoid COME group through D-erythronate 2-dehydrogenase (denD; Figure 5C ). The level of spermidine synthase (speE) was lower in the mucoid COME group compared to its level in the other groups.

In comparison to the control and serous COME groups, acylphosphatase (K01512; acyP) was identified to be a significantly different gene family in the adenoid microbiome of the mucoid COME group (p < 0.001; Figure 4B ). The normalized abundance of acyP was lowest in the mucoid COME group (p < 0.001; Figure 5D ), which may be related to the reduced production of acetate. This was associated with the reduced proportion of E. sulci, an indicator for the control group, in the adenoid microbiome of the mucoid COME group.