Species- and strain-specific microbial modulation of interferon, innate immunity, and epithelial barrier in 2D air–liquid interface respiratory epithelial cultures

Abstract

Background

The microbiome regulates the respiratory epithelium’s immunomodulatory functions. To explore how the microbiome’s biodiversity affects microbe-epithelial interactions, we screened 58 phylogenetically diverse microbes for their transcriptomic effect on human primary bronchial air–liquid interface (ALI) cell cultures.

Results

We found distinct species- and strain-level differences in host innate immunity and epithelial barrier response. Strikingly, we found that host interferon, an antiviral response, was one of the most variable host processes. This variability was not driven by microbial phylogenetic diversity, bioburden, nor by the microbe’s ability to stimulate other innate immunity pathways.

Conclusions

Microbial colonization differentially stimulates host gene expression with variations observed across phylogenetically diverse microbes and across different strains of the same species. Our study provides a foundation for understanding how the respiratory microbiome’s biodiversity affects epithelial, and particularly antiviral, innate immunity.

Background

The respiratory epithelium is the first line of defense against potentially damaging environmental stimuli, such as allergens, toxins, and microbial pathogens. In addition to being a highly specialized physical barrier, regulating what can enter systemic circulation and removing foreign particles, it can directly contribute to innate immunity [1, 2]. For example, the respiratory epithelium produces antimicrobial compounds (e.g., lysozymes [3], defensins [4, 5], and cathelicidins [6]) and numerous cytokines, including interferons which have direct antiviral functions [7] and other cytokines that recruit immune cells [8,9,10]. Therefore, understanding the factors that modulate homeostasis and respiratory epithelial response to diverse stimuli is critical to developing new therapies that address infectious or chronic respiratory diseases.

A major modulator of the epithelium and its immunomodulatory functions is the microbiome, the diverse bacteria, fungi, and viruses that colonize the human body. Microbial colonization in the respiratory tract varies significantly between different regions of the respiratory tract due to differing environmental conditions. For example, the nares more closely resembles the skin in terms of species’ prevalence, with Staphylococcus (S.), Corynebacterium, and Cutibacterium as some of the predominant genera [11, 12]. The microbes colonizing the oropharynx resemble oral microbes, with genera like Streptococcus (Str.), Neisseria, Prevotella, and Rothia [13, 14]. Finally, the lung microbiome consists of predominantly anaerobic genera such as Prevotella, Streptococcus, and Veillonella [13]. The phylogenetic biodiversity found in the respiratory tract is significant because the effects of microbial colonization can be species- and strain-dependent, such that different microbial isolates of the same genus or species can have divergent effects on the host [15,16,17].

However, we currently have a poor understanding of how microbial phylogenetic diversity affects host-microbiome interactions. Studies to date [15,16,17,18,19,20,21] have only examined a handful of microbial isolates, which offers a limited recapitulation of the biodiversity seen in the respiratory tract and limits our understanding of the preventative or synergistic role of the microbiome in infection response as well as potential development of microbiome-based therapies that could enhance immune response or barrier integrity. One barrier to systematically investigating microbe-epithelial interactions is that few models can accommodate microbial diversity. For example, murine models, while possessing a full repertoire of innate and adaptive immune elements, are relatively low-throughput and possess endogenous microbiota that can differ significantly from the species found in human microbiomes [22]. While enabling throughput, in vitro human cell culture screens using epithelial monolayers cannot replicate the cellular complexity as well as the inherently spatially constrained interactions with the microbiome that is observed in the native respiratory epithelium.

Here, we leverage 2D air–liquid interface cell cultures (ALI) to provide a more physiologically relevant model with sufficient throughput to allow us to test a range of diverse microbiota and more thoroughly interrogate how our phylogenetically diverse cohort of microbes differentially modulates host gene expression. ALIs are a flat, multi-cell type, multi-layered epithelium, derived from epithelial progenitor cells, with the apical surface exposed to the air. They are morphologically [23, 24] and transcriptionally [25] more similar to in vivo epithelium than monolayer cultures and, importantly, provide the necessary throughput to examine the effects of colonization of a diversity of microbiota. We colonized ALIs with 58 phylogenetically diverse, patient-derived respiratory microbes with the overarching goal of reconstructing how transcriptional activation of host processes, especially antiviral innate immunity, differed based on microbial phylogeny. Strikingly, we identified notable species- and strain-level differences in epithelial barrier integrity and host innate immunity, particularly interferon induction, a canonically antiviral immune response. Taken together, our study provides the largest dataset, thus far, investigating how microbial phylogenetic diversity affects respiratory microbiome modulation of host antiviral innate immunity, antibacterial innate immunity, and barrier responses in a 2D ALI respiratory epithelial model.

Results

Microbiome cultivation

As part of our pipeline to cultivate respiratory microbiota for screening, we collected samples from 47 patients (ages ranging from 46–83 years) who were recruited as part of a broader study investigating immunological changes in lung cancer. The patients were sampled at the tongue dorsum and nares as non-invasive approximations of the respiratory tract [14, 26] (Fig. 1A, Additional file 1: Table S1). While we emphasize that our goal in this study was not to categorically identify differences in lung cancer disease state, we characterized the metagenomes of these patient samples (Additional file 1: Table S2). In the nares, a few genera (i.e., Cutibacterium, Corynebacterium, and Staphylococcus) comprised the majority of identified microbes, which was concordant with other characterizations of healthy adult nasal microbiomes [11, 12] (Fig. 1B). Tongue dorsum samples were more diverse, including Rothia mucilaginosa, Prevotella spp., Streptococcus spp., Veillonella spp., Neisseria subflava, and Haemophilus parainfluenzae, which have previously been reported as highly abundant at this oral site [11, 27] (Fig. 1C).

[IMAGE OMITTED: SEE PDF]

Following metagenomic analyses, we cultured 1154 microbes as previously described [28] from the samples of 11 randomly selected patients (Fig. 1A). Our goal was to isolate a diversity of microbial species for our ALI host-interaction screen. As anticipated, common microbiota such as S. epidermidis and Streptococcus spp. were isolated with greater frequency, allowing potential investigation of strain-specific effects (Fig. 1D, Additional file 1: Table S3). Most species were cultured fewer than 20 times in both the oral and nasal samples, demonstrating potential patient-specific effects. To note how representative cultivars were of the predicted metagenome, we compared the number of uniquely identified microbes using metagenomics versus those identified using culturomics at both a genus and species level (Fig. 1E). Unsurprisingly, we saw more unique genera/species identified using metagenomics than using culturomics at both sites. While metagenomics identified 76–91% of the total genera/species detected by either method, culturomics only identified 28–53%, with oral samples tending to recover fewer genera/species than the nasal samples.

Microbial selection and colonization of ALIs

We sought to rationally subset the 1154 isolated respiratory microbes for follow-up investigations in ALIs, and at the same time gain insights into the diversity of potential phenotypes of interest generated by different microbial phylogenies. Because secreted microbial metabolites are immunogenic and microbial immunogenicity can be species- and strain-specific, we examined the inflammatory potential of microbial metabolites on lung cell cultures as a high throughput means of screening the isolated microbes. We collected conditioned bacterial supernatant from each isolated microbe, applied the supernatants (20% v/v) to A549 lung epithelial cells overnight, and measured the concentration of six human inflammatory cytokines (IL-8, IL-6, IL-1, IL-12, IL-10, and TNF) using cytokine bead array (CBA). Interestingly, bacterial conditioned supernatant only caused significant changes to the secretion of IL-8, a neutrophil chemoattractant [29] (Additional file 1: Table S4). Notably, modulation of IL-8 was not coupled to phylogeny (Additional file 2: Fig. S1A). At the genus level, the median log2 fold change for nearly every genus indicated decreased IL-8 secretion but within-genus and even to the strain level, there was notable isolate specificity with a large log2 fold change (FC) range in IL-8 secretion (Additional file 2: Fig. S1B).

We ultimately selected 56 microbes for colonization of ALIs, prioritizing based on genetic diversity (different species), maximizing the effect on IL-8 secretion (increased or decreased secretion) (Fig. 2A), and representation of the predominant genera present in the respiratory tract (e.g., Streptococcus, which is one of the main genera in the lower respiratory tract and Staphylococcus, one of the main genera in the upper respiratory tract). We also included a GFP-tagged S. epidermidis strain Tü3298 [30, 31] for a visualization and colonization control, a skin-derived S. epidermidis strain 3F3 [32, 33] as an additional colonization control, and water as a vehicle negative control. 107 CFUs of each microbe (colony-forming units, a measure of live bacteria) or vehicle were applied to the apical surface of mature tracheobronchial ALIs for 18 or 48 h (Fig. 2B, Additional file 2: Fig. S2-S3) to approximate spatial aspects of physiological conditions. At each time point, total eukaryotic RNA was extracted from the tissue and processed for RNA-seq (Additional file 1: Table S5), and basal media was collected for IL-8 quantification (Additional file 1: Table S6), as we did not find significant amounts of IL-8 to be secreted apically with bacterial colonization. In addition, IL-8 was the only cytokine with notable changes in our initial screen. At each timepoint, we also washed the apical surface to quantify CFUs to approximate bacterial growth/load (Additional file 2: Fig. S4, Additional file 1: Table S7). CFUs were also quantified from the inoculum for comparison (timepoint 0 h). We observed several growth patterns: (1) a decline in CFUs at 18 h followed by increasing CFUs by 48 h, (2) continuous bacterial growth at both timepoints, (3) a continuous decrease in CFUs at both timepoints, or (4) CFUs at a steady state, and these patterns were genus, species, and strain-specific. For example, most Staphylococcus spp. experienced an initial decline followed by growth equilibrating at ~ 107 CFUs by 48 h, though the magnitude of the decline and ultimate equilibration differed by species (e.g., S. haemolyticus and S. lugdunensis vs. S. aureus). Similarly, Streptococccus spp. equilibrated around ~ 105 CFUs with varying trajectories, with some isolates (Str. parasanguinis, Str. salivarius, and Str. vestibularis) decreasing through 48 h. We note that isolates may differ in their adhesive properties, thus this analysis is an approximation of microbial growth and load.

[IMAGE OMITTED: SEE PDF]

Microbial modulation of host gene expression is not phylogenetically conserved

First, we sought to determine if bacterial growth/load was the primary driver of the transcriptional changes observed, as it is conceivable that excessive growth on the ALIs would trigger a broad epithelial response. We used a linear mixed effects model (MaAsLin2 [34]) to determine if CFUs at the time of harvest were a confounder for differential gene expression. We observed that only 68 genes (of 42,609 genes expressed in our model) were driven by CFUs, none of which were in our pathways of interest for subsequent analyses.

Next, we examined global differences in the effects of microbial colonization on host gene expression as a function of time. A qualitative comparison of 18 and 48 h timepoints with principal components analysis (PCA) showed that most of the microbial treatments at both 18 and 48 h clustered together and that the outliers were not phylogenetically similar (Fig. 2C). This suggested that species-specific epithelial responses to colonization could be primarily compartmentalized to certain pathways such as innate immune response or barrier integrity. The number of upregulated differentially expressed genes (DEGs) at each time point was similar, with substantial overlap between timepoints. However, there were markedly more downregulated genes at 48 h (Fig. 2D). Examining the top 20 gene sets with the highest gene ratio (percent of genes enriched in the gene set) from either time point when all microbial treatments were combined (i.e., the most universally enriched across treatments), we observed an overwhelming number of enriched gene sets were from the 48 h time point, with low gene ratios at 18 h (Fig. 2E, Additional file 1: Table S8-S9). Interestingly, we saw enrichment in several antiviral gene sets including Type 1 interferon (IFN-I), which is a canonical antiviral response, and few studies have investigated microbiome stimulation of epithelial interferon. As these results suggested that transcriptional differences accumulate and are most discernible at 48 h, we focused on data from the 48 h timepoint for the remainder of analyses.

To investigate if the effects on epithelial gene expression at 48 h were phylogenetically driven, we compared the number of upregulated and downregulated DEGs between different microbial treatments (Fig. 2F lost the phylum/genus count), categorizing largely by phyla, with genus-level designations Staphylococcus and Streptococcus spp. given their predominance and Candida given there was only one isolate. Interestingly, the number of DEGs did not correlate with the group size as one might anticipate; for example, though Streptococcus had the largest number of treatments, its effect on the epithelium was moderate, as measured by the total number of upregulated DEGs. Strikingly, Streptococcus, despite having a comparable number of unique upregulated DEGs, shared few upregulated DEGs with Staphylococcus and other Firmicutes. Staphylococcus shared more upregulated DEGs with Proteobacteria than it did with Streptococcus. Regardless of how many groups were compared, Actinobacteria shared few upregulated DEGs and Candida even fewer. The strong similarity between Staphylococcus' upregulated DEGs and those of other Firmicutes is unsurprising, given the closer phylogeny. However, Streptococcus’ and Actinobacteria’s few upregulated DEGs shared with other groups suggests that microbial induction of host epithelial genes is not highly phylogenetically conserved and can vary in both effect and quantity. Contrastingly, with the exception of the fungus Candida, which only included one treatment and is phylogenetically most removed from the bacterial isolates, there were more shared downregulated DEGs, indicating that aspects of microbial inhibition of host epithelial genes are likely conserved.

We then examined if phylogenetically more similar microbes induced similar gene expression changes in the 25 gene sets that were enriched in the most microbial treatments. Most of the enriched gene sets, irrespective of phylogeny, were related to general inflammation, such as neutrophil migration/chemotaxis and cytokine binding/activity. Other gene sets included antiviral response (response to virus and IFN-I). Notably, these enrichments crossed phylogenetic lines—i.e., almost all gene sets were upregulated by at least one member of each phylogenetic group—but surprisingly, no gene sets were upregulated by all members of a phylogenetic group. This further supports the notion that microbial induction of host gene expression is not phylogenetically conserved.

Microbial induction of interferon is not phylogenetically determined

We then further investigated the surprising enrichment in IFN-I gene sets, because in the respiratory epithelium IFN-I (as well as Type 3 interferon (IFN-III) which is also expressed in the respiratory epithelium) is canonically an antiviral immune response that results in the activation of interferon-stimulated genes (ISGs). To quantitate bacterial stimulation of interferon, we manually curated a list of 61 ISGs from the literature [35] representing IFN-I and IFN-III signaling, which have almost identical signaling cascades (Additional file 1: Table S10). Interferon stimulation at 48 h was strikingly dichotomous (Fig. 3A) with approximately half of the microbes strongly stimulating interferon (i.e., strongly upregulating ISGs), while the other half had a minimal effect. For each microbial treatment, we calculated the median log2FC for the ISGs to create an ISG score (Fig. 3A, Additional file 1: Table S11), which we confirmed using previously published 6- [36] and 30-gene [37] ISG scores (Additional file 2: Fig. S5). Based on the hierarchical clustering and ISG score, we then classified the microbes as ISG non-stimulators versus simulators (Fig. 3B, P= 9.9e − 11, two-sided Mann–Whitney U test with Bonferroni correction). We compared the effect of microbial colonization on interferon genes to determine whether the microbes were inducing IFN-I or IFN-III; however, we saw no significant differences (data not shown).

[IMAGE OMITTED: SEE PDF]

Strikingly, the ability of a microbe to stimulate ISG was not predicted by phylogeny, even at the strain level. Nearly every genus included both non-stimulators and stimulators with a wide range of stimulatory potential (Fig. 3C). Interestingly, we even observed strain-level differences, e.g., Rothia dentocariosa (Fig. 3D). Strains could moreover differ temporally in their ability to activate interferon; a Paenibacillus lactis strain and a Str. parasanguinis strain induced interferon at 18 h while the other strains of those species (and the vast majority of ISG stimulators) induced interferon only at 48 h (Additional file 2: Fig. S6). Finally, to confirm that bacterial growth/load did not solely explain the difference in interferon stimulation, we compared CFUs at 48 h between the ISG non-stimulator and stimulator microbial treatments and found no significance (Additional file 2: Fig. S7, P = 0.12, two-sided Mann–Whitney U test).

To experimentally validate interferon activation, we selected three microbial treatments (Str. vestibularis, L. rhamnosus, and GFP S. epidermidis Tü3298) and validated MX1 expression, a direct antiviral mediator protein known to be induced by interferon, by immunofluorescence. These three isolates represented the extremes of transcriptional effects on colonization; S. epidermidis elicited the most upregulation, Str. vestibularis had the most downregulated genes and L. rhamnosus had a surprisingly minimal effect on any genes despite robust growth on ALIs (Fig. 3E). MX1 expression was quantified via immunofluorescence. While none of the microbial treatments caused statistically significant changes, colonization with S. epidermidis distinctly increased expression of MX1 while minimal changes in MX1 expression were observed for the other two treatments (Fig. 3F, G). S. epidermidis’ and L. rhamnosus’ effect on MX1 expression matched their respective stimulator and non-stimulator classification. However, Str. vestibularis’ lack of MX1 expression was discordant with the RNA-seq results. One explanation is that while Str. vestibularis can induce gene expression of MX1, expression does not pass the threshold for protein expression. Taken together, the degree of species and strain specificity in microbial induction of interferon suggests that small variations or changes to the respiratory microbiome may have significant effects on antiviral immunity.

Microbial induction of interferon is independent of antibacterial innate immune response and epithelial barrier quality

Host innate immunity includes a wide array of antibacterial responses to resist genetically diverse bacteria. Our 58 microbes included a mix of canonical commensals and opportunistic pathogens, and thus we anticipated that some microbes would be more immunogenic while others may be immunosuppressive. To investigate if our microbial treatments differentially activated those antibacterial responses, we manually curated a list of 36 genes including non-interferon cytokines and antimicrobial peptides from MSigDB [38] with gene ontologies [39, 40] relating to cytokines and antimicrobial immune response, which we termed antibacterial innate immune response (Fig. 4A, Additional file 1: Table S10). An antibacterial score was generated based on median log2FC, identifying two classes of microbial induction: weak versus strong stimulator, again identified by hierarchical clustering and having significantly different antibacterial scores (Fig. 4B, Additional file 1: Table S11). Like interferon response, microbial phylogeny did not predict whether the microbe would be a strong or weak effector of antibacterial immunity (Fig. 4C). However, unlike interferon, induction was not dichotomous and was instead more continuous within and between microbial clades. To experimentally validate these broad categorizations, we measured the concentration of IL-8 (CXCL8) secreted into the basal media following microbial colonization for 18 and 48 h. Nearly all microbial treatments at least slightly increased secretion of IL-8 at 18 h and at 48 h (Fig. 4D), and secretion was correlated with transcriptional levels (Additional file 2: Fig. S8, ρ = 0.72, P < 2.2e − 16, Spearman’s correlation coefficient, Additional file 1: Table S6). Unlike interferon, when we compared CFUs at 48 h between strong and weak antibacterial inducers, we found a significant difference (Additional file 2: Fig. S9, P = 4.2e − 06, two-sided Mann–Whitney U test). While none of the individual antibacterial genes were predominantly driven by CFUs per our multivariate model, grouping these genes together increased the effect to be statistically significant. This was in striking contrast to interferon’s induction, which was independent of microbial load.

[IMAGE OMITTED: SEE PDF]

Epithelial barrier integrity affects the penetrance of microbes and microbial products and thus immunogenicity. Some components of the epithelial barrier, such as mucus secretion, are often modulated by microbial colonization [1] while other components, such as tight junctions, are common targets of microbial pathogens [41]. Based on this, we also examined epithelial barrier genes. We curated a gene set of keratin, mucin, and junction (inclusive of tight junctions, gap junctions, and adherens junctions) genes (Additional file 2: Fig. S10, Additional file 1: Table S10), similar to ISGs and antibacterial gene lists, from the literature [42, 43] and HGNG (HUGO Gene Nomenclature Committee) [44]. We experimentally validated a subset of the transcriptional changes observed in the epithelial barrier genes. We measured the transepithelial electrical resistance (TEER), a measurement of barrier permeability, of ALIs treated with GFP S. epidermidis Tü3298, Str. vestibularis, and L. rhamnosus at 18 h and 48 h (Fig. 4E, Additional file 1: Table S12). L. rhamnosus, which was an ISG and antibacterial non-stimulator, showed increased TEER and therefore decreased permeability, while S. epidermidis, an ISG stimulator and antibacterial strong stimulator, decreased TEER/increased permeability. Str. vestibularis, an ISG stimulator and weak antibacterial stimulator, did not have a significant effect on TEER. Altogether, these results aligned with each microbe’s degree of antimicrobial stimulation and not ISG stimulation. This suggests there may be a connection between microbial induction of antibacterial pathways and epithelial barrier genes.

To further generalize these potential links between epithelial barrier, antibacterial response, and ISG induction, we determined if the same microbial species were inducing interferon and antibacterial innate immunity, and therefore, were more generally immunogenic. Notably, there was little overlap in a microbe’s ability to induce interferon versus antibacterial innate immunity; some microbes induced both, one of the two, or neither, again with no discernable phylogenetic driver (Fig. 4F). In addition, ISG expression did not have a strong correlation with antibacterial, mucin, or keratin gene expression. However, ISG expression did correlate with junction genes. Furthermore, antimicrobial innate immunity gene expression was strongly negatively correlated with keratin and downregulated junction gene expression and strongly positively correlated with mucin and upregulated junction gene expression. The weak correlation between mucin/keratin genes and junction genes suggests that they are regulated separately (Fig. 4G, Spearman’s correlation coefficient). Taken together, we concluded that microbial induction of host interferon is independent from modulation of host antibacterial response and most components of the epithelial barrier, which are more closely coupled in the respiratory epithelium.

In silico identification of potential immunomodulatory microbial pathways

A recurring theme in our examination of epithelial response to microbial colonization was that stimulation of different host pathways was largely independent of phylogeny. It would then follow that there could be numerous mechanisms by which microbes can stimulate pathways such as interferon. We conducted an in silico analysis of the microbial genomes used in this study to identify potential microbial genes/pathways that may stimulate host interferon, comparing gene content differences of ISG non-stimulators and stimulators. While we expected microbial induction of interferon to be poorly conserved, given the lack of phylogenetic similarity between ISG stimulators, we hypothesized that some broad pathways may be conserved.

First, we performed a global analysis, annotating KEGG orthologs using the Functional Mapping and Analysis Pipeline (FMAP)’s database [45] which consists of the UniProt [46] reference 90 cluster filtered for bacteria, fungi, and archaea sequences possessing a KEGG functional ortholog classification. We filtered for orthologs that were absent from either all ISG non-stimulators or stimulators, as a means of rationally subsetting the several thousand hits (Additional file 2: Fig. S11, Additional file 1: Table S13). However, we observed no discernible functional drivers from this broad comparison, likely due to the phylogenetic breadth of our genomes. Based on the ortholog annotation, we looked for known ISG-stimulator metabolites, such short-chain fatty acids [47, 48]; however, we found no associations between metabolite presence and ISG stimulation. Thus, we performed targeted analyses of virulence factors and secondary metabolites, bioactive molecules that are not required for microbial growth and development, because other studies have demonstrated that virulence factors (i.e., outer membrane/cell wall components [49, 50] and antibiotic resistance [51]) and secondary metabolites (i.e., deaminotyrosine [52]) can modulate interferon induction in mononuclear phagocytotic cells.

We quantified the proportion of ISG non-stimulators/stimulators that contained at least one virulence factor from several broad functional categories of virulence factors (via VFDB [53]): adherence, invasion, effector delivery system, motility, exotoxin, exoenzyme, immune modulation, biofilm, nutritional/metabolic factor, stress survival, post-translational modification, antimicrobial activity/competitive advantage, regulation, and others (Fig. 5A, Additional file 1: Table S14). There were few notable differences between ISG non-stimulators and stimulators. However, only ISG stimulators possessed exotoxins and exoenzymes virulence factors, secreted molecules that modulate many host responses including immunity [53, 54]. Furthermore, several of these isolates belonged to genera that also included ISG non-stimulators that did not possess the virulence factor, supporting this factor’s particular association with ISG stimulator rather than phylogeny. While motility was observed in the genomes of ISG stimulators, only one of the ISG stimulators was motile, making it difficult to conclude that motility was associated with ISG stimulation. The differential presence of exotoxins and exoenzymes between ISG non-stimulators and stimulators suggests that at least some microbes induce interferon through secreted factors.

[IMAGE OMITTED: SEE PDF]

Next, we identified broad categories of secondary metabolites, using antiSMASH [55] (fungiSMASH for Candida glabrata) and quantified the proportion of ISG non-stimulator or stimulator genomes that contained each secondary metabolite category (Fig. 5B, Additional file 1: Table S15). Similar to virulence factor annotation, most secondary metabolites were comparably present in both ISG non-stimulators and stimulators. However, there were some secondary metabolite categories that were identified in solely ISG non-stimulators or stimulators. For example, opine-like zincophores were identified in almost all ISG stimulator Staphylococcus spp. but not in the non-stimulator isolates. Zincophores are small molecules that bacteria use to uptake zinc, a necessary cofactor for many enzymes [56]. tRNA-dependent cyclodipeptides synthases, which was only identified in ISG stimulator S. lugundensis, synthesize cyclodipeptides which is a broad category for biodiverse molecules with quorum sensing, antibacterial, and antiviral functions [57]. While the presence of cyclic lactone autoinducer peptides was highly variable between ISG non-stimulators and stimulators, its presence tracked more closely with phylogeny (i.e., if a non-stimulator had a cyclic lactone autoinducer, there were ISG stimulators within the same genera that also possessed it).

A major challenge of this analysis, as we noted, is the phylogenetic breadth which makes it challenging to differentiate genes that differ functionally vs. phylogenetically. We sought to narrow the phylogenetic breadth by examining the Streptococcus genus, the most specific taxonomic category to contain both ISG non-stimulators and stimulators. Repeating the above analyses, we confirmed that few orthologs were shared between several isolates, limiting the utility of the KEGG analysis (Fig. 5C, Additional file 1: Table S16). Exotoxins and exoenzymes were only identified in ISG stimulators and other virulence factors were comparable between ISG stimulators and non-stimulators (Fig. 5D). We also identified several additional secondary metabolite categories that were differentially present between streptococcal ISG non-stimulators and stimulators (Fig. 5E), including lasso peptides, which have antibacterial, antiviral, and protein inhibition effects [58]—which were only found in a Streptococcal ISG stimulator. Given the degree of genetic variability within genera, we also compared the KEGG orthologs of two Rothia dentocariosa strains, one of which was an ISG stimulator and the other a non-stimulator. However, there were no discernible functional differences that would account for the variable stimulation patterns between these two strains.

Taken together, we identified several broad categories of virulence factors and secondary metabolites that are potential contributors to microbial stimulation of ISGs. While the broad phylogenetic diversity that we sampled limited our ability to confidently identify specific potential mechanisms, the variability in bacterial pathways associated with ISG stimulation emphasizes the diversity of potential interaction mechanisms with the epithelium.

Discussion

Here, we present the most extensive analysis to date interrogating the effects of microbial phylogenetic diversity on host-microbe interactions in the respiratory tract, and the most comprehensive analysis of diverse microbial stimulation of epithelial antiviral interferon response. Using a complex 2D epithelial cell culture model, we characterized the effect of colonization with one of 58 phylogenetically diverse, patient-derived bacteria and fungal isolates on epithelial gene expression. We emphasize that this study was designed to examine phylogenetic diverse microbes rather than systematically investigate cohort-specific origins of the microbes (lung cancer patients and some healthy individuals), or other potential variables related to microbial origin, e.g., smoking status (Additional file 1: Table S1). This was because, in our experience, interpersonal genetic variation in the microbiome at the strain level can supersede cohort-specific effects, and our data generally suggested that at the species level, cancer patients (from whom samples were collected prior to systemic therapy) had a similar prevalence of common species as healthy age-matched adults. However, this does not preclude the possibility that the microbes sampled from this cohort differ from those characteristic of matched controls. An additional limitation of our study is potential limited generalizability of the results. Given our goal of characterizing a large microbial panel, we used epithelial cells derived from a single donor to control for donor effects. Future studies are needed to investigate donor effects on microbial response.

Notably, we identified species- and strain-level differences in the modulation of host gene expression, supporting the importance of studying phylogenetically diverse microbes. Strikingly, we also identified a broad, dichotomous ability of phylogenetically diverse microbes to stimulate host interferon that likely functions through different mechanisms. Differential microbial induction of interferon has implications for antiviral immunity/susceptibility microbiome-based therapies for respiratory viral infections, and other immune responses like the antitumor response.

IFN-I and IFN-III, the two types of interferon present in epithelium, are critical for antiviral immunity [7] and, thus, have predominantly been studied in that context. Most studies that have investigated bacterial induction of interferon focused on peripheral blood mononuclear cells, especially dendritic cells and monocytes [50, 51, 59,60,61,62,63]. However, some studies have also demonstrated pathogenic bacteria can stimulate interferon in respiratory epithelium [49, 64, 65]. In this study, we found that half of our microbial treatments (bacterial and fungal) induced epithelial interferon following 48 h of colonization. We note that while we categorized our microbes as interferon non-stimulators versus stimulators, it is likely that this activation is not entirely binary. A few microbes were able to induce interferon by 18 h—in a mechanism unrelated to microbial growth rate—and all stimulator microbes by 48 h. Therefore, it is possible that there are also differences in how quickly a microbe can induce interferon, which has clinical relevance in the ability and rapidity of the microbiome to modulate antiviral immunity. Another striking finding was that a microbe’s potential to stimulate interferon was not phylogenetically driven; we observed species- and strain-level differences, which aligns with other studies [50, 51, 59, 60]. Furthermore, a microbe’s broad categorization as a typical commensal versus pathogen also did not correspond with their ability to induce interferon. For example, both S. aureus isolates, common respiratory pathogens, were interferon stimulators while Klebsiella pneumoniae, another major respiratory pathogen, was a non-stimulator.

In addition, we also observed that a microbe’s potential to stimulate interferon was not dependent upon its potential to stimulate antibacterial immunity, i.e., a microbe that was an interferon stimulator was not necessarily an antibacterial stimulator, and vice versa. This was consistent with another study which found that bacteria that induced IFN-I interferon were poor NF-κB inducers, and vice versa [61]. Induction of antibacterial immunity was moreover driven in part by microbial load, which was unsurprising in that more microbes led to more stimulation of pattern recognition receptors and consequently, an elevated immune response. In contrast, interferon induction did not correspond to microbial load—microbes that grew exponentially and microbes that consistently appeared to die or not grow on the ALIs both induced interferon. Together, this leads us to believe that microbes can induce host interferon through different mechanisms than used to induce antibacterial immunity.

Our results are notable because there has traditionally been a divide between bacterial pattern recognition receptors, which recognize extracellular microbes and initiate antibacterial immune responses, and viral pattern recognition receptors, which target intracellular microbes and induce antiviral immune responses [66]. Interferon, as an antiviral immune response, has canonically been associated with intracellular viral pattern recognition receptors. Therefore, microbial induction of epithelial interferon may be a previously uncharacterized function of traditional antibacterial immune response, a previously uncharacterized means of stimulating viral pattern recognition receptors, or an uncharacterized mechanism. Several different mechanisms of bacterial induction of interferon have been demonstrated. For example, outer membrane polysaccharide A has been shown to induce IFN-I through TLR4 [50]. Bacterial DNA can induce IFN-I through TLR9 [62], and bacterial nucleic acid can enhance IFN-I through TLR3, TLR9 [59], cGAS/STING [64], and RIGI-like receptors [61]. Secreted metabolites such as desaminotyrosine (DAT) [52] and secondary bile acids [63] also have been shown to promote IFN-I signaling. However, all of these mechanisms were demonstrated through microbial interactions with phagocytotic immune cells and may not translate to epithelial cells. Furthermore, given that IFN stimulators were not necessarily phylogenetically similar, it is plausible that there are numerous mechanisms through which microbes induce epithelial interferon.

While we underscore that the genetic heterogeneity of our microbial set is so broad that identifying ISG-specific genetic associations is difficult, we identified genes present in only ISG-stimulators, including exotoxins, exoenzymes, opine-like zincophores, and tRNA-dependent cyclodipeptide synthase. This heterogeneity complicated comparative efforts even down to the genus level; thus, we purport that an even higher resolution analysis, such as examining strain variation within a species, will be needed to further distinguish genetic differences between ISG stimulators and non-stimulators. While our strain-level investigation did not identify functional differences that would account for the differential stimulation of interferon, it is likely that this is due to insufficient power at the strain level, given our overarching focus on species-level diversity.

In conclusion, we demonstrated that microbial induction of interferon is highly heterogeneous, not phylogenetically driven, and occurs independently of microbial induction of more general inflammatory processes. This dataset builds a foundation for understanding how respiratory microbiome biodiversity affects microbe-epithelial interactions, particularly in the context of antiviral innate immunity. Subsequent studies will further interrogate the biological consequences of these transcriptional changes, particularly in the context of viral infection, and the degree of conservation of the bacterial mechanisms of stimulating changes in host antiviral innate immunity.

Conclusions

We showed that colonization of the respiratory epithelium with patient-derived microbes differentially stimulates host gene expression. When looking at innate immunity and epithelial barrier genes, we saw variation in host gene expression following colonization with phylogenetically diverse microbes and following colonization with different strains of the same species. Strikingly, some of the microbes stimulated host expression of antiviral interferon. However, this host response was not driven by microbial phylogeny nor a microbe’s propensity to stimulate canonical antibacterial innate immunity. Finally, we identified some potential mechanisms behind microbial stimulation of interferon. This study significantly expands our understanding of how the respiratory microbiome’s biodiversity drives and modulates epithelium-microbiome interactions, especially host innate immune responses.

Methods

Lung cancer patient recruitment

We recruited 47 patients (46–83 years old) with established evidence of adenocarcinoma lung cancer. Written informed consent was obtained from all participants and the study protocol was approved by the Human Subjects Institutional Review Board at Hartford Hospital, Hartford, USA, and performed in the context of U19AI142733 grant at the Jackson Laboratory.

Sample collection

Nares and oral samples were obtained using sterile swabs (Puritan™ PurFlock Ultra, #22–029–506, Guilford, ME, USA). For the nares, two swabs pre-moistened in nuclease-free water (Qiagen, Hilden, Germany) were inserted 2 cm into one nostril and rotated against the anterior nasal mucosal epithelium for up to 10 s. The tongue dorsum was gently rubbed with two dry swabs for up to 30 s.

One swab each was stored in 350 μl Tissue and Cell lysis solution (Lucigen, #MTC096H, Middleton, WI) lysis buffer with 100-μl glass beads (0.1 diameter, BioSpec Products, Cat No. 11079101, Bartlesville, OK, USA). The other swab was placed and stored in R2A (Innovation Diagnostics, LAB203-A) for future microorganism recovery. All samples were stored at − 80 °C until shipping. Samples were shipped on dry ice to The Jackson Laboratory for Genomic Medicine, Farmington, Connecticut, USA, and stored at − 80 °C until DNA extraction or cultivation.

mWGS DNA extraction

Metagenomic whole genome shotgun sequencing (mWGS) allows for comprehensive sampling of all genes in all organisms present in a given complex sample and subsequent identification of the bacteria, viruses, and fungi present. Genomic DNA was extracted from frozen samples using the GenElute Bacterial DNA Isolation kit (Sigma Aldrich, NA2110-1KT, St. Louis, MO, USA) with the following minor modifications as previously described [67]: Briefly, 5 μl of Lysozyme (10 mg/ml, Sigma-Aldrich, L6876, St. Louis, MO, USA), 1 μl of Lysostaphin (5000 U/ml, Sigma-Aldrich, L9043, St. Louis, MO, USA), and 1 μl mutanolysin (5000 U/ml, Sigma-Aldrich, M9901, St. Louis, MO, USA) were added to each sample for a digest at 37 °C, 30 min. Samples were mechanically disrupted 2 × for 3 min at 30 Hz (TissueLyser II, Qiagen, Hilden, Germany). 5 μl of proteinase K (20 mg/ml, Sigma-Aldrich), and 300 μl of Solution C (55 °C, 30 min) was added. The samples were precipitated by adding 300 μl of 100% EtOH (Fisher Scientific, Fair Lawn, NJ, USA), and the lysates were loaded on the GenElute columns. Subsequent steps were carried out according to the manufacturer’s instructions.

mWGS DNA sequencing

Libraries were prepared using the Nextera XT DNA Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions, however using one-quarter of each reaction volume. Whole genome sequencing (WGS) was carried out using a 2 × 150 bp (paired-end) sequencing protocol for the Illumina NovaSeq 6000 Sequencing System (Illumina, San Diego, CA, USA) according to the manufacturer’s manual. Sequencing was conducted at the Genome Technologies core facility at the Jackson Laboratory for Genomic Medicine, Farmington, USA.

Positive and negative controls

Per patient, one air sample (negative control) was collected by waving a moist swab through the air, to collect potential environmental contaminants, and extracted as described above. Samples that yielded a Qubit measurement were processed for library preparation and sequenced as described. Per extraction round, one sample of a defined, in-house mock community (25 diverse Gram-positive and Gram-negative bacteria and fungi) and a negative reagent control (nuclease-free water, Qiagen, Hilden, Germany) were included. For library generation, a negative control (nuclease-free water, Qiagen, Hilden, Germany) was included as well. One mock community sample was added to each sequencing run and a library/extraction negative control was sequenced if a library product was measurable on the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA) or identified on the 4200 TapeStation System (Agilent Technologies, Santa Clara, CA) using the High Sensitivity D1000 ScreenTape Assay (Agilent Technologies, Santa Clara, CA).

mWGS data processing for taxonomic classification

The mWGS data analysis includes both taxonomic and functional profiling from complex microbial communities. After demultiplexing, the fastq files were processed with MetaPhlAn 4 [68] using the flags –add_viruses and –unclassified_estimation. Samples with no reads mapping to any taxa were excluded.

Relative proportions were used for all analyses. All taxonomic features at the species level with a mean relative abundance < 0.0001% (denoise function [69]) across the dataset were removed from the dataset to reduce potential false positives and allow for multiple hypothesis correction (except for the comparison with the culturomics data).

Isolate culturing and identification

Oral and nares samples that would be used for culturing were thoroughly vortexed and then diluted 1:100 and 1:1000 in R2A, to increase the chance of recovering single colonies. Fifty microliters from each dilution was then spread on half of an agar plate (R2A, LB, TSA blood agar, Chocolate agar, aerobic only: SaSelect) for each cultivation condition using a sterile spread tool (Thomas Scientific, 229,616) as previously described [28]. Briefly, for agar cultivation, the following media were sourced from Fisher Scientific: Luria Broth (LB) agar (BP1425500), R2A agar (R454372), Tryptic Soy Broth (TSB) (DF0370-17–3), Chocolate agar (B21169X), Tryptic Soy Agar (TSA) with 5% sheep blood agar (B21261X). SaSelect agar plates (63748) were obtained from BioRad. The anaerobic atmosphere consisted of 5% hydrogen, 5% carbon dioxide, and 90% nitrogen (Airgas, Z03NI9022000008). Aerobic cultures were conducted in ambient atmosphere.

Microbial isolates were identified as previously described [28] using matrix assisted laser desorption ionization-time of flight (MALDI-TOF, Bruker Daltonics, Germany) mass spectrometry. Rapid DNA extraction from microbial isolates was adapted from Köser et al. [70] as previously described [28] with the following modifications. A 2-mL overnight culture was centrifuged at 20,000 × g for 1 min until a bacterial pellet was formed. The pellet was then resuspended in 150 µL of 10 mM Tris and transferred to a 2-mL bead-beating tube containing 150 µL of 10 mM Tris and 100 µL 0.5-mm diameter glass beads (BioSpec Products, NC0417355). The microbial DNA was then sequenced as previously described [28] on the NovaSeq 6000 Sequencing System (Illumina, San Diego, CA, USA), 2 × 150 bp paired end to approximately 100X coverage per genome. Following sequencing, genomic reads were preprocessed and dereplicated using PRINSEQ lite [71] (-derep 1) and trimmed using trimmomatic [72]. Reads were assembled into contigs using SPAdes [73]. QUAST [74] was used to assess contig quality. All assembly steps used default parameters unless otherwise noted. A reference genome was used for Candida glabrata (GCF_000002545.3). Taxonomic classification was assigned to the contigs and the contigs were placed in a phylogenetic tree using GTDB [75, 76] using default parameters. The phylogenetic tree was visualized using the R package ggtree [77].

Statistics and sample comparisons for metagenomics and culturomics data

Overlap of metagenomics and culturomics data was calculated extracting the uniquely identified species of genera in the culturomics data and metagenomics data. Data was analyzed and visualized using the following R packages: reshape2 [78], ggplot2 [79], tidyverse [80], ggpubr [81], dplyr [82], plyr [83], and ggvenn [84].

Comparative genomic analysis

Gene coding sequences of bacterial genomes were predicted using Prokka [85] and Funannotate [86] was used to sort and predict gene coding sequences for the Candida glabrata genome, all with default parameters.

KEGG orthologs were annotated by blasting the microbial genomes against Functional Mapping and Analysis Pipeline (FMAP)’s database [45] (downloaded in 2020), which consists of the UniProt [46] reference cluster filtered for bacteria, fungi, and archaea sequences possessing a KEGG functional classification, using UBLAST [87]. A minimum e-value of 10−9 was required. KEGG ortholog hits were subsequently filtered for a minimum percent identity of 80%, a minimum of 80% coverage of the query and target sequence, differential ortholog presence between the microbial genomes, and absent from either all ISG stimulator or non-stimulator genomes.

Microbial genomes were annotated for virulence factors by blasting them against VFDB’s protein database A (verified virulence factors, downloaded in 2023) [53] using UBLAST [87] with a e-value of 10−9. Hits were subsequently filtered for a minimum percent identity of 80% and a minimum of 80% coverage of the query and target sequence.

Microbial files were annotated for secondary metabolites using antiSMASH [55] (fungiSMASH for the fungal isolate) using strictness “relaxed.” antiSMASH annotates secondary metabolite gene clusters by identifying experimentally characterized proteins and filters for gene clusters that include the minimal core components of each gene cluster. –genefinding-tool was set to use Prodigal [88] if any gene was lacking an annotation. Defaults were used for other parameters.

Data was analyzed and visualized using the following R packages: ggplot2 [79], ggpubr [81], dplyr [82], rstatix [89], and RColorBrewer [90], tidyverse [80], and pheatmap [91].

Air–liquid interface cell culture cultivation

9 mm primary normal human bronchial epithelial air–liquid interface (ALI) tissue cultures were obtained from MatTek Corporation (EpiAirway AIR-100, Gothenburg, Sweden). The ALI cultures had been switched to antibiotic-free media 3 days prior to receipt to facilitate microbial growth. All batches were grown using cells from MatTek Corporation’s EpiAirway AIR-100 standard donor, a healthy adult male. ALI cultures were cultured according to the manufacturer’s directions. Briefly, upon arrival, the ALI cultures were placed in 6-well plates with 1 mL of warmed antibiotic-free EpiAirway AIR-100 Maintenance Media (MatTek Corporation, Gothenburg, Sweden) or the equivalent EpiAirway AIR-100 Assay Media (MatTek Corporation, Gothenburg, Sweden) basally per well. The basal media was replaced daily, and the ALI cultures were kept at 37 °C with 5% CO2.

ALI treatments

For each microbial isolate, a single colony was grown overnight in sterile 1X Tryptic Soy Broth (TSB, Becton, Dickinson, and Company, #211,825, Sparks, MD) with 0.1 mg Vitamin K (Sigma Aldrich, #95,271, St. Louis, MO) and 5 mg heme / 1L TSB. For isolates that could not be grown from a single colony, a single colony was patched, and that patch was used to start a liquid culture. ODs were obtained by measuring absorbance at 600 nM in a 96-well plate using Cytation Station 5 (BioTek Agilent Technologies, Santa Clara, CA). 108 colony-forming units (CFUs) were taken from each liquid culture and washed with ultrapure water (Fisher Scientific, #AAJ71786AP, Hampton, NH) being resuspended in 100 µL of ultrapure water (Fisher Scientific, #AAJ71786AP, Hampton, NH), with a final concentration of 107 CFUs per 10 µL.

Just prior to treatment, the ALI cultures were quickly washed twice with transepithelial electrical resistance buffer to remove excess mucus, where the wash buffer was immediately removed without an incubation period. Each wash consisted of the addition and immediate removal of 400 µL of the buffer. ALI cultures were then dosed with 10 µL of microbial isolate/vehicle (ultrapure water, Fisher Scientific, #AAJ71786AP, Hampton, NH). Microbial isolates were dosed at 107 CFUs. Dosed ALI cultures were incubated for 18 or 48 h prior to harvest. Extra inoculum was serially diluted in sterile PBS (MatTek Corporation, Gothenburg, Sweden) and grown on TSA with 5% sheep’s blood (Fisher Scientific, #221,261, Hampton, NH) to determine the number of microbes added to the ALI cultures. To minimize the risk of confounding batch effect, replicates for each microbial treatment were spread across batches to distribute potential variance. Therefore, sample size indicates the number of independent experiments.

At harvest, 200 µL of transepithelial electrical resistance buffer (MatTek Corporation, Gothenburg, Sweden) was added to the apical surface of each ALI culture, pipette mixed, and removed. One hundred forty microliters of Buffer RLT (Qiagen, Hilden, Germany) + 1% beta-mercaptoethanol was added to each ALI culture, dissolving the tissue. The Buffer RLT-tissue solution was frozen at − 80 °C until RNA extraction. S. epidermidis Tü3298-GFP-colonized ALI were visualized under blue light for a qualitative examination of colonization. The wash was serially diluted in PBS (MatTek Corporation, Gothenburg, Sweden) and plated on TSA with 5% sheep’s blood (Fisher Scientific, #221,261, Hampton, NH) or 1X TSB (Becton, Dickinson, and Company, #211,825, Sparks, MD) with 1X Bacto Dehydrated Agar (Fisher Scientific, #214,010, Hampton, NH) for CFU counts. Basal media was collected and frozen at − 80 °C for cytokine bead array assays.

RNA extraction and RNA-seq

All RNA extraction and sequencing library preparation steps were performed in a sterile tissue culture hood. RNA was extracted using RNeasy 96 QIAcube HT kit (Qiagen, Hilden, Germany) according to the manufacturer’s directions. Samples were eluted in nuclease-free water (Qiagen, Hilden, Germany) and frozen at − 80 °C until sequencing preparation. RNA quality was evaluated using the 4200 TapeStation System (Agilent Technologies, Santa Clara, CA) with the High Sensitivity RNA ScreenTape Analysis (Agilent Technologies, Santa Clara, CA) or RNA ScreenTape Analysis (Agilent Technologies, Santa Clara, CA). RINs ranged from 1.7 to 9.9 with an average of 8.4 and a median of 8.9. RNA quantity was measured on the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA). The sequencing libraries were prepared with NEBNext rRNA Depletion Kit v2 (New England Biolabs, Ipswich, MA) and NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA) following the manufacturer’s directions. Library quality was evaluated using the 4200 TapeStation System (Agilent Technologies, Santa Clara, CA) with the High Sensitivity D1000 ScreenTape Assay (Agilent Technologies, Santa Clara, CA). Library quantity was measured on the Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA). Samples were sequenced using Illumina NovaSeq targeting 30 million reads. Reads per sample ranged from 76 thousand to 318 million with an average and median of 29 million reads.

Transcriptional profiling, gene set enrichment analysis, and CFU confounder analyses

Raw RNA-seq reads were processed with trimmomatic 0.39 [72] to remove low-quality reads. Reads were then mapped to T2T-CHM13v2.0 reference genome using STAR 2.5.3a [92]. The raw read counts were calculated using featureCounts (from Subread1.6.4) [93]. The raw read counts were normalized by RUVg [94] within each comparison group of microbe-treated samples versus vehicle controls using 5000 of the most stably expressed genes. Differentially expressed genes were identified and expression level changes of each gene were computed using DESeq2 [95] with shrinkage by rlogTransformation and DEG threshold of |log2FC|> 1 and adjusted p-value < 0.05. Gene set enrichment analysis was performed on pre-ranked gene lists (sign of log2FC * 2|log2FC| * − log10(p-value)) calculated from output of DESeq2 on ALI transcriptomes by GSEA4.3.2 [96, 97] with gene ontology (GO) database or clusterProfiler [98, 99] with the gseGO method, a p-value cut-off of 0.05, and FDR p-value adjustment. CFU confounder analyses were generated by MaAsLin2 [34] with TSS (total sum-scaling) normalization. Fixed effects included CFUs (at the time of harvest), ALI batch, and timepoint. CFUs were log10 transformed prior to MaAsLin2. An effect was considered a confounder if the FDR adjusted P-value < 0.05.

Gene list generation and determination of responsive genes in RNA-seq data

ISG gene lists were adapted from the literature: total ISG list [100], 6-gene ISG list [36], and 30-gene ISG list [37]. The antibacterial gene list was derived from MSigDB [38] with gene ontologies [39, 40] (GO:0005125, GO: 0061844, GO: 0019730, GO: 0050832, GO: 0009620). The keratin gene list was adapted from the literature [43]. The junction and mucin gene lists were derived from HGNC (HUGO Gene Nomenclature Committee) [44]. The junction gene list was further adapted from the literature [42]. Within each gene list, genes were further processed by removing genes with raw read count lower than 5. The remaining genes were further subselected with hierarchical clustering using the seaborn Python package [101]. Clusters of genes that were enriched with increased expression, as determined by high log2FC, were defined as responsive genes. The degree of increased expression depended on the gene list.

Statistics and sample comparisons for RNA-seq data

Data was analyzed and visualized using the following Python packages: seaborn [101], scikit-learn [102], statannot [103], pandas [104], numpy [105], and matplotlib [106]; and R packages: ComplexHeatmap [107], reshape2 [78], and ggplot2 [79].

Cytokine measurement by Cytometric Bead Array (CBA)

Each cultured microbe was grown overnight in TSB with 0.1 mg vitamin K and 5 mg heme / 1 L TSB in a 96-well deep well plates (#503,501, Southern Labware, Cumming, GA, USA). ODs were obtained by measuring absorbance at 600 nM in a 96-well plate using Cytation Station 5 (BioTek Agilent Technologies, Santa Clara, CA). The liquid cultures were centrifuged at 3500 rpm for 5 min (Sorvall Legend X1R M20 rotor, Thermo Scientific, Waltham, MA, USA) to separate the bacterial pellet and bacterial supernatant. The bacterial supernatant was collected and sterile filtered (0.2 μm) using a filter plate (#MSGVS2210, Sigma Aldrich, St. Louis, MO, USA) at 3500 rpm for 5 min (Sorvall Legend X1R M20 rotor, Thermo Scientific, Waltham, MA, USA). The sterile bacterial supernatants were stored at − 80 °C until subsequent use. Sterile bacterial supernatants and ALI basal media were shipped at − 80 °C to Washington University School of Medicine, St. Louis, MO.

For conditioned bacterial supernatant suppression/induction of cytokines, lung epithelial cell line A549 was stimulated overnight with sterile (0.2 μm filtered) bacterial supernatants. Briefly, 1 × 105 cells per well were with 20% volume/volume of each supernatant overnight. The conditioned A549 cell media was collected for CBA analysis.

CBA analysis was conducted using six bead populations with distinct fluorescence intensities that were coated with capture antibodies specific for IL-8, IL-6, IL-10, IL-1, TNF, and IL-12p70 proteins (Becton, Dickinson, and Company BioSciences, #551,811, Sparks, MD). They were incubated together with the samples (sterile bacterial supernatant or microbially colonized ALIs’ basal media), negative controls (sterile bacterial growth media or basal media from vehicle treated ALI), or recombinant standards, and PE-conjugated detection antibodies, to form sandwich complexes and were detected by flow cytometry. Results were generated in graphical and tabular format using the BD CBA Analysis Software (Becton, Dickinson, and Company BioSciences, Sparks, MD). The standard curve for each protein covers a defined set of concentrations from 20 to 5000 pg/ml.

Data was analyzed and visualized using the following R packages: ggplot2 [79], ggpubr [81], dplyr [82], rstatix [89], tidyverse [80], forcats [108], and RColorBrewer [90].

Immunofluorescence staining

ALI cultures were washed with 1X PBS (MatTek Corporation, Gothenburg, Sweden), then embedded in OCT (Sakura Finetek USA, Torrance, CA, USA) and snap-frozen at − 80 °C. Frozen sections were cut at 8 µm, air dried on Superfrost plus slides, fixed with 4% paraformaldehyde (#28,906, Thermo Fisher Scientific, Waltham, MA, USA) for 15 min, then permeabilized with 1X PBS/0.1% Triton-X (#HFH10, Thermo Fisher Scientific, Waltham, MA, USA) for 15 min. Tissue sections were treated with Fc Receptor Block (#NB309, Innovex Bioscience, Richmond, CA, USA), followed by Background Buster (#NB306, Innovex Bioscience, Richmond, CA, USA). The sections were then stained with anti-MX1 primary antibody (polyclonal Rabbit N2C2, #GTX110256, GeneTex, Irvine, CA, USA) for 1 h followed by secondary antibody (anti-rabbit IgG AF555, #406,412, BioLegend, San Diego, CA, USA) for 30 min at room temperature in 1X PBS/5% BSA/0.05% saponin. Then, sections were washed three times with 1X PBS for 15 min. Finally, sections were counterstained with 1 µg/ml of 4',6-diamidino-2-phenylindole (DAPI, #D1306, Thermo Fisher Scientific, Waltham, MA, USA) then mounted with Fluoromount-G (#00–4958-02, Thermo Fisher Scientific, Waltham, MA, USA), acquired using a Leica SP8 confocal microscope (Leica Microsystems, Wetzlar, Germany) for high-resolution images and Thunder widefield microscope (Leica Microsystems, Wetzlar, Germany) for quantification both with Leica LAS X software (Leica Microsystems, Wetzlar, Germany) and analyzed using Imaris software (Bitplane, Oxford Instruments, Abingdon, United Kingdom). Following Imaris image analysis, the intensity mean of each marker was quantified using histocytometry (FlowJo).

Data was analyzed and visualized using the following R packages: ggplot2 [79], ggpubr [81], dplyr [82], rstatix [89], and RColorBrewer [90]. All results are expressed as a mean of the replicates, unless specified. All comparisons were made between infection conditions with time point-matched, uninfected controls.

Transepithelial electrical resistance (TEER) measurements

TEER measurements were taken using EVOM Manual for TEER Measurement (#EVM-MT-03–01, WPI, Sarasota, FL, USA) and STX4 EVOM Electrode (#EVM-EL-03–03-01, WPI, Sarasota, FL, USA) following 18 and 48 h of colonization. The electrodes were equilibrated following the manufacturer’s instructions. ALIs were transferred to a 12-well plate. One milliliter of media was added to the basal side and 300 µL of TEER buffer was added to the apical surface. Three readings were taken of each insert, evenly spread across the insert. In between each reading, the electrodes were washed in 70% ethanol and then washed in TEER buffer.

Data was analyzed and visualized using the following R packages: ggplot2 [79], ggpubr [81], dplyr [82], rstatix [89], and RColorBrewer [90].

LDH quantification

LDH was quantified from ALI apical media, basal media, and cell lysate using CyQuant LDH Cytotoxicity Assay (Thermo Fisher Scientific, #C20301, Waltham, MA) following the manufacturer’s directions but with the following modifications. ALI were colonized for 48 h as described above. Following colonization, 200 µL of TEER buffer was added to the apical surface of each ALI. ALI were subsequently incubated at 37 °C for 10 min. Following incubation, the TEER buffer was pipetted up and down 5 times prior to removal to maximize the amount of bacteria removed and minimize confounding signal. The apical wash was saved to quantify apical LDH signal. An aliquot of the basal media was also saved to determine basal LDH secretion. Three hundred microliters of 5X lysis buffer was then added to the apical surface of each ALI and ALI were incubated for at 37 °C for 1.5 h. In the middle of the incubation, a pipette tip was used to mechanically disrupt the ALI to ensure thorough cell lysis. LDH was quantified in the basal media, apical wash, and ALI cell lysate. Samples were run in a 384-well plate format and the signal was read using Cytation Station 5 (BioTek Agilent Technologies, Santa Clara, CA). Apical wash and basal media were diluted 1:2 and cell lysate 1:20 in TEER buffer.

Additional ALI, following the apical wash, were each transferred to a 1.5-mL Eppendorf tube with ~ 200 µL of TEER buffer and 100 µL of 1 mm Zironia/Silica beads (BioSpec Products, #11079110z, Barlesville, OK) and were mechanically homogenized for 3 min at 30 Hz using TissueLyser II (Qiagen, Hilden, Germany). The cell lysates were then serially diluted and plated on TSA to quantify live, adherent bacteria for each microbial treatment.

To normalize the cell lysate LDH signal for any confounding bacterial LDH that was released during lysis, fresh TSB cultures were grown for each microbial treatment. Bacterial pellets were washed with sterile water to remove excess TSB, and then each microbial treatment was serially diluted to create to a standard curve ranging 104–109 CFUs. The bacterial standard curves were incubated for 1.5 h at 37 °C in 100 µL of 10X lysis buffer and then LDH from the bacterial lysate was measured. A higher concentration of lysis buffer than was used for the ALI was used for the bacteria to determine the maximal contribution of bacterial LDH to the cell lysate LDH.

Bacterial LDH signal was subtracted from the cell lysate LDH, based on the CFUs determined from the cell homogenate. Replicates were averaged together. Data was analyzed and visualized using the following R packages: ggplot2 [79], dplyr [82], rstatix [89], and RColorBrewer [90].

Data availability

All data can be accessed through ImmPort accession number: SDY2281. Metagenomic sequence files can also be accessed in the NCBI Sequence Read Archive BioProject PRJNA981523 and RNA-sequencing files can be accessed in NCBI database of Genotypes and Phenotypes (dbGaP) BioProject phs003627. Analyzed gene tables and code are available on our GitHub (ohlab).

Abbreviations

ALI:

Air-liquid interface cell culture

S. :

Staphylococcus

Str. :

Streptococcus

spp.:

Species

CBA:

Cytokine bead array

CFU:

Colony-forming unit

PCA:

Principal component analysis

DEGs:

Differentially expressed genes

IFN-I:

Type 1 interferon

IFN-III:

Type 3 interferon

ISG:

Interferon stimulated gene

Log2FC:

Log2 fold change

L. :

Lactobacillus

TEER:

Transepithelial electrical resistance

FMAP:

Functional Mapping and Analysis Pipeline

VFDB:

Virulence Factor Database

TLR:

Toll-like receptor

DAT:

Desaminotyrosine

TSA:

Tryptic soy agar

TSB:

Tryptic soy broth

GSEA:

Gene set enrichment analysis

NES:

Normalized enrichment score

A. :

Actinomyces

B. :

Bacillus

C. :

Corynebacterium

D. :

Dermacoccus

E. :

Enterococcus

K. :

Klebsiella

R. :

Rothia

Tü:

Tü3298GFP

References

Whitsett JA, Alenghat T. Respiratory epithelial cells orchestrate pulmonary innate immunity. Nat Immunol. 2015;16:27–35.

Leiva-Juárez MM, Kolls JK, Evans SE. Lung epithelial cells: therapeutically inducible effectors of antimicrobial defense. Mucosal Immunol. 2018;11:21–34.

Dajani R, Zhang Y, Taft PJ, Travis SM, Starner TD, Olsen A, et al. Lysozyme Secretion by Submucosal Glands Protects the Airway from Bacterial Infection. Am J Respir Cell Mol Biol. 2005;32:548–52.

Bals R, Wang X, Wu Z, Freeman T, Bafna V, Zasloff M, et al. Human beta-defensin 2 is a salt-sensitive peptide antibiotic expressed in human lung. J Clin Invest. 1998;102:874–80.

Goldman MJ, Anderson GM, Stolzenberg ED, Kari UP, Zasloff M, Wilson JM. Human beta-defensin-1 is a salt-sensitive antibiotic in lung that is inactivated in cystic fibrosis. Cell. 1997;88:553–60.

Bals R, Wang X, Zasloff M, Wilson JM. The peptide antibiotic LL-37/hCAP-18 is expressed in epithelia of the human lung where it has broad antimicrobial activity at the airway surface. Proc Natl Acad Sci. 1998;95:9541–6.

Lazear HM, Schoggins JW, Diamond MS. Shared and distinct functions of type I and type III interferons. Immunity. 2019;50:907–23.

Denney L, Byrne AJ, Shea TJ, Buckley JS, Pease JE, Herledan GMF, et al. Pulmonary epithelial cell-derived cytokine TGF-β1 Is a CRITICAL Cofactor for enhanced innate lymphoid cell function. Immunity. 2015;43:945–58.

Subauste MC, Proud D. Effects of tumor necrosis factor-α, epidermal growth factor and transforming growth factor-α on interleukin-8 production by, and human rhinovirus replication in, bronchial epithelial cells. Int Immunopharmacol. 2001;1:1229–34.

Allakhverdi Z, Comeau MR, Jessup HK, Yoon B-RP, Brewer A, Chartier S, et al. Thymic stromal lymphopoietin is released by human epithelial cells in response to microbes, trauma, or inflammation and potently activates mast cells. J Exp Med. 2007;204:253–8.

The Human Microbiome Project Consortium. Structure, function and diversity of the healthy human microbiome. Nature. 2012;486:207–14.

Escapa IF, Chen T, Huang Y, Gajare P, Dewhirst FE, Lemon KP. New insights into human nostril microbiome from the expanded human oral microbiome database (eHOMD): a resource for the microbiome of the human aerodigestive tract. mSystems. 2018;3. https://doi.org/10.1128/msystems.00187-18.

Natalini JG, Singh S, Segal LN. The dynamic lung microbiome in health and disease. Nat Rev Microbiol. 2023;21:222–35.

Bassis CM, Erb-Downward JR, Dickson RP, Freeman CM, Schmidt TM, Young VB, et al. Analysis of the upper respiratory tract microbiotas as the source of the lung and gastric microbiotas in healthy individuals. mBio. 2015;6:e00037.

Han X, Lee A, Huang S, Gao J, Spence JR, Owyang C. Lactobacillus rhamnosus GG prevents epithelial barrier dysfunction induced by interferon-gamma and fecal supernatants from irritable bowel syndrome patients in human intestinal enteroids and colonoids. Gut Microbes. 2019;10:59–76.

Kim HJ, Jo A, Jeon YJ, An S, Lee K-M, Yoon SS, et al. Nasal commensal Staphylococcus epidermidis enhances interferon-λ-dependent immunity against influenza virus. Microbiome. 2019;7:80.

Miao P, Jiang Y, Jian Y, Shi J, Liu Y, Piewngam P, et al. Exacerbation of allergic rhinitis by the commensal bacterium Streptococcus salivarius. Nat Microbiol. 2023;8:218–30.

Gulraiz F, Rellinghausen C, Bruggeman CA, Stassen FR. Haemophilus influenzae increases the susceptibility and inflammatory response of airway epithelial cells to viral infections. FASEB J. 2015;29:849–58.

Ratner AJ, Lysenko ES, Paul MN, Weiser JN. Synergistic proinflammatory responses induced by polymicrobial colonization of epithelial surfaces. Proc Natl Acad Sci. 2005;102:3429–34.

Ortiz Moyano R, Raya Tonetti F, Tomokiyo M, Kanmani P, Vizoso-Pinto MG, Kim H, et al. The Ability of Respiratory Commensal Bacteria to Beneficially Modulate the Lung Innate Immune Response Is a Strain Dependent Characteristic. Microorganisms. 2020;8:E727.

Wang J, Li F, Sun R, Gao X, Wei H, Li L-J, et al. Bacterial colonization dampens influenza-mediated acute lung injury via induction of M2 alveolar macrophages. Nat Commun. 2013;4:2106.

Beresford-Jones BS, Forster SC, Stares MD, Notley G, Viciani E, Browne HP, et al. The Mouse Gastrointestinal Bacteria Catalogue enables translation between the mouse and human gut microbiotas via functional mapping. Cell Host Microbe. 2022;30:124–138.e8.

Jiang D, Schaefer N, Chu HW. Air–liquid interface culture of human and mouse airway epithelial cells. In: Alper S, Janssen WJ, editors. Lung innate immunity and inflammation. New York: Springer New York; 2018. p. 91–109.

Leung C, Wadsworth SJ, Yang SJ, Dorscheid DR. Structural and functional variations in human bronchial epithelial cells cultured in air-liquid interface using different growth media. Am J Physiol-Lung Cell Mol Physiol. 2020;318:L1063–73.

Pezzulo AA, Starner TD, Scheetz TE, Traver GL, Tilley AE, Harvey B-G, et al. The air-liquid interface and use of primary cell cultures are important to recapitulate the transcriptional profile of in vivo airway epithelia. Am J Physiol Lung Cell Mol Physiol. 2011;300:L25–31.

Huffnagle GB, Dickson RP, Lukacs NW. The respiratory tract microbiome and lung inflammation: a two-way street. Mucosal Immunol. 2017;10:299–306.

Wilbert SA, Mark Welch JL, Borisy GG. Spatial Ecology of the Human Tongue Dorsum Microbiome. Cell Rep. 2020;30:4003-4015.e3.

Fleming E, Pabst V, Scholar Z, Xiong R, Voigt AY, Zhou W, et al. Cultivation of common bacterial species and strains from human skin, oral, and gut microbiota. BMC Microbiol. 2021;21:278.

Baggiolini M, Clark-Lewis I. Interleukin-8, a chemotactic and inflammatory cytokine. FEBS Lett. 1992;307:97–101.

Tu3298 Genome. NCBI GenBank. 2016. http://identifiers.org/bioproject:PRJEB11651.

Moran JC, Horsburgh MJ. Whole-Genome Sequence of Staphylococcus epidermidis Tü3298. Genome Announc. 2016;4:e00112-e116.

Staphylococcus epidermidis Genome sequencing and assembly. NCBI GenBank. 2020. http://identifiers.org/bioproject:PRJNA559376.

Zhou W, Spoto M, Hardy R, Guan C, Fleming E, Larson PJ, et al. Host-Specific Evolutionary and Transmission Dynamics Shape the Functional Diversification of Staphylococcus epidermidis in Human Skin. Cell. 2020;180:454–70.e18.

Mallick H, Rahnavard A, McIver LJ, Ma S, Zhang Y, Nguyen LH, et al. Multivariable association discovery in population-scale meta-omics studies. PLOS Comput Biol. 2021;17:e1009442.

Schoggins JW, Rice CM. Interferon-stimulated genes and their antiviral effector functions. Curr Opin Virol. 2011;1:519–25.

Rice GI, Forte GMA, Szynkiewicz M, Chase DS, Aeby A, Abdel-Hamid MS, et al. Assessment of interferon-related biomarkers in Aicardi-Goutières syndrome associated with mutations in TREX1, RNASEH2A, RNASEH2B, RNASEH2C, SAMHD1, and ADAR: a case-control study. Lancet Neurol. 2013;12:1159–69.

Reynolds JA, Briggs TA, Rice GI, Darmalinggam S, Bondet V, Bruce E, et al. Type I interferon in patients with systemic autoimmune rheumatic disease is associated with haematological abnormalities and specific autoantibody profiles. Arthritis Res Ther. 2019;21:147.

Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417–25.

Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene Ontology: tool for the unification of biology. Nat Genet. 2000;25:25–9.

The Gene Ontology Consortium, Aleksander SA, Balhoff J, Carbon S, Cherry JM, Drabkin HJ, et al. The Gene Ontology knowledgebase in 2023. Genetics. 2023;224:iyad031.

Gao N, Rezaee F. Airway Epithelial Cell Junctions as Targets for Pathogens and Antimicrobial Therapy. Pharmaceutics. 2022;14:2619.

Yuksel H, Turkeli A. Airway epithelial barrier dysfunction in the pathogenesis and prognosis of respiratory tract diseases in childhood and adulthood. Tissue Barriers. 2017;5:e1367458.

Jacob JT, Coulombe PA, Kwan R, Omary MB. Types I and II Keratin Intermediate Filaments. Cold Spring Harb Perspect Biol. 2018;10:a018275.

Seal RL, Braschi B, Gray K, Jones TEM, Tweedie S, Haim-Vilmovsky L, et al. Genenames.org: the HGNC resources in 2023. Nucleic Acids Res. 2023;51:D1003-9.

Kim J, Kim MS, Koh AY, Xie Y, Zhan X. FMAP: Functional Mapping and Analysis Pipeline for metagenomics and metatranscriptomics studies. BMC Bioinformatics. 2016;17:420.

UniProt Consortium. UniProt: a hub for protein information. Nucleic Acids Res. 2015;43(Database issue):D204-212.

Antunes KH, Singanayagam A, Williams L, Faiez TS, Farias A, Jackson MM, et al. Airway-delivered short-chain fatty acid acetate boosts antiviral immunity during rhinovirus infection. J Allergy Clin Immunol. 2023;151:447-457.e5.

Niu J, Cui M, Yang X, Li J, Yao Y, Guo Q, et al. Microbiota-derived acetate enhances host antiviral response via NLRP3. Nat Commun. 2023;14:642.

Martin FJ, Gomez MI, Wetzel DM, Memmi G, O’Seaghdha M, Soong G, et al. Staphylococcus aureus activates type I IFN signaling in mice and humans through the Xr repeated sequences of protein A. J Clin Invest. 2009;119:1931–9.

Stefan KL, Kim MV, Iwasaki A, Kasper DL. Commensal Microbiota Modulation of Natural Resistance to Virus Infection. Cell. 2020;183:1312-1324.e10.

Peignier A, Planet PJ, Parker D. Differential Induction of Type I and III Interferons by Staphylococcus aureus. Infect Immun. 2020;88:e00352-e420.

Steed AL, Christophi GP, Kaiko GE, Sun L, Goodwin VM, Jain U, et al. The microbial metabolite desaminotyrosine protects from influenza through type I interferon. Science. 2017;357:498–502.

Liu B, Zheng D, Zhou S, Chen L, Yang J. VFDB 2022: a general classification scheme for bacterial virulence factors. Nucleic Acids Res. 2022;50:D912–7.

Sastalla I, Monack DM, Kubatzky KF. Editorial: Bacterial Exotoxins: How Bacteria Fight the Immune System. Front Immunol. 2016;7:300.

Blin K, Shaw S, Augustijn HE, Reitz ZL, Biermann F, Alanjary M, et al. antiSMASH 7.0: new and improved predictions for detection, regulation, chemical structures and visualisation. Nucleic Acids Res. 2023;51:W46-50.

Morey JR, Kehl-Fie TE. Bioinformatic Mapping of Opine-Like Zincophore Biosynthesis in Bacteria. mSystems. 2020;5. https://doi.org/10.1128/msystems.00554-20.

Mishra AK, Choi J, Choi S-J, Baek K-H. Cyclodipeptides: An Overview of Their Biosynthesis and Biological Activity. Mol J Synth Chem Nat Prod Chem. 2017;22:1796.

Cheng C, Hua ZC. Lasso Peptides: Heterologous Production and Potential Medical Application. Front Bioeng Biotechnol. 2020;8:571165.

Kawashima T, Kosaka A, Yan H, Guo Z, Uchiyama R, Fukui R, et al. Double-Stranded RNA of Intestinal Commensal but Not Pathogenic Bacteria Triggers Production of Protective Interferon-β. Immunity. 2013;38:1187–97.

Yang X-L, Wang G, Xie J-Y, Li H, Chen S-X, Liu W, et al. The Intestinal Microbiome Primes Host Innate Immunity against Enteric Virus Systemic Infection through Type I Interferon. mBio. 2021;12. https://doi.org/10.1128/mbio.00366-21.

Gutierrez-Merino J, Isla B, Combes T, Martinez-Estrada F, Maluquer De Motes C. Beneficial bacteria activate type-I interferon production via the intracellular cytosolic sensors STING and MAVS. Gut Microbes. 2020;11:771–88.

Parker D, Prince A. Staphylococcus aureus induces type I IFN signaling in dendritic cells via TLR9. J Immunol Baltim Md. 1950;2012(189):4040–6.

Winkler ES, Shrihari S, Hykes BL, Handley SA, Andhey PS, Huang YJS, et al. The Intestinal Microbiome Restricts Alphavirus Infection and Dissemination through a Bile Acid-Type I IFN Signaling Axis. Cell. 2020;182:901-918.e18.

Parker D, Martin FJ, Soong G, Harfenist BS, Aguilar JL, Ratner AJ, et al. Streptococcus pneumoniae DNA initiates type I interferon signaling in the respiratory tract. mBio. 2011;2:e00016-00011.

Parker D, Cohen TS, Alhede M, Harfenist BS, Martin FJ, Prince A. Induction of type I interferon signaling by Pseudomonas aeruginosa is diminished in cystic fibrosis epithelial cells. Am J Respir Cell Mol Biol. 2012;46:6–13.

Annunziato F, Romagnani C, Romagnani S. The 3 major types of innate and adaptive cell-mediated effector immunity. J Allergy Clin Immunol. 2015;135:626–35.

Voigt AY, Emiola A, Johnson JS, Fleming ES, Nguyen H, Zhou W, et al. Skin microbiome variation with cancer progression in human cutaneous squamous cell carcinoma. J Invest Dermatol. 2022;142:2773-2782.e16.

Blanco-Míguez A, Beghini F, Cumbo F, McIver LJ, Thompson KN, Zolfo M, et al. Extending and improving metagenomic taxonomic profiling with uncharacterized species using MetaPhlAn 4. Nat Biotechnol. 2023:1–12.

Arumugam M, Raes J, Pelletier E, Le Paslier D, Yamada T, Mende DR, et al. Enterotypes of the human gut microbiome. Nature. 2011;473:174–80.

Köser CU, Ellington MJ, Peacock SJ. Whole-genome sequencing to control antimicrobial resistance. Trends Genet. 2014;30:401–7.

Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinforma Oxf Engl. 2011;27:863–4.

Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

Prjibelski A, Antipov D, Meleshko D, Lapidus A, Korobeynikov A. Using SPAdes De Novo Assembler. Curr Protoc Bioinforma. 2020;70:e102.

Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 2013;29:1072–5.

Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH. GTDB-Tk v2: memory friendly classification with the genome taxonomy database. Bioinformatics. 2022;38:5315–6.

Parks DH, Chuvochina M, Rinke C, Mussig AJ, Chaumeil P-A, Hugenholtz P. GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res. 2022;50:D785–94.

Xu S, Li L, Luo X, Chen M, Tang W, Zhan L, et al. Ggtree: A serialized data object for visualization of a phylogenetic tree and annotation data. Meta. 2022;1:e56.

Wickham H. Reshaping Data with the reshape Package. J Stat Softw. 2007;21:1–20.

Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2016.

Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, et al. Welcome to the Tidyverse. J Open Source Softw. 2019;4:1686.

Kassambara A. ggpubr:`ggplot2` Based Publication Ready Plots. R package version 0.3.0. 2020.

Wickham H, François R, Henry L, Müller K, Vaughan D. dplyr: a grammar of data manipulation. 2023.

Wickham H. The Split-Apply-Combine Strategy for Data Analysis. J Stat Softw. 2011;40:1–29.

Linlin Y. ggvenn: Draw venn diagram by “ggplot2.” 2023.

Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinforma Oxf Engl. 2014;30:2068–9.

Palmer JM, Stajich J. Eukaryotic genome annotation (v1.8.1). Zenodo. 2020. https://doi.org/10.5281/zenodo.4054262.

Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.

Hyatt D, Chen G-L, Locascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119.

Kassambara A. rstatix: Pipe-friendly framework for basic statistical tests. 2023.

Neuwirth E. RColorBrewer: ColorBrewer Palettes. 2022.

Kolde R. pheatmap: Pretty Heatmaps. CRAN. 2019.

Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinforma Oxf Engl. 2013;29:15–21.

Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinforma Oxf Engl. 2014;30:923–30.

Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. 2014;32:896–902.

Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102:15545–50.

Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003;34:267–73.

Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation. 2021;2:100141.

Guangchuang Y, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS. 2012;16:284–7.

Schoggins JW. Interferon-Stimulated Genes: What Do They All Do? Annu Rev Virol. 2019;6:567–84.

Waskom ML. seaborn: statistical data visualization. J Open Source Softw. 2021;6:3021.

Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine Learning in Python. J Mach Learn Res. 2011;12:2825–30.

Charlier F, Weber M, Izak D, Harkin E, Magnus M, Lalli J, et al. Statannotations (v0.5). Zenodo. 2022.

The pandas development team. pandas-dev/pandas: Pandas(v2.1.0). Zenodo. 2023.

Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D, et al. Array programming with NumPy. Nature. 2020;585:357–62.

Hunter JD. Matplotlib: A 2D graphics environment. Comput Sci Eng. 2007;9:90–5.

Gu Z. Complex heatmap visualization. iMeta. 2022;1:e43.

Wickham H. forcats: tools for working with categorical variables (factors). 2023.