We recruited 110 participants, including 20 healthy controls, 30 subjects with AR, 30 subjects with asthma and 30 subjects with asthma and comorbid AR. Subjects with asthma were diagnosed based on GINA (Global Initiative for Asthma) clinical criteria. Subjects diagnosed with AR or asthma and comorbid AR were confirmed using the Allergic Rhinitis and its Impact on Asthma (ARIA) guideline (Bousquet et al., 2019). Exclusion criteria were as follows: (1) age less than 18 or greater than 80 years; (2) respiratory infection within 4 weeks; (3) antibiotic use within 3 months; (4) systemic steroid therapy within 5 months; (5) pulmonary disease other than asthma; (5) cancer; (6) currently pregnant or lactating.

Smoking status, comorbidities, relevant medical history, medications, first degree relatives with a history of AR or asthma, and type of allergen were recorded for all participants. The severity of AR was quantified using the total nasal symptom score (TNSS). The Asthma Quality of Life Questionnaire (AQLQ) was used to assess the life quality for patients with asthma (Juniper et al., 1999a). The Asthma Control Questionnaire (ACQ-7) was used to evaluate asthma control status (Barnes et al., 2014). Controlled or partially controlled asthma was defined as ACQ less than 1.5, uncontrolled asthma was defined as ACQ more than 1.5 (Juniper et al., 1999b).

This study was conducted at Ningbo First Hospital with approval from the ethics committee of Ningbo First Hospital (approval 2020-R145). All subjects provided written informed consent before participating in the study.

All participants performed pulmonary function testing according to the American Thoracic Society (ATS) guidelines (Graham et al., 2019). All measurements were undertaken at least three times using the same spirometer (Jager, MasterScreen, Höchberg, Germany) (Hankinson et al., 1999). Pulmonary function was recorded as a percentage of predicted forced expiratory volume in one second (FEV₁), forced vital capacity (FVC), maximal mid-expiratory flow (MMEF75/25), and peak expiratory flow (PEF).

Nasal lavage was performed according to previously described methods (Hirvonen et al., 1999; Frischer et al., 2000). Participants’ heads were briefly held downwards to avoid fluid entering the nasopharynx and allow the liquid to drip into a sterile basin. The nasal passage was slowly instilled with 10 ml of sterile saline, and the nasal lavage was immediately recovered into 15 ml in sterile tubes and placed on ice. The tubes were frozen and stored at -80°C until analysis.

Total DNA was extracted from 5ml of nasal lavage fluid using DNeasy PowerWater Kit (Qiagen, Hilden, Germany) following instructions provided by the manufacturer. DNA yield and integrity were identified using a Qubit Fluorometer (Qubit, Invitrogen, USA) and 1% agarose gel electrophoresis, respectively. The V3-V4 regions in bacteria 16S rRNA gene of the nasal lavage fluid were amplified with PCR primers 341F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’- GGACTACHVGGGTWTCTAAT-3’). The primers were removed before processing. Subsequently, the PCR amplicons were purified with Agencourt AMPure XP magnetic beads, dissolved in Elution Buffer (pH = 8.0) and labelled. The insert fragment size was estimated with an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA), then the sample was sequenced on an Illumina HiSeq2500 platform (Illumina, Inc., San Diego, CA, USA) using the PE300 module. Low quality and ambiguous bases were removed using Cutadapt (V2.6) according to the following: (a) raw reads with average quality less than 20, (b) final length less than 75% of their original sequence, (c) reads with an ambiguous bases (N), (d) reads with low complexity (repeats of length more than 10 bases). Amplicon sequence information in the study is available at the Sequence Read Archive (SRA) under BioProject Accession Number PRJNA793600.

The overlap and paired-end reads were processed by the Fast Length Adjustment of Short reads program (FLASH, v1.2.11) (Magoč and Salzberg, 2011) to obtain sequence data. The amount of raw sequence data was more than 50,000 for each sample. Next, sequences were clustered into operational taxonomic units (OTUs) based on a 97% similarity threshold using UPARSE software (v7.0.1090) (Edgar, 2013). OTU representative sequences were classified by Ribosomal Database Project (RDP) Classifier (RDP release version 11.5; release date, 2017.12.0); 0.8 was set as the minimum confidence threshold (Cole et al., 2014). Singleton OTUs and low abundant OTUs were not included in further analyses, additionally, any OTUs not classified as bacteria were removed. Sequences identified as mitochondria were also eliminated. All samples were processed simultaneously in the same research laboratory (BGI, Shenzhen, China) to control batch variation.

R (v3.2.1) was applied to determine α-diversity (the evenness and richness of bacteria taxonomic diversity), which was expressed by the Chao, Shannon, and Simpson index. β-diversity (distance between samples, based on the difference in OUT in each sample) was performed by Quantitative Insights into Microbial Ecology (QIME, v1.80) (Caporaso et al., 2010); it was evaluated by the principal coordinates analysis (PCoA) and the Partial least-squares discrimination analysis (PLS-DA). PLS-DA was conducted with R package mixOmics (Rohart et al., 2017). False discovery rate (FDR) was calculated according to the Benjamini-Hochberg correction. The heatmap was plotted with Graph Pad Prism 8 (Graph Pad Software, USA). Linear discriminant analysis of effect size (LEfSe) combined with linear discriminant analysis (LDA) was adopted to explore the differences in taxonomic composition (Segata et al., 2011); LDA > 3.6 was presented in this study (default: 2.0).

The function of the nasal microbiome was analyzed via Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) and further categorized into Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (Langille et al., 2013). The associations between variables of bacteria load and metabolic function were evaluated by Pearson correlation coefficient (r).

Basic characteristics were performed using SPSS software (version 21; SPSS, Inc., Chicago, IL, USA). Categorical variables were compared using the chi-squared test or the Fisher exact test. A two-tailed P value less than 0.05 was considered statistically significant.

A total of 110 participants were enrolled, including 90 patients with disease and 20 healthy controls ( Figure 1 ). The subjects with disease were divided into three groups, 30 subjects with AR, 30 subjects with asthma and 30 subjects with asthma and comorbid AR. Detailed demographic and baseline characteristics are listed in Table 1 . No significant differences in age, BMI, sex or smoking status were observed among the groups. Pulmonary function tests were performed for all participants. Family history and allergic of study were listed in Supplementary Table 1 . Participants The subjects with asthma and those with asthma and comorbid AR were subdivided into groups according to their ACQ score: controlled or partially controlled (defined by ACQ score < 1.5) and uncontrolled (defined by ACQ score > 1.5). The characteristics of the asthmatic patients with or without comorbid AR are summarized in Supplementary Tables 2, 3 .

After filtering for low-quality reads, the average number of clean reads was 61140.73 ± 1538.64 (SD) for subjects for AR, 57586.10 ± 4623.94 for asthma, 55341.50 ± 5747.72 for combined asthma + AR, and 59230.30 ± 6256.84 in the healthy controls, detailed data and species accumulation curves was listed in Supplement File 1 .

The researchers used three metrics to compare the α-diversity: the Chao, Shannon, and Simpson index. As shown in Figure 2A , the Chao index tended to be lower in subjects with AR (503.38 ± 101.70, P = 0.001), asthma (511.01 ± 106.75, P = 0.001), and combined asthma + AR (504.99 ± 121.14, P = 0.001) when compared with healthy controls (632.26 ± 141.98, Figure 2A ). Similarly, a lower Shannon index was found in subjects with asthma (3.73 ± 0.68, P = 0.013) and combined asthma + AR (3.61 ± 0.75, P = 0.004) in comparison to the control subjects (4.21 ± 0l58, Figure 2A ). However, no significant difference in either the Shannon or Simpson index of nasal microbiome communities was found between healthy controls and patients with AR ( Figure 2A ). Asthma may have a stronger effect on α-diversity than AR in the upper airway.

PLS-DA analyses was applied to study the structure of nasal microbiome communities (beta-diversity) among the four groups. Furthermore, PLS-DA showed distinct microbiome composition among all groups ( Figure 2B ). Overall, the results exhibit the specific structure in nasal bacterial community composition among subjects with AR, asthma, combined asthma + AR and healthy controls.

Following initial evaluations for differences in overall nasal bacterial community composition, we investigated the relative abundance of specific taxa in subjects with disease and healthy controls. Evaluation of the general landscape of the nasal microbiome revealed similar bacterial communities at phylum and genus levels ( Figure 2C ). We then evaluated the relative abundance of the bacterial communities at phylum and genus levels within the four groups. Firmicutes were the most predominant phylum in all those four groups, followed by Bacteroidetes, Proteobacteria, Actinobacteria ( Figure S1 ). Compared with healthy controls, Firmicutes were enriched in the samples of subjects with AR (mean relative abundances, 53.55% vs 68.73%, FDR adj P-value = 0.005) but depleted in subjects with asthma with or without AR (asthma, 53.55% vs 45.37%, FDR adj P-value < 0.001; Asthma + AR, 53.55% vs 41.97%; FDR adj P-value < 0.001). No significant difference was found in the relative abundance of Actinobacteria or Proteobacteria among the four groups with disease ( Figure S2 and Supplementary Table 4 ).

To further confirm these findings, we used the LEfSe analysis (Segata et al., 2011); this found marked differences in the nasal bacteria community among the four groups ( Figure S3 ). Subjects with asthma had a high proportion of Actinomycetaceae, Listeriaceae, and Neisseriaceae at the family level. The samples of subjects with combined asthma + AR were enriched with Provotellaceae, Sphingobacteriaceae, Rhodocyclaceae, Aeromonadaceae, and Leptotrichiaceae (LAD > 2.0, P < 0.05, Figure S3A ). At the genus level, the top 10 most abundant nasal microbiome were selected for further analysis. We identified seven taxa that differed among the four groups: Streptococcus, Prevotella, Faecalibacterium, Neisseria, Lactobacillus, Haemophilus and Clostridium_XlVa (FDR adj P-value < 0.05 in each analysis, Figure S4 and Supplementary Table 5 ). The relative abundance of Streptococcus was higher in the asthma (P < 0.01) and combined asthma + AR (P < 0.01) groups compared with healthy controls, respectively. However, the relative abundance of Faecalibacterium, Lactobacillus, and Clostridium_XlVa were lower in the asthma and combined asthma + AR groups in relation to healthy controls (P < 0.05).

LDA analyses were adopted to determine the specific bacterial taxa in each group (LDA > 3.6, P < 0.05; Figure S3B ). The genera of Prevotella and Kineococcus characterized the samples of the healthy control subjects. Additionally, 13 genera of the nasal microbiome were differentially abundant in subjects with AR, including Faecalibacterium, Lactobacillus, Escherichia and Clostridium_IV genera; 84.6% (11/13) belong to the Firmicutes phylum ( Figure S3B ). The genera of Neisseria and Rothia dominated the microbiome of subjects with asthma; whereas, the genera of Pelomonas, Alloprevotella, Leptotrichia, and Granulicatella dominated the microbiome of subjects with asthma and comorbid AR ( Figure S3B ). We used heatmaps to visualize the relative abundance of the dominant taxonomic communities at the genera level of the nasal microbiome in different groups. A higher quantity corresponds to a deeper red color. The heatmap displayed a similar abundance of bacteria in subjects with asthma relative to those with asthma and comorbid AR ( Figure S5 ).

At the genus level, there was no significant difference in taxonomic distribution between the subjects with asthma and those with asthma and comorbid AR ( Figure 3 ). However, a significant difference was found in the genera composition between the samples of subjects with AR, asthma, combined asthma + AR and healthy controls ( Figure 3 ). The relative abundance [median (interquartile range, IQR)] of Prevotella was lower in the AR group compared with healthy controls [0.78% (1.99%) vs. 2.23% (6.58%), FDR adj P-value = 0.016, Figure 3A ]. Conversely, the relative abundance of Faecalibacterium was greater in the AR group [8.67% (9.71%)] relative to healthy controls [3.91% (6.25%), FDR adj P-value = 0.011], subjects with asthma [0.50% (4.20%), FDR adj P-value < 0.001], and subjects with asthma and comorbid AR [0.24% (2.66%), FDR adj P-value < 0.001, Figure 3B ). Similarly, the genera of Lactobacillus, Escherichia, Clostridium_IV, Blautia and Butyricicoccus in the AR group had the highest relative abundance of all groups ( Figure 3C-G ). In contrast, a high level of Neisseria, Rothiawas, Pelomonas, and Alloprevotella was detected in subjects with asthma and combined asthma + AR when compared to the healthy controls ( Figure 3I-L ). However, the most of these genera have very low relative abundance across all samples.

We further assessed whether the predicted function of the nasal bacterial dominant taxa in different groups differed distribution across the groups. In total, 161 metabolites were analyzed by KEGG pathway analysis. The Pearson correlation heat map in Figure 4 summarized the clustering of enriched functions among the three groups with disease and the healthy control group. Nasal bacterial communities enriched in AR had multiple functional features related to the pathway of pentose phosphate and galactose metabolism. Whereas, decreased activities in the pathway of primary bile acid biosynthesis were observed in the genus of Rothia (Pearson r = - 0.48; P <0.001, Figure 4 ) and Neisseria (r = - 0.44; P <0.001, Figure 4 ), which enriched in subjects with asthma ( Figure 4 ). The genera of Pelomonas were highly abundant in patients with combined asthma + AR; it was involved in the pathway of lipid metabolism, especially with upregulated fatty acid degradation (Pearson r = 0.62; P <0.01, Figure 4 ). Additionally, the genera of Pelomonas associated with the metabolism of glutathione (r = 0.726; P <0.05, Figure 4 ).

To further investigate the relationship between nasal microbiome composition and disease control, we divided subjects with asthma into two groups: 23 (76.7%) of them with partially or well-controlled, and 7 (23.3%) subjects with uncontrolled asthma. No significant differences in microbiome evenness or diversities were found according to the Shannon, Chao, and Simpson indexes ( Figure 5A ). β-diversity, based on the PLS-DA analyses, was performed to compare microbial structure. It is notable that there was a clear separation between the two groups from the PLS-DA analysis ( Figure 5B ).

The general features of the nasal microbiome at phylum and genus levels are shown in Figure 5C . Five dominant phyla were found in nasal lavage fluid samples: Firmicutes, Proteobacteria, Actinobacteria, Bacteroidetes and Fusobacteria ( Supplementary Figure 6 ). Similarity in the abundance at phylum and genus levels was observed between the two groups ( Supplementary Table 6 ). However, a higher proportion of Cyanobacteria was found in subjects with uncontrolled asthma, compared to subjects with controlled or partially controlled asthma ( Supplementary Table 6 ). In contrast, a lower proportion of Firmicutes was observed in subjects with uncontrolled asthma ( Supplementary Table 6 ). The LDA analyses were performed to validate the specific genera microbiome according to level of asthma control (LDA > 2.0, P < 0.05, Figure 5D ). The three top genera of nasal taxa in patients with controlled or partially controlled asthma were Flavonifractor (phylum Firmicutes), Pseudoflavonifractor (phylum Firmicutes) and Anaerofilum (phylum Firmicutes). While Pedobacter (phylum Bacteroidetes), Jeotgalicoccus (phylum Firmicutes), and Janthinobacterium (phylum Proteobacteria) were enriched in the samples of the uncontrolled asthma group ( Figure 5D ).

It was speculated that the bacterial genus enrichment in the upper airway might lead to metabolic changes in subjects with asthma. Here, metabolomics analysis was conducted on nasal bacteria with respect to the enrichment in different degrees of asthma severity. The data presented above suggest that patients with uncontrolled asthma are enriched with the genus of Janthinobacterium and Pedobacter. Conversely, the genus of Collinsella, Anaerofilum, Novosphingobium, Flavonifractor, Pseudoflavonifractor, and Vibrionimonas were enriched in patients with controlled or partially controlled asthma. Both the Lipoic acid metabolism (r = 0.47, P =0.001) and Vibrio cholerae infection (r = 0.40, P = 0.004) pathways exhibited the positive correlation with the bacterial loading of Janthinobacterium in subjects with uncontrolled asthma ( Figure S7 ). Interestingly, a negative correlation was observed between the pathway of Lipoic acid metabolism and the genera of Flavonifractor (r = - 0.53, P < 0.001), which was enriched in patients with controlled and partially controlled asthma.

Among 9 (30%) subjects with uncontrolled asthma and combined asthma + AR, a lower Shannon index (P = 0.040) was observed in nasal microbiome composition. This suggests reduced evenness and diversity of nasal bacteria relative to those 21 (70%) subjects with controlled or partially controlled combined asthma + AR. A similar trend was observed in the Chao and Simpson indexes, although the difference was not statistically significant ( Figure 6A ). To determine variation between different asthma control groups in participants with asthma and comorbid AR, PLS-DA analyses were performed. When compared with the uncontrolled asthma and comorbid AR groups, significant separation of PLS-DA was also found between these two groups ( Figure 6B ).

We next assessed the landscape of the nasal microbiome at the phylum and genus levels in all subjects with asthma and comorbid AR ( Figure 6C ). Proteobacteria were enriched in the samples of subjects with uncontrolled asthma and comorbid AR; however, the difference was not statistically significant after being adjusted for FDR ( Figure S8 and Supplementary Table 9 ). In nasal lavage samples, the relative abundance of five genera: Pseudoflavonifractor (phylum Firmicutes), Eubacterium (phylum Firmicutes), Sporobacter (phylum Firmicutes), Intestinimonas (phylum Firmicutes), and Mycobacterium (phylum Actinobacteria) were enriched in the samples of subjects with controlled or partially controlled in asthma and comorbid AR group ( Figure 6D ).

We further conducted function analysis to explore nasal bacteria in the development of allergic inflammation. The genera of Eubacterium, Pseudoflavonifractor, and Sporobacter were enriched in patients with controlled or partially controlled combined asthma + AR; they were adapted to assess the correlation between taxa loading and metabolic pathway. The load of Eubacterium is closely related to the active pathway of NOD-like receptor signaling (r =0.63, P < 0.001) and RNA polymerase (r=0.66, P < 0.001) ( Figure S9 ).