For this study, participants were recruited from the COPSAC clinic. All participants and their parents are part of the COPSAC₂₀₀₀ (Bisgaard, 2004) or COPSAC₂₀₁₀ (Bisgaard et al., 2013) cohorts, and have provided written informed consent. A total of 16 patients were included. These include 8 healthy patients, 4 patients with asthma, and 4 patients with PCR-validated COVID-19 infection. The diagnosis of asthma was physician-diagnosed and performed solely in the research unit according to a previously validated algorithm (Bisgaard et al., 2011). Participants were invited to acute care visits in the research unit during the COVID-19 pandemic, whenever they had airway infections. Here, nasopharyngeal aspirates were obtained using a soft suction catheter passed through the nose into the nasopharynx, following established procedures (Bisgaard et al., 2007). The aspirates were diluted in 1 mL of sterile 0.9% NaCl solution, transported to the laboratory, and stored at −80 °C prior to DNA extraction.

The DNA extraction procedure was conducted using MasterPure™ Complete DNA and RNA Purification Kit (Epicentre), following the manufacturer’s guidelines. Control reagents were extracted and used as a negative control. The extracted DNA was stored at −20 °C until further analysis through shotgun metagenome sequencing. The DNA samples were shipped to Novogene (Cambridge, UK). Novogene conducted DNA library preparation and sequencing using an Illumina NovaSeq 6,000 system, with 2 × 150 bp paired-end sequencing. We did not detect any DNA from the negative control. Moreover, the negative control was determined to have insufficient concentration for library preparation and was subsequently excluded from further processing.

The process of quality filtering and trimming of the raw sequence reads was carried out using KneadData v0.12.0.¹ Trimmomatic v0.39 (Bolger et al., 2014) was used to eliminate Illumina sequencing adaptors and conduct initial quality filtering on raw shotgun metagenomic reads by excluding low-quality reads with a Phred score below 20. Subsequently, the reads were aligned against the reference genome (hg37) using Bowtie2 v.2.4.5 (Langmead and Salzberg, 2012) with default parameters to eliminate human host genome contamination. The filtered reads were then used for further analyses.

To analyze the structure and diversity of bacterial communities, we use Kraken2 v2.1.2 with the Kraken2 PlusPFP database containing genome references from archaeal, bacterial, viral, protozoa, fungi, plant, and human sources.² This software is used to classify shotgun metagenomic reads (Wood et al., 2019). Following this classification, Bracken v2.6.0 (Lu et al., 2017) was used to estimate species abundances. Subsequently, the resulting bacterial and fungal abundance tables were collapsed to genus-level read counts.

We used the ARGs-OAP (v3.0) (Yin et al., 2023) to identify and quantify antibiotic resistance genes (ARGs). Clean reads were aligned with the structured ARG database (SARG v3.0) using BLAST+ (v2.12.0) (Camacho et al., 2009), applying the criteria of a similarity threshold of 80%, an e-value of ≤1e−7, and a query coverage of at least 75%. The abundance of ARGs was expressed as the ratio of ARG copies to copies of the 16S rRNA gene, as determined by the ARGs-OAP framework.

We employed a variety of binning techniques, such as Maxbin2 v2.2.7, MetaBAT2 v2.12.1, and CONCOCT v1.1.0, (Alneberg et al., 2014; Kang et al., 2019; Wu et al., 2016) to generate metagenome-assembled genomes (MAGs). DASTool v1.1.1 was utilized to identify MAGs with the highest quality from all binning tools (Sieber et al., 2018). The quality assessment of MAGs, regarding completeness and contamination levels, was conducted utilizing CheckM v1.0.13 (Parks et al., 2015). Only medium-quality MAGs with completeness of more than 50% and contamination levels below 10%, were retained for comparison of metabolic capabilities (Bowers et al., 2017). Dereplication of metagenome-assembled bacterial genomes was done using dRep v2.2.3 (Olm et al., 2017), resulting in a nonredundant set of genomes. Taxonomical information for each MAG was obtained through the Genome Taxonomy Database Toolkit (GTDB-Tk), and phylogenetic trees were constructed using PhyloPhlAn (Asnicar et al., 2020; Chaumeil et al., 2020). The generated MAGs were subsequently analyzed using ABRicate,³ to align the contigs with the Comprehensive Antibiotic Resistance Database (CARD) (Alcock et al., 2020) and the Virulence Factor Database (VFDB) (Liu B. et al., 2022) and to detect the presence of ARGs and virulence factors. We defined the presence of antibiotic resistance genes and virulence factors in the genome if the sequences have at least 70% similarity and 70% gene coverage. Abundance profiles of each MAG were calculated using CoverM v0.4.0⁴ with the option -m rpkm, to determine MAG abundance as mapped reads per kilobase per million reads (RPKM).

Statistical analysis and graphical representation were performed in R v4.1.2 using the R packages Phyloseq v1.38.0 and vegan v2.6.4 (McMurdie and Holmes, 2013; Oksanen et al., 2007). Microbial alpha and beta diversity analyses were conducted using normalized data obtained through subsampling to the minimum number of reads. An analysis of variance (ANOVA) was used to identify significant differences (p < 0.05) in the diversity of microbial communities and antimicrobial resistance genes (ARGs). Group comparisons were made using a t-test for multiple comparisons with Benjamini–Hochberg false discovery rate correction (FDR adjusted p < 0.05). For beta diversity analysis, Bray–Curtis dissimilarity matrices were constructed from normalized datasets. The betadisper function in the vegan package v. 2.6.4 was employed to calculate the distance of each sample from the corresponding group centroid based on the Bray–Curtis distance matrix. These matrix distances were subsequently subjected to permutational multivariate analysis of variance using the Adonis2 function and Bonferroni multiple testing correction (p < 0.05) to evaluate the significant effects of disease state on the microbial community and resistome compositions. Lastly, biomarkers at the bacterial genus level using a linear discriminant analysis effect size (LefSe) analysis (Segata et al., 2011) for both healthy individuals and those with asthma. We defined differentially abundant taxa as those with an LDA score greater than 3 and an FDR-adjusted p value less than 0.1.

Shotgun metagenomic sequencing was conducted on nasopharyngeal aspirate samples collected from 16 participants (10 males and 6 females), including individuals with asthma, COVID-19 infection, and healthy controls (see Table 1). Among these, 11 participants were from the COPSAC₂₀₀₀ cohort and 5 from the COPSAC₂₀₁₀ cohort. Samples were collected at acute clinical visits arranged due to symptoms suggestive of COVID-19 infection during the first year of the pandemic. Participants’ ages ranged from 9 to 19 years at the time of sample collection. The average age of individuals with asthma was 10 years, those with COVID-19 infection averaged 12 years, and healthy controls had an average age of 15 years.

A total of 1,167,675,803 metagenome reads (min: 46,205,750; max: 95,141,275; average: 72,979,737) were generated and analyzed from all samples (Supplementary Table S1). The analysis of the metagenomes using Kraken2 and Bracken indicated that, on average, 99.1% (range: 97.4–99.6%) of the reads were categorized as human reads, while bacterial and fungal reads accounted for, on average, 0.46% (range: 0.002–1.883%) and 0.015% (range: 0.0003–0.148%), respectively. A total of 7,558 bacterial taxa and 87 fungal taxa were identified from the nasopharyngeal aspirate samples. Our analysis indicates that, although the shotgun metagenome reads are relatively shallow, the rarefaction curve suggests that the current sequencing depth is sufficient to capture the overall diversity present in the samples (Supplementary Figure S1). The removal of low-abundance taxa has been recommended to minimize the risk of false positives in shotgun metagenomic data (Edwin et al., 2024; Garrido-Sanz et al., 2022). We therefore excluded rare bacterial genera (relative abundance <0.01%) and used the resulting data for further analysis, which identified 332 bacterial species and 50 fungal species.

The dominant bacteria in the nasopharyngeal aspirates samples belonged to Actinomycetia (mean of relative abundance 33.8%), Gammaproteobacteria (26.2%), Alphaproteobacteria (14.4%), and Bacilli (9.1%) whereas the dominant fungi in the samples were Malasseziomycetes (42.9%), Saccharomycetes (35.3%), and Sordariomycetes (11.6%). At genus level, several bacterial genera—including Cutibacterium, Pseudomonas, Moraxella, Paracoccus, and Prevotella—and fungal genera—including Aspergillus, Fusarium, Malassezia, and Saccharomyces —were identified in at least 80% of the samples. These taxa accounted for an average of 46.7% of the total bacterial reads and 87.4% of the total fungal reads. At a more detailed taxonomic level, the bacterial community profiles were primarily dominated by several species, including Cutibacterium acnes, Pseudomonas brenneri, Moraxella osloensis, Prevotella nigrescens, Corynebacterium tuberculostearicum, Mycobacterium canetti, and Staphylococcus epidermidis. These species were consistently detected across all samples. Healthy individuals were dominated by bacterial genera Cutibacterium (19.2%), Moraxella (10.7%), Paracoccus (9.6%), Pseudomonas (6.4%), Corynebacterium (3.7%), and Staphylococcus (3.6%). A similar profile was also observed in the individuals with COVID-19. While the bacterial community composition among healthy individuals was relatively similar; a high variation in composition was noted among individuals with asthma. For instance, Dolosigranulum, a bacterial taxon typically present in low abundance in healthy individuals, was found in significantly higher abundance in an individual with asthma (Figure 1A). LEfSe analysis identified several taxa as biomarkers that are enriched in healthy individuals relative to those with asthma (Supplementary Table S2). Specifically, Paracoccus, Moraxella and Brevundimonas, were significantly associated with the healthy cohort.

Differences in bacterial community richness and structure were observed between healthy individuals and individuals with asthma. A significantly lower number of bacterial species (S) (p = 0.038) and Shannon diversity index (H′) (p = 0.016) were observed in individuals with asthma (S = 80, H′ = 2.2) compared to healthy individuals (S = 133, H′ = 3.3; Figures 1B,D). In contrast, no significant difference was detected between the groups when diversity was estimated using the Chao1 index (Figure 1C). In contrast, no significant differences in alpha diversity were observed between the healthy individuals without asthma and individuals with COVID-19. Beta diversity analysis revealed that the disease state explained 20.5% of the bacterial community variation. Multiple pairwise comparisons revealed significantly different clustering of bacterial community structures between healthy individuals and those with asthma (R² = 0.204, P_adj = 0.039), but no significant differences were observed when comparing healthy individuals with those diagnosed with COVID-19 (Supplementary Figure S2). Furthermore, beta dispersion analysis indicated that the variability in bacterial community structure was greater in individuals with asthma compared to their healthy counterparts (Figure 1E).

Compared to bacterial community profiles, fungal community variation did not exhibit a clear distinction between healthy individuals and those with asthma or COVID-19. For example, the fungal species Malassezia restricta was the predominant member of the fungal community in a subset of individuals (n = 8), irrespective of their health status (Figure 2A). On the other hand, Saccharomyces was the predominant member of the fungal community in other individuals (n = 6). A higher number of fungal genera and an increased Shannon diversity index were observed in individuals with asthma compared to healthy individuals and those with COVID-19. However, statistical analysis on fungal data did not reveal any differences in fungal diversity. No differences were observed in the number of fungal genera (p = 0.105, Figure 2B), Chao1 index (p = 0.707, Figure 2C), and Shannon diversity index (p = 0.158, Figure 2D). Similarly, beta diversity and dispersion analysis did not indicate significant differences in fungal community structures (p = 0.972) among individuals with various disease states (healthy vs. asthma vs. COVID, Figure 2E). Our findings suggest that the composition of the airway fungal community is more individualized than significantly influenced by the disease state. Consequently, further analyses were concentrated on the bacterial community data.

From all samples in the dataset, antibiotic efflux (63.9%) and antibiotic inactivation (24.6%) were the primary mechanisms of resistance observed in the nasopharyngeal microbiome regardless of the disease state (Figure 3A). Collectively, the identified resistance determinants encompassed 23 distinct drug classes. The majority of the ARGs confer resistance to multidrug (45.5%), followed by those genes that confer resistance to aminoglycosides (15.9%), tetracyclines (14.8%), polymyxin (5.3%), beta-lactam (4.9%) and macrolide-lincosamide-streptogramin (3.8%). Furthermore, we assigned taxonomically each read that was annotated as ARGs with Kraken2 (Wood et al., 2019) and identified that a diverse range of taxa contributed to ARG diversity in healthy individuals (Figure 3B). Interestingly, based on the taxonomic classification of the short reads identified as antimicrobial resistance, a high proportion of ARGs (healthy—13.3%, asthma—21.1%, and COVID-19–15%, Figure 3B) were associated with various Pseudomonas species, including P. brenneri, P. putida, P. aeruginosa, and P. fluorescens. This suggests that this widespread environmental genus was the main contributor to the nasopharyngeal resistome.

We subsequently constructed and analyzed metagenome-assembled genomes (MAGs) to identify potential genes associated with the pathogenicity of taxa linked to disease states. From 16 metagenome samples, we reconstructed 19 bacterial MAGs with a minimum of 30% genome completeness. Additionally, the same bacterial species, such as Micrococcus luteus and Cutibacterium acnes, were detected across multiple sample (Supplementary Table S3). Of those MAGs, a total of 10 non-redundant MAGs were retained, with at least 50% completeness and less than 10% contamination (Figure 3C). The MAGs were classified within the following taxonomic groups: Pseudomonas brenneri, Pseudomonas sp., Prevotella nigrescens, Prevotella oris, Stenotrophomonas sp., Paracoccus sp., Cutibacterium acnes, Micrococcus luteus, Fusobacterium sp., and Parvimonas sp. Notably, 7 of the MAGs were categorized as high-quality, exhibiting at least 90% completeness and less than 5% contamination.

We identified several taxa that harbour multiple antibiotic resistance genes (ARGs) and virulence factors. Most of these taxa belong to the Gammaproteobacteria class, namely Pseudomonas and Stenotrophomonas. The majority of the ARGs detected in these taxa were antibiotic efflux. Pseudomonas MAGs that carried multiple ARGs were predominantly detected in individuals with asthma and COVID-19. We conducted further analysis of ARGs and virulence factors within these genomes namely Pseudomonas brenneri and Stenotrophomonas maltophilia (Supplementary Tables S4, S5). We observed that the MAG identified in this study exhibits profiles similar to those found in host-associated (plant and animal) samples, as opposed to environmental samples. This suggests a potential association with adaptability to the host environment. We also compared the ARG and virulence profile associated with Stenotrophomonas, the second taxon containing multiple ARGs and virulence factors, with those of clinical and environmental strains of Stenotrophomonas maltophilia, a species known to be an opportunistic pathogen (Supplementary Table S5; Lira et al., 2017). Our analysis revealed a comparatively lower number of detected ARGs and virulence factors in the MAGs identified in this study relative to S. maltophilia, suggesting a lower potential pathogenicity of the Stenotropomonas characterized herein. We also detected ARGs in MAGs belonging to Cutibacterium acnes and Micrococcus luteus; however, their numbers were significantly lower than those observed in the aforementioned taxa.

Some MAGs were exclusively detected nasopharyngeal in individuals with asthma. We identified two MAGs belonging to two different species of Prevotella: Prevotella nigrescens and Prevotella oris. Notably, these MAGs were exclusively found in individuals with asthma (Figure 2C). Furthermore, one of the species namely P. nigrescens may exhibit increased virulence due to the presence of the tet37 gene, which confers resistance to tetracycline, along with the virulence factors gmd and neuB. Virulence factors were also detected from MAGs belonging to Pseudomonas and Stenotrophomonas. Genes encode different structural parts of the flagellar apparatus, i.e., flgC, flgF, flgG, flgI, and flgG, and twitching motility phenotype, i.e., pilG, pilH, and pilJ were detected from the MAGs.

The number of detected ARGs in the nasopharyngeal metagenome samples was different between individuals with different disease states. In our study, nasopharyngeal samples of healthy individuals tend to have a higher number of ARGs in comparison to individuals with asthma and COVID-19 (p = 0.093, Supplementary Figure S3A). ARG composition did not show significant variation among individuals with different disease states (p = 0.466, Supplementary Figure S3B). The Mantel test indicated a moderate correlation between bacterial composition and the resistome profiles (Mantel test—p = 0.016, r = 0.487). There were no statistical differences (p = 0.387) in ARG load normalized with copy of 16S rRNA genes between different individuals (Supplementary Figure S3C). Consistent results were obtained using RPKM and PPM normalization, with no significant differences in ARG community composition between the tested groups (RPKM—p = 0.617; PPM—p = 0.677). However, two samples that belong to the asthma group had a relatively lower abundance of ARGs (normalized abundance 0.09 and 0.04) whereas one sample from COVID-19 group had substantially higher ARG load (normalized abundance 5.6) in comparison to other samples (median of normalized abundance 0.56). Subsequent correlation analyses revealed a significant positive association between the number of bacterial taxa belonging to Pseudomonadales and the overall ARG load (Pearson’s correlation coefficient—p = 0.005, r = 0.66), as well as multidrug resistance genes (p < 0.001, r = 0.78). As mentioned above, our study highlights the significant role of Pseudomonadales as the primary contributor to the nasopharyngeal resistome.

We specifically focused on the detection and analysis of genes conferring macrolide resistance. The presence of high-abundance genes conferring resistance to macrolides was previously reported in airway and stool samples from individuals with asthma (Ivan et al., 2024; Wilson et al., 2024) who are often prescribed these medications. However, in our study, this ARG class did not show significant differences among the different disease groups (p = 0.636). The aforementioned results indicated that a specific subset of the taxa contributes to the ARG abundance in the nasopharynx. These findings led us to further investigate the correlation between ARG abundance and the diversity and abundance of specific taxa. Subsequent correlation analyses revealed a significant positive association between the number of bacterial taxa belonging to Pseudomonadales and the overall ARG load (Pearson’s correlation coefficient—p = 0.005, r = 0.66), as well as multidrug resistance genes (p < 0.001, r = 0.78). As mentioned above, our study highlights the significant role of Pseudomonadales as the primary contributor to the nasopharyngeal resistome.