Particulate Matter Exposure and Allergic Rhinitis: The Role of Plasmatic Extracellular Vesicles and Bacterial Nasal Microbiome

Particulate matter (PM) exposure is linked to the worsening of respiratory conditions, including allergic rhinitis (AR), as it can trigger nasal and systemic inflammation. To unveil the underlying molecular mechanisms, we investigated the effects of PM exposure on the release of plasmatic extracellular vesicles (EV) and on the complex cross-talk between the host and the nasal microbiome. To this aim, we evaluated the effects of PM10 and PM2.5 exposures on both the bacteria-derived-EV portion (bEV) and the host-derived EVs (hEV), as well as on bacterial nasal microbiome (bNM) features in 26 AR patients and 24 matched healthy subjects (HS). In addition, we assessed the role exerted by the bNM as a modifier of PM effects on the complex EV signaling network in the paradigmatic context of AR. We observed that PM exposure differently affected EV release and bNM composition in HS compared to AR, thus potentially contributing to the molecular mechanisms underlying AR. The obtained results represent the first step towards the understanding of the complex signaling network linking external stimuli, bNM composition, and the immune risponse.


Introduction
Exposure to particulate matter (PM) has been linked to the worsening of several respiratory allergies [1][2][3][4][5], including allergic rhinitis (AR), which affects more than 400 million people worldwide [6]. Although the underlying mechanisms linking PM exposure and allergy remain a matter of debate [7][8][9], it is becoming increasingly clear that PM can trigger local and systemic inflammatory responses involving innate and adaptive immunity, thus potentially enhancing the exacerbation of allergic symptoms [10,11]. In this context, plasmatic extracellular vesicles (EVs), which have a critical role in the cross-talk between cells [12][13][14], are a good candidate as mediators of the airway allergic process enhancement [15,16]. Moreover, EVs strongly contribute to the modulation of the inflammatory response after PM exposure, thus enforcing this hypothesis [17][18][19]. Interestingly, EVs represent a conserved communication mechanism between domains of life; indeed they can be released by human (i.e., hEVs) and bacterial cells (i.e., bEVs) into the bloodstream [13,14,20]. Mounting evidence suggests that the bEVs released by the microbiome communities inhabiting the host mucosae play a crucial role in microbiome-host interactions also influencing the host immune response [21,22] and potentially impacting the exacerbation of allergies [23,24]. It is worthy to note that the bacterial nasal microbiome (bNM), which plays a role at the interface between the environment and the host, and has a critical role in the mediation of this interaction, inlfuences both pathogenesis and exacerbation of airway allergies, including AR [25][26][27].
In this context, we hypothesize that the effects of PM exposure on the EV cross-talk and the nasal microbiome composition might affect AR pathogenesis. To this aim, in the present study, we evaluate the effects of PM 10 and PM 2.5 exposures on both bEV and hEV release and on bNM features in 26 AR patients and 24 matched healthy subjects (HS). In addition, we assess the role exerted by the bNM as a modifier of PM effects on the complex EV-signaling network to better elucidate the systemic effects induced by PM exposure in the paradigmatic context of AR.

Study Population
This study included 26 adult patients with AR and 24 healthy subjects matched for age, gender, and body mass index (BMI), recruited between February and June 2019 at the Allergology Unit, Fondazione IRCCS Ca' Granda-Ospedale Maggiore Policlinico, Milan, Italy. Each participant signed a written informed consent form, as approved by the Comitato Etico ("Milano Area 2" of the Fondazione IRCCS Ca' Granda Ospedale Maggiore Policlinico, 20,122 Milan, Italy (approval number 157_2019bis)), following the Helsinki Declaration principles. AR was diagnosed according to the Allergic Rhinitis and its Impact on Asthma (ARIA) guidelines [28]. Inclusion criteria included subjects with AR with a positive standard battery of the skin prick test for at least one inhalant allergen and/or eosinophil count. Exclusion criteria included diabetes, hypertension, autoimmune diseases, cancer, or other major chronic health conditions (e.g., chronic and/or allergic asthma), pregnancy, and history of illicit drug use. After signing a detailed informed consent form, all participants were asked to donate blood and nasal swab samples. In addition, each participant gave his/her written informed consent to participate in the study and filled out a standardized questionnaire about demographics and lifestyle information (e.g., smoking habits and occupation). Clinical history was also collected for each of the AR patients.

Particulate Matter Exposure Level Attribution
Air pollution data (PM 10 and PM 2.5 ) were retrieved from the Open Data Lombardy Region (https://www.dati.lombardia.it (accessed on 18 March 2020)) database, which contains daily estimates of municipal aggregate values calculated by the Regional Environmental Protection Agency (ARPA Lombardy). The assessment of the pollutant concentrations was based on the ARIA Regional Modelling (www.aria-net.it (accessed on 3 February 2020)), which is a chemical-physical model of air quality that simulates the dispersion and chemical reactions of atmospheric pollutants. It integrates the data measured from the monitoring stations of the ARPA Lombardy air quality network as well as meteorological data, emissions, concentrations at the beginning of the simulation period, and trends in adjacent areas, covering the whole Lombard territory with a grid of 1 × 1 km 2 cells and providing daily mean estimates available from the website at municipality resolution. Each subject was assigned the daily PM 10 and PM 2.5 concentration from (1) the municipality of residence in the 6 days preceding nasal sampling (i.e., day -1 to day -6) and from (2) the Municipality of Milan (i.e., where the biological sampling was performed) for the day of recruitment (day -0).

Sample Collection and Processing
Nasal swabs were collected and stored at −80 • C following standard porcedures. DNA was extracted using QIAamp ® UCP Pathogen Mini (Qiagen, Hilden, Germany) following the manufacturer's guidelines. The extracted DNA was stored at −20 • C and later shipped in cold packs to the sequencing service facility Personal Genomics Srl (Verona, Italy) to perform qualitative and quantitative checks, PCR amplification, and second-generation sequencing analysis. In particular, DNA quantification was performed by the Qubit dsDNA BR assay kit (Thermofisher Scientific, Waltham, MA, USA). Libraries were generated by following Illumina 16S Metagenomic Sequencing Library Preparation (Illumina, San Diego, CA, USA). Bacterial communities were investigated through amplicon sequencing analysis of the 16S rRNA gene hypervariable regions V3-V4, amplified with the primer pair Pro341F (5 -CCTACGGGNBGCASCAG-3 ) and Pro805R (5 -GACTACNVGGGTATCTAATCC-3 ). The obtained libraries were evaluated by Labchip DNA High Sensitivity (Perkin Elmer, Waltham, MA, USA), quantified by the Qubit dsDNA BR assay kit (Thermofisher Scientific), and then sequenced through the Illumina MiSeq platform (Illumina) using a paired-end library of 300 bp insert size.

Upstream Analyses and Operational Taxonomic Unit (OTU) Clustering
Raw reads' quality and statistics were checked using FastQC v0.11.2 and then trimmed at the 3 end using Trimmomatic v0.32 to improve the following read joining. The Fastqjoin.py tool was applied to join the raw reads and quality filtered using a minimum base quality of 20 (Phred-scale) over 5 bases' sliding windows, and were then analyzed using the default settings for QIIME 1.9.1 [29]. After chimeric reads removal, the resulting reads were clustered using 97% similarity, applying the open-reference OTU pipeline using USEARCH61 [30]. Taxonomic assignment was carried out with the RDP classifier [31] through the comparison of representative reads against the Greengenes v13.8 database using standard options. The PyNast method and default settings suggested in the QIIME pipeline were applied to align sequences [32,33]. The resulting OTU table was successively filtered, removing singleton and low abundance OTUs (i.e., <1% relative abundance) to perform downstream and statistical analysis. The final OTU table contained 1,941,731 reads (sequences length (mean ± standard deviation [SD]) = 408 bp ± 37 bp).

Downstream Analysis
Downstream analyses were carried out using QIIME v1.9.1 to analyze the abovedescribed OTU table. Taxonomic values within each sample and group were assigned to each OTU from the phylum to the genus level. OTUs, which failed genus attribution, were tagged as "Unassigned" followed by the specific family label. Before diversity analysis, all samples were rarefied to 10,000 sequences with a seed of 10 to avoid the influence of different sequencing depths, as this number of sequences was the minimum identified in the OTU table. Alpha-diversity richness, evenness, and genetic distance were calculated using observed species, Shannon, PD whole tree, Chao 1, and Simpson indices. Beta-diversity was examined by applying the weighted normalized UniFrac distance measure.

EVs Analysis
The isolation, purification, and characterization of EVs were performed by following MISEV 2018 guidelines [34].

Isolation and Purification of EVs
Blood samples were collected in ethylenediamine tetra-acetic acid (EDTA) tubes and processed within 2 h from the phlebotomy as described below. For isolation of plasma EVs, two aliquots of 3 mL of plasma for each subject were subsequently centrifuged at 1000, 2000, and 3000× g for 15 min at 4 • C. The obtained pellets were discarded to remove cell debris. EVs were then isolated from supernatants by ultracentrifugation at 110,000× g for 94 min at 4 • C in polypropylene ultracentrifuge tubes (Beckman Coulter; Brea, CA, USA) filled with phosphate-buffered saline (PBS) previously filtered through a 0.10 µm pore-size polyethersulfone filter (StericupRVP, Merck Millipore; Burlington, MA, USA). To carry out the Nanoparticle tracking and flow cytometry analyses, the EV-rich pellet was resuspended in 1 mL of triple-filtered PBS (pore size 0.1 µm). The methods described here are further detailed in the supplementary information (Supplementary Table S1).

Nanoparticle Tracking Analysis (NTA) of EVs
NTA analysis was carried out with the Malvern NanoSight NS300 system (Malvern Panalytical Ltd., Malvern, UK) used to visualize the EVs by laser light scattering. For each sample, five 30 s records were registered. NTA output was then analyzed with integrated NTA software (Malvern Panalytical Ltd.), providing high-resolution particle size distribution profiles and EV concentration measurements. NTA EV data were expressed as 10 6 for 1 mL of plasma.

Statistical Analysis
Standard descriptive statistics were performed on all variables. Categorical data are presented as frequencies and percentages. Continuous variables are expressed as the mean ± standard deviation (SD).
The differences between the two groups of AR and HS subjects were evaluated through Pearson's chi-square test or Fisher's exact test for categorical data, or t-test or the Mann-Whitney U-test for continuous variables, as appropriate.
To evaluate the association between PM 10 and PM 2.5 exposure on the different EVs subtypes, in AR and HS, we applied negative binomial regression for over-dispersed count observations. We tested the presence of over-dispersion based upon the Lagrange Multiplier (LM) test. Under the hypothesis that a short-term mechanism underlies variations in EV levels, we chose to investigate the exposure at the day of recruitment (defined as day -0). Models were adjusted for age, gender, smoking habit, mean temperature, mean humidity, case-control condition, and interactions between case-control and PM exposure.
Each model was tested for normality and linearity. All potential confounders were included in the multivariate model after verifying the presence of an association in a univariate model. Best model selection was based on the minimization of the Akaike infor-mation criterion and the maximization of the explained variance of the model. Estimated effects are reported as β and 95% confidence intervals (CI) associated with an increase of 10 units in each pollutant. To verify the feasibility of stratifying the analyses for allergy status (AR and HS), we tested the interaction between allergy status and exposure.
We analyzed the differences in AR and HS in terms of the distribution of vesicle marginal mean concentrations for each EV size. For each EV size: (1) we estimated marginal EVs' mean concentration and 95% CI at each size, with multivariate negative binomial regression models adjusted for age, gender, smoking habit, mean temperature, mean humidity, case-control condition, PM 10 (day -0) , and interactions between case-control and PM exposure; (2) we compared the EV mean differences between AR and HS; and (3) we calculated q-FDR values using the multiple comparison methods based on Benjamini-Hochberg False Discovery Rate (FDR) that takes into account the high number of comparisons, with a threshold of 0.10 to detect significance. Results were reported as a series graph for EV mean concentrations of each size and vertical bar charts to represent the size-specific p-values and q-values (i.e., FDR p-values) obtained comparing AR and HS. In the graphs, the x-axis corresponds to the size of EVs.
To assess differences between the OTU abundances and the alpha-diversity indices in the AR and HS groups, the t-test was applied.
Beta-diversity was examined by applying the weighted normalized UniFrac distance measure. To compare the tightness of the clustering distance within all the samples in the two groups comparing the state of disease, we performed a non-parametric permutational multivariate analysis of variance (NPMANOVA). This is used to compare groups of subjects and tests the null hypothesis that the centroids and dispersion of the groups, as defined by measure space, are equivalent for all groups. The test is based on the prior calculation of the distance between any two subjects and the rejection of the null hypothesis means that the centroid of the subject is different between the groups. Furthermore, to visualize and interpret the result of the applied distance measure, principal coordinate analysis (PCoA) was performed on the produced distance matrix and plotted using Emperor, and the adonis function in the R Vegan package was used to test the significance in dissimilarity matrices between AR and HS groups.
To evaluate the effects of PM 10 and PM 2.5 on α-diversity indices, for both the AR and HS groups, we applied multivariable linear regression models, adjusting for age, gender, smoking habit, mean temperature, mean humidity, allergy status, and interactions between allergy status and PM exposure. The association between PM and α-diversity indices was evaluated from the sixth (day -6) to the day of sampling (day -0), running several regression equations separately but reporting those on day -5, which were the only statistically significant associations we observed. Therefore, we hypothesized that day -5 of PM exposure might influence bNM composition and this, in turn, might affect the plasmatic EV release in response to acute PM exposure (i.e., day -0).
Factor analysis was applied to identify a small number of underlying independent variables and microbiomes factors that explain most of the variance that was observed in the much larger number of manifest variables of the microbiome composition at the genus level.
We excluded a priori 16 genera (i.e., "unassigned", and "other") because they did not provide any interpretable results and 21 genera with a very low relative abundance (mean ≤4% and median ≤0.3%). Next, we analyzed the correlation matrix of the logtransformed and standardized variables to investigate how genera correlated to each other and we took into account existing relationships among genera to avoid over-representation of single genera and the subsequent artificially higher correlations. Since three genera did not correlate (p-value < 0.05) with any other genera and correlation coefficients were less than |0.30|, they were not included in the factor analysis. We evaluated whether the correlation matrix of the relative abundances of the genera was factorable through both visual inspections of the matrix and statistical procedures, including Bartlett's test of sphericity, the Kaiser-Meyer-Olkin measure, and individual measures of sampling adequacy. Bartlett's test of sphericity tests the hypothesis that the correlation matrix is an identity matrix, which would indicate that genera variables are unrelated and therefore not suitable for structure detection. A significant p-value indicates that factor analysis may be useful with our data. The Kaiser-Meyer-Olkin measure of sampling adequacy indicates the proportion of variance in the genera variables that might be caused by underlying factors. The overall and individual measure of sampling adequacy ranges between 0 and 1. We considered an overall KMO ≤ 0.50 for the factor analysis and genera with a measure of sampling adequacy <0.30 as unacceptable. Thus, we excluded 17 genera and again verified the method assumption on the correlation matrix considering the remaining 24 genera. The new correlation matrix was factorable but five genera were excluded because of their low communality, i.e., they explained less than 15% of the variance each. In the last correlation matrix, all the assumptions were satisfied and brought the results obtained (Supplementary Table S2).
We carried out an explanatory factor analysis on the correlation matrix of 19 selected microbiome genera to derive a smaller set of uncorrelated underlying factors and to obtain the microbiome patterns (Supplementary Table S3).
We chose two factors to be included in the analysis considering the following criteria: factor eigenvalues >1, scree-plot construction, and factor interpretability [36]. We applied a varimax rotation to the factor-loading matrix to obtain a simpler loadings structure and make the factor matrix more interpretable. We calculated the factor scores using the mean principal components with iteration. Genera with an absolute rotated factor loading ≥0.63 on a given factor were used to name the factor and are indicated as "dominant genera" hereafter [37]. The factor scores, calculated for each subject, indicate the degree to which each subject's microbiome corresponded to the two identified patterns. To assess the reliability of microbiome patterns and internal consistency of genera that load more than |0.40| on any factor, we calculated Cronbach's coefficient alpha for each factor and coefficient alpha when the item was deleted.
We investigated if the two microbiome patterns identified by factor analysis could act as a modifier of the association between PM exposure and the different EVs subtypes. Thus, we investigated the role of microbiome patterns in the relationship between PM exposure and different EVs subtypes, applying multiple negative binomial regression models adjusted for age, gender, smoking habit, mean temperature, mean humidity, allergy status, PM 10 (day -0) , and the interactions between microbiome factors (Factor1 or Factor2, evaluated in separate models) and PM exposure. We reported the association between PM and different EVs subtypes adapted for the selected values of each factor (i.e., the 5th quantile, median, and the 95th quantile) and the other covariates.
All statistical analyses were performed by using SAS 9.4 statistical software (SAS Institute Inc., Cary, NC, USA) and R version 4.0.2.(R Core Team, Vienna, Austria)

Study Population and PM Exposure Data
The study population included 50 subjects (26 AR patients and 24 HS). No differences were found between AR (i.e., patients with allergic rhinitis) and HS (i.e., healthy matched controls) groups in any of the socio-demographic characteristics factors (Table 1). Patients' mean age was 41 years (±14.6 years) and the mean BMI was 23.8 kg/m 2 . Age, BMI, and gender values were matched in the two study groups. In total, 38% of enrolled subjects were smokers. Median PM 10 and PM 2.5 values in the week before the enrolment were 22.9 µg/m 3 and 14.9 µg/m 3 , respectively. A description of the PM 10 and PM 2.5 fractions for each considered day (from day -0 to day -6) is reported in Supplementary Figure S5.

EV Characterization in the AR and HS Groups
The median total number of EVs between AR patients and HS (i.e., total EVs within Table 2) was investigated by NTA analysis and AR patients were characterized by a higher concentration of plasmatic EVs than HS (EV median in AR: 1267 × 10 6 /mL vs. EV median in HS: 9688 × 10 6 /mL; p = 0.048).
Moreover, we characterized plasmatic EVs, focusing on both bEVs and hEVs subpopulations by flow cytometry analysis. Consistent with the NTA results, the median concentrations of all EVs subtypes were higher in AR patients compared to those of the HS group ( Table 2). The median concentrations of bEVs released from eosinophils (CD294 + EVs) and basophils (CD203C + EVs) were significantly increased in AR patients compared to the HS group.

Association between PM Exposure and Plasmatic EV Concentrations in AR and HS Group
We investigated the association between PM 10 and PM 2.5 exposure and both bEV and hEV plasmatic concentrations in the analyzed groups (Table 3). Considering the HS group, we observed a negative association between PM 10 (day -0) levels and the plasmatic concentration of LTA + EVs (β = −0.38, p = 0.013), LPS + EVs (β = −0.38, p = 0.023), and CD14 + EVs (β = −0.44, p = 0.013). No associations were found for AR patients. We further compared the AR patients and HS in terms of the distribution of marginal mean vesicle concentrations for each size (Figure 1). AR patients were characterized by a higher concentration of EVs ranging between 30 and 35 nm, 92-111 nm (p < 0.05; q < 0.10), with a peak at 99 nm. . EV concentrations were calculated as marginal means from a negative binomial model adjusted for age, gender, smoking habit, mean temperature, mean humidity, case-control condition, PM10 day -0, and interactions between case-control and PM exposure.

Compositional Overview of the Nasal Microbiota between AR and HS Participants
A sufficient yield for amplicon-based analysis was obtained from 23 AR patients and 23 HS samples. Considering the entire study population as a whole, the bNM was dominated by the Actinobacteria (abundance range AR: 2.5-83.5%; abundance range HS: 18.6-92.7%), Firmicutes (AR: 6.1-63.5%; HS: 6.4-72.2%), and Proteobacteria (AR: 0.6-89.6%; HS 0.7-38%) phyla. We compared the mean relative abundance of each phylum in the two groups and observed that the Proteobacteria phylum (mostly represented by Actinobacteria and the Firmicutes) was higher in the AR group compared to HS, although this evidence did not reach the significance level (Figure 2a). Interestingly, among the AR patients, six subjects were characterized by a very high presence of the Proteobacteria phylum. The remaining 17 AR patients seemed to exhibit a more heterogeneous composition with high relative abundances for the Actinobacteria and the Firmicutes phyla (Figure 2b). No differences were observed in terms of clinical, demographic, and lifestyle information between the six AR patients with a high prevalence of Proteobacteria phylum and the remaining 17 AR patients (data not shown).
We identified 81 genera represented in the study population. The most represented genera were Corynebacterium (AR: 2. Stratifying by the AR and HS group, the AR group showed a tendency of increased relative abundance for the Moraxella genus (belonging to Proteobacteria phylum), while Corynebacterium and Staphylococcus genera abundances were higher in the HS group (Figure 2c,d).

Compositional Overview of the Nasal Microbiota between AR and HS Participants
A sufficient yield for amplicon-based analysis was obtained from 23 AR patients and 23 HS samples. Considering the entire study population as a whole, the bNM was dominated by the Actinobacteria (abundance range AR: 2.5-83.5%; abundance range HS: 18.6-92.7%), Firmicutes (AR: 6.1-63.5%; HS: 6.4-72.2%), and Proteobacteria (AR: 0.6-89.6%; HS 0.7-38%) phyla. We compared the mean relative abundance of each phylum in the two groups and observed that the Proteobacteria phylum (mostly represented by Actinobacteria and the Firmicutes) was higher in the AR group compared to HS, although this evidence did not reach the significance level (Figure 2a). Interestingly, among the AR patients, six subjects were characterized by a very high presence of the Proteobacteria phylum. The remaining 17 AR patients seemed to exhibit a more heterogeneous composition with high relative abundances for the Actinobacteria and the Firmicutes phyla (Figure 2b). No differences were observed in terms of clinical, demographic, and lifestyle information between the six AR patients with a high prevalence of Proteobacteria phylum and the remaining 17 AR patients (data not shown).

Figure 2. bNM phylum and genus-profile composition. (a)
Average relative abundance of the bNM structure in the AR and HS groups; (b) relative abundance of the bNM of each enrolled subject; (c) average relative abundance in the AR and HS group; and (d) genus-profile structure of the bNM for each recruited subject.In addition, αdiversity analyses were performed to investigate compositional species richness and evenness within the study groups. The AR group was characterized by less diverse bNM compared to the HS group, as the bNM of the AR patients expressed a reduced mean richness and phylogenetic diversity compared to the HS group, as shown by the observed species (p = 0.009), Chao1 (p = 0.023), and the PD whole tree (p = 0.009) indices (Table 4). (c) average relative abundance in the AR and HS group; and (d) genus-profile structure of the bNM for each recruited subject. In addition, α-diversity analyses were performed to investigate compositional species richness and evenness within the study groups. The AR group was characterized by less diverse bNM compared to the HS group, as the bNM of the AR patients expressed a reduced mean richness and phylogenetic diversity compared to the HS group, as shown by the observed species (p = 0.009), Chao1 (p = 0.023), and the PD whole tree (p = 0.009) indices (Table 4) and Pseudomonas (AR: 0.002-11.6%; HS: 0.006-5.4%). Stratifying by the AR and HS group, the AR group showed a tendency of increased relative abundance for the Moraxella genus (belonging to Proteobacteria phylum), while Corynebacterium and Staphylococcus genera abundances were higher in the HS group (Figure 2c,d).
Distances between groups of subjects were additionally analyzed through the weighted normalized UniFrac distance metric. Intra-group distance distribution was significantly greater in the AR patients than in the HS (non-parametric test = 3.77, p = 0.01; Supplementary Figure S6a). To visualize this intra-group variability, a PCoA plot was generated (Supplementary Figure S6b). The first two axes, PC1 and PC2, were found to explain 66.9% of the bNM variability across the study population. The PCoA confirmed the greater dispersion of the AR subjects compared to the HS. However, no differences between the considered groups were identified (R 2 = 0.045 p = 0.10).

Exploratory Factor Analysis
An exploratory factor analysis was conducted to determine which latent variables (i.e., factors) modify the relationship between PM exposure and EV release. Factor analysis was applied to identify a small number of underlying independent variables, namely microbiomes factors, that explain most of the variance that is observed in the much larger number of manifest variables of the microbiome composition at the genus level.
The correlation matrix of the 19 selected genera (Figure 3) was suitable for factor analysis. Supplementary Table S3 reports the results on statistical procedures for checking matrix factorability. Bartlett's test of sphericity was significant (p < 0.001). The overall measure of the sampling adequacy was equal to 0.62, indicating that the sample size was sufficient, as compared to the number of genera under consideration. In addition, the individual measures of the sampling adequacy were satisfactory. Table 5 shows the factor-loading matrix for the two retained microbiome patterns, namely the corresponding communality estimates and the proportion of explained variance. The two retained factors explained about 66% of the total variance of the 19 genera, accounting for about 44% and 22%, respectively. The greater the loading of a given genus to one factor, the higher the contribution of that genera to the factor. The first factor, named Factor1, had the highest contribution from Streptococcus, Actinomyces, Rothia, and Fusobacterium. The second factor, named Factor2, was characterized by the greatest positive loadings on Anaerococcus, Acinetobacter, Peptoniphilus, and Staphylococcus. All the examined genera had at least onefactor loading greater than |0.39|, thus proving an important role of all the genera included in the analysis.  Estimated from a factor analysis performed on 19 genera. The magnitude of each loading indicates the importance of the corresponding genera to the factor. Loading ≥ 0.63 were shown in bold typeface; loading < |0.10| were suppressed.

PM Effects on bNM Composition
We evaluated the effects of PM 10 and PM 2.5 exposure on the α-diversity (from day -6 to day -0) for both the AR and HS groups. As shown in Table 6, in the AR subjects, positive associations were identified only between the fifth day preceding sampling PM 2.5 (Day-5) levels and observed species (PM 2.5 β = −34.2, p = 0.01), the PD whole tree (PM 2.5 β = 1.8, p = 0.036), and Shannon (PM 2.5 β = 1.2, p = 0.04) indices. On the contrary, negative associations were observed in HS among PM 2.5 and observed OTUs (PM 2.5 β = −22.5, p = 0.032), Shannon (PM 2.5 β = −0.97, p = 0.041), Chao1 (PM 2.5 β = −25.3, p = 0.097), and Simpson (PM 2.5 β = −0.2, p = 0.012) diversity indexes. No significant associations were observed with PM 10 exposure. Table 6. Association between PM exposure (day -5) and α-diversity indices. Adjusted for age, gender, smoking habit, mean temperature, mean humidity, case-control condition, and the interactions between casecontrol and PM exposure. β represents the change in the α-diversity index for a 10 µg/m 3 increase in PM concentration. Bold p-values denote a significant difference (p-value < 0.05).

Association between PM Exposure and Plasmatic EV Concentration in the Whole Population, Stratified According to Factors
We investigated the role of bNM in terms of microbiome patterns as a possible modifier effect of plasmatic EV release after PM exposure at day -0. The group with a bNM dominated by high levels of Factor1 (including Streptococcus, Actinomyces, Rothia, and Fusobacterium) showed negative associations with LTA+, LPS+, and CD14 + EVs when considering PM 10 (day -0) exposure. For PM 2.5 , we observed a negative association between high levels of Factor1 and CD14 + EVs ( Table 7). The group with a bNM dominated by high levels of Factor2 (including Anaerococcus, Acinetobacter, Peptoniphilus, and Staphylococcus) showed negative associations with LPS + when considering PM 10 (day -0) exposure No associations were observed for PM 2.5 (Table 8).

Discussion
The prevalence of respiratory diseases is rising dramatically worldwide, turning respiratory allergies into an "epidemic" phenomenon [37][38][39]. In particular, severe AR has been associated with significant impairments in quality of life, work performance, and sleep [39,40]. Increasing epidemiological studies conducted both on adults and children have been reporting an associative relationship between air pollution, including PM, and increased incidence of AR [7,9,[41][42][43][44][45][46]. Indeed, allergens and environmental agents trigger reactive nasal inflammatory conditions such as AR via the involvement of host intrinsic factors, including the innate and adaptive immune system, as well as the bNM membership [47,48]. In order to highlight the molecular mechanisms linking air pollution exposure and AR, we investigated both plasmatic EVs and bNM as potential actors in this complex scenario.
We assessed plasmatic EV release in AR compared to HS and the higher plasmatic EV concentrations found for the AR group could be linked to pro-inflammatory activation, as it was previously reported [49][50][51]. Regarding the hEVs fraction, the increased concentrations of both mast cells and neutrophils-derived EVs may denote the enhanced IgE-based immune response, a specific signature of allergic conditions [52,53]. We observed a tendency of increased concentration also for monocyte and platelet-derived EVs, which release was reported to be increased during proinflammatory stimuli, also mediated by environmental exposure [54][55][56]. When considering the PM effect, the AR and HS groups showed an opposite tendency between day -0 PM 10 exposure and EVs released by both gram-positive and gram-negative bacteria (i.e., LTA+ and LPS + EVs, respectively), although only the negative associations in HS reached the statistical significance level. This evidence suggests that PM exposure influences bNM cross-talk with the host, affecting bEV release [57]. In particular, in-vitro studies emphasized that bEVs can trigger inflammation through the activation of the toll-like receptors 2 and 4 [58,59].
The complex interaction between the microbiome and the host is at the basis of human health. Microbiome dysbiosis can be caused by exposure to several environmental pollutants, including PM, and has been associated with numerous diseases. Nonetheless, the role of the microbiome in the multifaceted host-environment interplay is far from being fully understood [21]. Both HS and AR bNMs were dominated by the Actinobacteria, Proteobacteria, and Firmicutes phyla according to previous studies [21]. Interestingly, the high Proteobacteria-relative abundance, which characterized the bNM of the AR patients, has been linked to several respiratory diseases, such as chronic obstructive pulmonary disease (COPD), tuberculosis, and asthma [60][61][62].
The importance of a diverse microbiome has been widely documented regarding maintaining mucosa integrity and an effective immune system [21,63]. Thus, the less diverse bacterial community observed in the AR group compared to the HS group may be linked to the disease condition. In particular, reduced diversity in the bNM was identified also in other respiratory diseases such as chronic rhinosinusitis [23,64], asthma [65], and granulomatosis with polyangiitis [66]. This evidence is further supported by the opposite behavior observed between PM 2.5 and α-diversity indices in HS and AR, possibly leading to an imbalanced or dysfunctional microbial community with an impact on disease severity and inflammatory response [21,57,[67][68][69][70][71].
Since bNM is a complex and dynamic community [72], we applied factor analysis to obtain a more comprehensive picture of the most representative bNM composition in the study group, thus considering a smaller number of independent factors able to predict microbiome composition at the genus level. The factorial analysis allowed us to further investigate whether bNM is able to modify the effects induced by PM exposure on both bEVs and hEVs concentrations. To this aim, we considered the two factors (i.e., Factor1 and Factor2) identified, tested the interaction between the two categories and the PM exposure in the multivariable regression model, and observed different associations depending on the factor. This evidence supports our hypothesis that bNM composition is influenced by PM and, in turn, this effect modifies the acute plasmatic EV release after PM exposure. In accordance with this, we recently reported that in healthy subjects, bNM membership was linked to the pattern of plasmatic EV release when PM was considered [73].
We need to acknowledge that the study has some limitations. The present observational study was conducted by taking into account a size-limited population, which, consequently, may not be sufficient to uncover all the variations induced by the PM effect, especially when either EV concentration or the bNM feature were considered. Moreover, we estimated PM exposures at the participants' residential addresses, without considering their daily activities and locations, nor the indoor PM levels, which may have resulted in misclassifications of the exposure attributions. However, the consistency of our findings with the literature makes us confident that the obtained results were not found by chance and help to understand the complex signaling network, thus linking external stimuli, bNM composition, and the immune system.

Conclusions
In the present study, we observed an association between PM exposure in the context of both bNM composition and plasmatic EV release in patients with AR. Moreover, we observed that PM exposure differently influences plasmatic EV release depending on the bNM composition. To the best of our knowledge, this is the first study specifically investigating this association. Longitudinal studies will help to clarify the link between PM exposure and both bNM modulation and plasmatic EV release in those subjects with the AR condition. Moreover, further functional studies are needed to better characterize the different responses observed in the AR and HS groups after PM exposure, as well as to investigate the specific cell subpopulations involved in this mechanism and the consequent gene expression variations.