Ethical approval was obtained from the COMIRB prior to the start of the study. Informed consent was obtained from study participants, including parents of children enrolled in the study. The IACUC of the Veterans Affairs Medical Center, San Diego, California granted approval for mouse studies.

Clinical data were obtained from 91 pediatric patients undergoing surgery for OM, with information on age, sex, self-reported ethnicity, family history, breastfeeding history, history of exposure to smoking, OM diagnoses and surgical technique ( Table 1 ). We also had clinical information and samples from 15 adult patients with OM, but these samples were removed from further analyses because of marked differences in expression and microbiota profiles due to age ( Figure 1 ). DNA samples were collected from the 91 pediatric patients with OM using the Oragene-DNA OGR-500 or OGR-575 kits (DNA Genotek, Ottawa, Ontario, Canada).

Saliva samples were also collected from pediatric patients with OM using Oragene-RNA RE-100 kits and sufficient RNA was isolated from 30 samples using the manufacturer’s protocol ( Figure 1 ). A total of 296 microbial samples were obtained from the ME (n=171) and NP (n=125) of 86 individuals, including 74 ME swabs, 86 ME aspirates, and seven ME mucosal tissue samples. Four ME cholesteatoma/granuloma tissue samples and 125 NP swabs were also collected ( Figure 1 ). Microbial DNA was isolated from 217 (73%) samples using the Epicentre Masterpure Complete DNA Purification Kit (Lucigen, Middleton, WI, USA); the rest of the samples from which no microbial DNA was isolated were excluded from further study ( Figure 1 ).

A variant in RASIP1 c.1801C>T (p.Arg601Cys; rs2287922) is in moderate LD with the FUT2 variant c.461G>A (r^{2 =} 0.82) and with the rs681343 variant (r^{2 =} 0.65) that was associated with childhood ear infections (Pickrell et al., 2016; Buniello et al., 2019). Sanger sequencing was performed for the FUT2 NM_000511.6:c.461G>A and RASIP1 NM_017805.3:c.1801C>T variants using DNA from saliva samples of pediatric patients with OM. Both variants were in Hardy-Weinberg equilibrium within the entire cohort and in each cohort used for RNA-seq and microbiota analyses ( Table 1 ).

Thirty salivary RNA samples (median RIN=7.1) were submitted for RNA-sequencing at the University of Colorado Denver Genomics and Microarray Core, as previously described (Larson et al., 2019). In summary, RNA samples were processed using the Nugen Trio RNA-Seq Kit (Tecan, Redwood City, CA, USA). Sequencing was performed on an Illumina HiSeq 4000 with an average of 31 million reads per sample. Reads were trimmed using the FASTX-Toolkit v0.0.13 and aligned using STAR v2.5.3a (Dobin et al., 2013). Principal components analysis (PCA) was performed on this dataset and one outlier sample was removed from further analyses due to not clustering with other samples ( Supplementary Figure 1 ). Transcript counts were summarized at the gene level and analyses included genes with an average read count >3. DE analysis was performed on 28 samples ( Figure 1 ) according to carriage of the FUT2 c.461G>A variant using the DESeq2 package in R (Love et al., 2014), with correction for age, sex and batch effects ( Supplementary Figure 1 ). Results were considered significant for genes with log₂-transformed fold change > ± 2 and false discovery rate (FDR)-adjusted p-value <0.05 using the Benjamini-Hochberg method.

FUT2, RASIP1 and DE genes were used as input in NetworkAnalyst for construction of a protein-protein interaction network using the IMEx interactome database (Xia et al., 2014; Xia et al., 2015). Pathway enrichment analysis was performed on the resulting network using the KEGG and PANTHER GO-slim BP databases in NetworkAnalyst (Kanehisa and Goto, 2000; Kanehisa, 2019; Mi et al., 2019; Kanehisa et al., 2021). Pathways with an FDR-adjusted p<0.05 were deemed significantly enriched.

A total of 171 ME and 125 NP samples were obtained from 86 Coloradan pediatric patients with OM and submitted for 16S rRNA sequencing. Microbial DNA isolation was performed using the Epicentre MasterPure™ Kit. In order to test for contaminating bacterial DNA in reagents or plastics, every batch of samples that was submitted for 16S rRNA gene PCR and sequencing included ≥3 negative process controls. Bacterial profiles were determined by broad-range PCR amplification and sequence analysis of the 16S rRNA gene V1V2 regions, as previously described (Santos-Cortez et al., 2018; Frank et al., 2021). Illumina paired-end sequencing was performed on MiSeq using the 600 cycle version 3 kit. Assembled and quality-filtered sequences were aligned and classified with SINA (1.3.0-r23838) using the 418,497 bacterial sequences in Silva 115NR99 (Pruesse et al., 2012; Quast et al., 2013). Operational taxonomic units (OTUs) were produced by clustering sequences with identical taxonomic assignments (median: 115,176 sequences/sample; interquartile range: 46,274.5 – 170,300.0). Goods coverage scores were ≥99.7% for all samples, indicating adequate depth of sequence coverage for all samples. Of the 296 microbial samples submitted for sequencing, 79 did not pass quality control (DNA concentration ≥10 ng/ul; 2500 reads after sequencing; Figure 1 ). Because it was not possible to determine whether the lack of microbial DNA is due to a relatively sterile ME or from a sample collection issue, these 79 samples were excluded. Bacterial alpha-diversity indices (richness, diversity, and evenness; Robertson et al., 2013) were tested for association with carriage of each of the FUT2 c.461G>A or RASIP1 c.1801C>T variants independently via Wilcoxon test and adjusted for ethnicity (Robertson et al., 2013). Associations of individual OTUs with FUT2 c.461G>A and RASIP1 c.1801C>T variants were assessed using linear regression with sample batch as a covariate. To minimize multiple-comparisons, only taxa with a prevalence >10% and relative abundance >1% were included in the analysis. Beta-diversity was determined via PERMANOVA using the Morisita-Horn dissimilarity index and adjusted for age, sex and batch effects. R software was used for data analyses and figure generation.

All animal experiments were performed according to the recommendations of the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health and carried out in strict accordance with an approved Institutional Animal Care and Use Committee (IACUC) protocol (A13-022) of the Veteran Affairs Medical Center (San Diego, CA). All animal experiments employed the best efforts for minimizing animal suffering under general anesthesia according to the NIH guidelines.

For gene array studies, wild-type (WT) C57/WB F1 hybrid mice were purchased from the Jackson Laboratory (Bar Harbor, ME USA). NTHi strain 3655 (non-typeable, biotype II, originally isolated from the ME of a child with OM in St Louis, MO USA) was cultured in defined liquid media (Coleman et al., 2003). To induce ME infection, mice were deeply anesthetized with an intraperitoneal injection of rodent cocktail (13.3 mg/ml ketamine hydrochloride, 1.3 mg/ml xylazine, 0.25 mg/ml acepromazine; at 0.1-0.2 ml per 25-30 g body weight of the mouse). The bullae were bilaterally exposed through soft tissue dissection via a ventral approach. A hole was made in the bulla with a 23 gauge needle, allowing approximately 5 μl of NTHi inoculum (~5x10⁴ CFU/mL) to be injected using a Hamilton syringe with a 30-gauge needle. After the injection of NTHi inoculum, the tympanic membranes were visually inspected and confirmed to be intact. The incision was then stapled and the mice were given normal saline and analgesics (buprenorphine at 0.05mg/Kg) subcutaneously while recovering on the heated mat. Following recovery from anesthesia the mice appeared healthy, with a clinical activity index ≤ 3 throughout the duration of OM.

Gene array data were generated as previously described (Hernandez et al., 2015). In summary, forty mice per time point were inoculated bilaterally with NTHi. Mucosal tissue and exudate were harvested from 20 mice at each of the following intervals – 0 hours (0h, no treatment), 3h, 6h, 1 day (1d), 2d, 3d, 5d and 7d after inoculation – then pooled. The tissue was homogenized in TRIzol (Life Technologies, Carlsbad, CA) and total RNA extracted, reverse transcribed and re-transcribed in vitro to generate biotinylated cRNA probes that were hybridized to 2 Affymetrix MU430 2.0 microarrays. Hybridization intensity data were median-normalized and differences in gene transcript expression levels evaluated using variance-modeled posterior inference (VAMPIRE) (Hsiao et al., 2005). Bonferroni multiple testing correction (αBonf < 0.05) was applied to identify only those genes with the most robust changes. The data were duplicated at each time point to obtain a second, independent biological replicate. Thus each data point represents 2 separate samples consisting of 20 mice each, and 4 Affymetrix arrays. A total of 3,605 genes, approximately 14.4% of the mouse genome, defined the signature of acute, NTHi-induced OM across time. Hybridization of RNA to specific gene probes was assessed at individual time points by comparison to uninfected MEs, after Bonferroni correction for multiple tests, using Genespring GX 7.3 (Agilent Technologies, Santa Clara, CA).

For single-cell RNASeq, the same ME inoculation protocol was followed, except that C57BL/6J mice (Jackson Labs) were employed. Single-cell samples for RNA-sequencing were generated from the entire contents of the mouse ME (Ryan et al., 2020). For each of three independent samples, tissue was harvested from both ears of six young adult C57-BL6 mice 6 hours after inoculation of the ME with NTHi. Single-cell libraries were generated using the 10X Genomics (Pleasanton, CA, USA) Chromium Single Cell 3’ Reagent Kit V2. cDNA synthesis, barcoding, and library preparation were then carried out on a 10X Genomics Chromium Controller according to the manufacturers’ instructions. After validating quality of cDNA library, sequencing was performed on an Illumina HiSeq 2500 (Illumina, San Diego, CA USA). Reads were demultiplexed and aligned to the murine reference genome (mm10 with annotations from Ensembl, release 84). 10X Genomics Cellranger aggr and Seurat were used to generate PCA clustering (Satija et al., 2015). The expression of well-recognized marker genes identified 24 distinct cell types (Ryan et al., 2020). Linearized relative expression levels of each gene examined in this study were log-transformed from single-cell mRNA copy numbers, normalized, and scaled for each cell type. Data were visualized in 10X Genomics cLoupe, with UMI numbers expressed colorimetrically for each cell.

Samples were collected from 91 pediatric patients with OM with ages ranging from 8.7 months to 14.9 years old (median 2.0 years; Table 1 ). Of these, 86 had sufficient DNA sample for Sanger sequencing of FUT2 c.461G>A ( Figure 1 ) and 83.3% are homozygous or heterozygous for the FUT2 variant. Carriage of the FUT2 variant was not associated with age, sex, ethnicity or OM diagnosis among children with OM ( Table 1 ). In the entire cohort and in each subset analyses, males were predominant (≥81%), which is a known phenomenon for OM (Paradise et al., 1997).

RNA-seq data from 28 pediatric patients (0.8 to 14.8 years old; Table 1 ) passed QC and were available for analysis according to FUT2 genotype. DE analysis was performed using FUT2 c.461G>A variant carriage as the classifier (5 wildtype and 23 variant carriers) and with adjustment for age, sex and batch effects ( Table 1 and Supplementary Figure 1 ). Five DE genes were significant, namely: FN1 (log-fold change = -3.7, FDR-adj-p=0.006); KMT2D/MLL2 (log-fold change = -3.8, FDR-adj-p=0.04); MUC16 (log-fold change = -4.3, FDR-adj-p=0.04); MTAP (log-fold change = +5.4, FDR-adj-p=0.006); and NBPF20 (log-fold change = -3.5, FDR-adj-p=0.04). In carriers of the FUT2 c.461G>A variant, FN1, KMT2D/MLL2, MUC16 and NBPF20 were downregulated whereas MTAP was upregulated ( Figure 2 ).

To further investigate how FUT2, FN1, KMT2D/MLL2, MUC16, MTAP, NBPF20 and RASIP1 are related, these genes were used as input for network analysis. RASIP1, FN1, KMT2D/MLL2 and MTAP were connected in a single protein-protein interaction network ( Figure 3A ). Pathway enrichment analysis of this network revealed 27 significant pathways in KEGG and 21 significant processes in PANTHER GO-slim BP, many of which overlap ( Figures 3B, C and Table 2 ). Among these are processes pertaining to viral and bacterial infection, cell cycle regulation, apoptosis, and endocytosis ( Table 2 ).

To further understand the role and interactions between FUT2, RASIP1 and DE genes, expression of orthologs Fut2, Fn1, Kmt2d, Muc16, Mtap and Rasip1 were measured by gene array in ME of wildtype mice at multiple time points (from 3 hours to 7 days) post-infection with NTHi ( Figure 4 and Tables 3 , 4 ). The NBPF gene family results from segmental duplication events in primate, thus an ortholog for NBPF20 is not present in mice (Vandepoele et al., 2005). Expression of Fut2, Rasip1 and Mtap were significantly increased after inoculation, with Fut2 and Mtap peaking around one day post-inoculation, and Rasip1 and Fn1 at 3 hours post-inoculation ( Figure 4 ). Additionally, expression of Muc16 was significantly decreased one day post-inoculation. For Fn1, Mtap and Muc16, DE was sustained through days 2-7 post-inoculation, including when OM is supposedly in recovery phase (Hernandez et al., 2015). Kmt2d showed no significant changes in ME expression at any point during the 7 days when compared to control mice ( Figure 4 ).

Single-cell RNA-sequence (scRNA-Seq) data were derived from the MEs of NTHi-infected mice six hours after inoculation ( Figure 5 and Table 5 ). In uninfected ME (time point 0h), Fut2 was expressed primarily in ciliated epithelial cells (Hydin+). Muc16 was expressed in most epithelial cells except basal epithelial cells (Krt14+). Rasip1 was expressed in most endothelial cells (Egfl7+) and Fn1 mostly in stromal cells (Col1a2+) and melanocytes (Mlana+). Mtap and Kmt2d were modestly expressed in all ME cell types. Six hours after ME inoculation with NTHi, when overall expression data was strongest ( Figure 5 ), Fut2 had increased expression in non-ciliated epithelial cells (Krt18/19+) and Muc16 in all epithelial cell types. Rasip1 continued to be expressed in endothelial cells, but was also observed in polymorphonuclear cells (PMNs) and monocytes (Csf1r+). Fn1 increased expression in stromal cells and monocytes and some endothelial cells ( Figure 5 ). Mtap and Kmt2d remained moderately expressed in all ME cell types except infiltrating PMNs and red blood cells. Level of gene expression per cell peaked at 1 day, and then declined ( Table 5 ). Taken together, the mouse ME expression profiles for Fut2, Rasip1 and DE genes support the findings of DE genes in OM patients using RNA-seq data from saliva ( Table 4 ), and also the overall expression of these genes in other human mucosal tissues ( Table 6 ).

A total of 296 microbial samples were collected from the NP and ME of 86 children ( Figure 1 ). For microbiota analyses, samples were filtered for: (1) those with >2500 16S rRNA sequencing reads; (2) one ME and one NP sample per individual where bilateral samples were collected (if bilateral, right-sided sample was used); and (3) available genotypes for FUT2 and RASIP1 variants ( Figure 1 ). No differences were identified between right and left NP or ME samples from the same individuals in PCA and PERMANOVA analyses (data not shown). After filtering, 16S rRNA sequence data from 34 ME and 65 NP samples were analyzed according to carriage of the FUT2 c.461G>A (p.Trp154*) variant.

In the ME, based on carriage of the FUT2 variant, Chao1 which denotes bacterial richness was significant when all ethnic groups were included (p=0.03); however, all alpha-diversity indices were not significant when only individuals of European descent were included in analyses ( Supplementary Table 1 ). Overall microbiota composition (i.e., beta-diversity) did not differ significantly by FUT2 variant according to PERMANOVA analysis with adjustment for age, sex, or batch effects ( Figure 6A ). Additionally, the relative abundances of Haemophilus (nominal p=0.03) and Moraxella (nominal p=0.02) were increased with wildtype FUT2 genotype, whereas increased Propionibacterium (nominal p=0.04) and Anoxybacillus (nominal p=0.02) were associated with presence (homozygous or heterozygous genotypes combined) of the variant ( Figure 6C and Supplementary Table 2 ). Performing these analyses by genotype had no overall effect on results ( Supplementary Figure 2 ).

In the NP, there were also no significant differences in alpha- or beta-diversity ( Supplementary Table 1 and Figure 6B ). Similar to ME, Propionibacterium had increased relative abundance in the NP (nominal p=0.01) among carriers of the FUT2 variant. In addition, the relative abundances of Actinobacillus (nominal p=0.03), Selenomonas (nominal p=0.03) and Candidate Division TM7 (Saccharibacteria; nominal p=0.0002) were increased in wildtype individuals ( Figure 6D , Supplementary Table 2 ). When individual taxa were tested for association by genotype, no taxa were significant ( Supplementary Figure 2 ). Note however that these FUT2-microbiota associations were nominal and were non-significant after FDR correction, with the exception of Candidate Division TM7 in the NP (FDR-adjusted p=0.009).

Sanger sequencing of DNA samples confirmed that the RASIP1 c.1801C>T and FUT2 c.461G>A variants are in moderate LD in our cohort as the genotypes for 57 of 71 (80.3%) individuals were identical. In the ME, similar to findings with the FUT2 variant, an increased relative abundance of Haemophilus (nominal p=0.04) was associated with wildtype genotype whereas increased Propionibacterium (nominal p=0.04) was associated with the RASIP1 variant ( Figure 7C and Supplementary Table 3 ). When analyzed by genotype, Haemophilus remained nominally associated with wildtype ( Supplementary Figure 3 ). In the NP, increased abundance of Propionibacterium (nominal p=0.006), chloroplast (FDR-adjusted p=0.05), Escherichia-Shigella (nominal p=0.04) and Staphylococcus (nominal p=0.04) was associated with carriage of the RASIP1 variant, whereas increased abundance of Candidate Division SR1 (FDR-adjusted p=0.05), Candidate Division TM7 (FDR-adjusted p=0.05), and Actinobacillus (nominal p=0.01) was associated with wildtype genotype ( Figure 7D and Supplementary Table 3 ).