Characterization of long-range transported bioaerosols in the Central Mediterranean.

Airborne bacteria were characterized over a 2-y period via high-throughput massive sequencing of 16S rRNA gene in aerosol samples collected at a background mountain European Monitoring and Evaluation Programme (EMEP) Network site (Monte Martano, Italy) located in the Central Mediterranean area. The air mass origin of nineteen samples was identified by air mass modelling and a detailed chemical analysis was performed. Four main origins (Saharan, North-western, North-eastern, and Regional) were identified, and distinct microbial communities were associated with these air masses. Samples featured a great bacterial diversity with Protobacteria being the most abundant phylum, and Sphingomonas followed by Acidovorax, Acinetobacter and Stenotrophomonas the most abundant genera of the dataset. Bacterial genera including potential human and animal pathogens were more abundant in European and in Regional samples compared to Saharan samples; this stressed the relevance of anthropic impact on bacterial populations transported by air masses that cross densely populated areas. The principal aerosol chemical characteristics and the airborne bacterial communities were correlated by cluster analysis, similarity tests and non-metric multidimensional scaling analysis, explaining most of the variability observed. However, the strong correlation between bacterial community structure and air mass origin hampered the possibility to disentangle the effects of variations in bacterial populations and in dust provenance on variations in chemical variables.

Introduction

The presence and diffusion of bioaerosols (bacteria, viruses, fungi, and other dead or living organisms including biological debris) in the Earth atmosphere impact ecosystems, climate, and human health (Burrows et al., 2009a; Burrows et al., 2009b; Fröhlich-Nowoisky et al., 2016; Pöschl and Shiraiwa, 2015). The biosphere directly emits bioaerosols into the atmosphere, which subsequently enables their dispersion and transport even at long distances (Després et al., 2012; Womack et al., 2010). In the course of atmospheric transport, bioaerosols may undergo further chemical and physical transformation, stress, and biological aging upon interaction with UV radiation, photo-oxidants, and various air pollutants like acids, nitrogen oxides, ozone, and aromatic compounds. All these processes can limit or even suppress the vitality of the living fraction of bioaerosol and therefore affect their capacity to diffuse and to colonize new ecosystems (Womack et al., 2010). Due to the above challenges, the present knowledge on the ability of viruses and bacteria to spread in the air and diffuse infections and more in general diseases is still immature and demands a wide spectrum of investigation (Middleton, 2017; Morawska and Cao, 2020; Polymenakou, 2012).

Most of the previous studies in the Mediterranean area have been limited to advections of air masses from the Sahara Desert only. The occurrence and impact of this type of air mass, very rich in desert dust, are frequent and well documented (Escudero et al., 2006; Formenti et al., 2011; Goudie and Middleton, 2001; Pey et al., 2013). Saharan dust particles can be transported over long distances towards Europe (Amato et al., 2016; Barnaba et al., 2017; Cusack et al., 2012) and America (Garrison et al., 2014; Prospero et al., 2005). However, although Saharan dust is recognized as one of the most relevant sources of atmospheric aerosol and bioaerosol on a global scale, the Central Mediterranean is a geographic area profoundly affected by the circulation of air masses of different origin and distinguished by nature, type, quality, and extent of contributions (Cusack et al., 2012; Kallos et al., 2014, Kallos et al., 2007; Petroselli et al., 2018a, Petroselli et al., 2018b). Due to the very different characteristics of the source areas, these air masses are expected to carry different bacterial populations and specific chemical markers and pollutants. Moreover only a few studies have used molecular-based approaches to investigate the relationships of different air masses with the bacterial communities in the Mediterranean area. In such studies, bioaerosol characterization was conducted by a low-throughput approach (cloning and sequencing of 16S rRNA gene), while High-Throughput Sequencing (HTS) approaches were used in an even smaller number of cases. Most of the previous studies on aerosol-associated microbial communities in the Mediterranean area have been focused on intense Saharan intrusions sampled in the proximity of the dust sources (Gat et al., 2017; Katra et al., 2015; Mazar et al., 2016; Polymenakou et al., 2008), or after a long-range transport over the Mediterranean basin (Federici et al., 2018; Rosselli et al., 2015; Sanchez De La Campa et al., 2013). Much less is known about the specific characteristic of the bacterial communities transported by air masses from continental Europe.

In this frame, the present study aims at defining the patterns of the bacterial communities of atmospheric aerosol from distinct geographic regions reaching the Mediterranean. The samples were collected during different long-range transport events towards a background monitoring site located on the top of a mountain range (1100 m asl), in Central Italy. The low background aerosol concentrations make the site well suited to characterize long-range air mass transports. Samples were selected based on their provenance, unambiguously individuated by a modelling of the air masses. An HTS approach, combined with a thorough chemical analysis, allowed us to build a significant dataset including bacterial diversity and chemical composition. The hypotheses to be tested in this work are two-fold: (i) bacterial community structure associated with long-range transported aerosol in the Central Mediterranean area is significantly different based on the air mass provenance; (ii) there is a correlation between the main aerosol chemical characteristics and the airborne bacterial communities. To test these hypotheses, we investigated the chemical and microbial datasets by cluster analysis, similarity tests, and non-metric multidimensional scaling analysis.

Material and methods

Aerosol sampling

All the aerosol samples analyzed in this work were collected at the EMEP regional background site of Monte Martano (MM) in Central Italy (42°48′19″N, 12°33′55″E). MM has been established in a relatively undisturbed location, near a television antenna, on the ridge of a small mountain chain (1100 m asl), above the timberline and facing a completely free horizon (Moroni et al., 2015). The site is equipped with aerosol, gaseous pollutants, and meteorological monitoring instrumentations (Moroni et al., 2015). Due to its elevation, the low background concentrations and the 360° free horizon, the site is particularly suited for the assessment of long-range transport events of atmospheric aerosol (Federici et al., 2018; Petroselli et al., 2018a, Petroselli et al., 2018b). The importance of the site for the monitoring of Saharan dust advections was recognized in 2013 when it joined the WMO SDS-WAS1 network. PM10 and PM2.5 are measured daily at MM with standard low volume impactors and the dust load is estimated following the procedure by (Escudero et al., 2007). Moreover, PM10 samples are also collected weakly by means of a high-volume sampler (HVS, TISCH, TE6001, 1140 L min−1) on quartz fiber filters (Whatman QMA, 20 × 25 cm) previously sterilized, together with the filter holder, under UV lamp for 50 min on both sides. In the period 2014–2015, for the present campaign, a set of additional HVS samplings were planned on a daily basis thanks to the forecast system based on DREAM8b model simulations run by the SDS-WAS Barcelona Supercomputing Center (BSC) (see Supporting Information, Fig. SM1) and by complementary forecast back trajectory (BT) calculations. Typically, we sampled 1 filter for 24 h for each advection event. In one case (30 Nov 2014–1 Dec 2014) we were able to take two successive 24 h samples representative of the same long-lasting event (Federici et al., 2018). The filters were collected immediately at the end of the sampling and kept stored in the freezer at −20 °C until the processing. A significant portion of the filter (20 × 10 cm) was dedicated to the microbiological analysis while two smaller slices (2 × 2 cm) were employed for ion chromatography and for elemental analyses, respectively. Laboratory blank and field blank filters have been characterized at the beginning and the end of the campaign for checking the sampling protocol for sterility and quantification of limit-of-detections for chemical analyses.

Successively, the provenance of every sample was verified using BTs based on reanalysis of meteorological fields. BTs were calculated with the HYSPLITv4 model (Stein et al., 2015) computed hourly for the period of interests, considering three endpoints located at 50, 500 and 1000 m a.g.l. and exploiting meteorological fields from GDAS (Global Data Assimilation Service) with 1 × 1 degree spatial resolution. The vertical structure in the atmosphere of the bacterial community is also of relevance (González-Toril et al., 2020) and the DREAM8b model allowed also to check the vertical distribution of dust over the MM site (see Fig. SM2).

Nineteen of these HVS samples (9 Saharian, 4 Regional, 4 Northwestern, 2 Northeastern), whose provenances were clearly identified, are considered in the present work.

Chemical analyses

The sampled filters underwent a thorough chemical characterization that included the investigation of both the inorganic and organic fractions of particulate matter.

Major ion composition was determined by ion chromatography (DIONEX 2100) after 30 min ultra-sonication in ultrapure water (18 MΩ). The quantified analytes were: Li+, Na+, NH4 +, K+, Mg2+, Ca2+, F−, HCOO−, MSA, Cl−, NO2 −, SO4 2−, oxalate, Br−, NO3 −, PO4 3−.

Concerning the elemental composition, the samples underwent a microwave-assisted acid-digestion process (CEM MARS-5) being treated with nitric acid (HNO3 65%, Millipore Suprapur) and hydrogen peroxide (H2O2 30–32%, Carlo Erba Reagents) in 4:1 proportion. The samples were then analyzed by ICP-AES (Horiba ULTIMA 2000) equipped with an ultrasonic nebulizer (CETAC 5000). The determined elements were: Zn, Cr, Ni, Cu, V, Mn, Co, Ca, Ti, Fe.

Thermal-Optical-Transmittance (TOT, Sunset Laboratory Inc.TM) was used to determine elemental and organic carbon concentrations (EC/OC), following the NIOSH protocol. Finally, gas chromatographic-mass spectrometric analysis (GC Chrompack 3800 coupled with ITD-MSn Saturn 2000-Varian) was performed for the characterization of PAH and n-alkane fractions (Cartechini et al., 2015; Federici et al., 2018).

DNA extraction and sequencing

Total DNA extraction

Filters were handled in aseptic conditions, i.e. under a biological laminar flow hood. A portion of the filter (6 × 6 cm) was suspended in 40 mL of sterile water, shaken for 1 h at maximum speed, centrifuged for 30′ at 10000 ×g and then at 11500 ×g for 15′ at 4 °C to recover bacteria (Radosevich et al., 2002). Supernatant was discarded and DNA was extracted from the pellet using the POWER SOIL DNA kit (MOBio) following the manufacturer's protocol.

16S rRNA gene fragment amplification and Illumina sequencing

Bacterial communities were studied with an HTS approach by Illumina MiSeq sequencing. The V5-V6 hypervariable regions of the 16S rRNA gene were amplified in 2 × 50 μL volume reactions with a Hot start Taq polymerase (Solis-Biodyne). The reaction included 1 μM each of primers 783F and 1027R (Huber et al., 2007; Wang and Qian, 2009). At the 5′ end of each primer one 6-bp barcode was also included to allow sample pooling and subsequent sequence sorting. The cycling conditions were: initial denaturation at 94 °C for 5 min; 29 cycles of 94 °C for 50 s, 47 °C for 30 s, and 72 °C for 30 s and final extension at 72 °C for 5 min. The amplified products were purified with the Wizard SV PCR purification kit (Promega Corporation) and DNA quantity and purity were spectrophotometrically evaluated by Qubit (Invitrogen). Groups of 9 purified amplicons bearing different barcode pairs were pooled to obtain a single library. Further library preparation with the addition of standard Nextera indexes (Illumina, Inc., San Diego, CA, USA) and sequencing were carried out at Parco Tecnologico Padano (Lodi, Italy). Multiplexed sequencing of all the pooled samples was performed in a single Illumina MiSeq run, using a paired-end 2 × 250 base-pair protocol and the 4.0 sequencing chemistry.

Sequence analysis

Each sequence was assigned to its original sample according to its index oligos and barcodes. After sorting the sequences, the reverse read of each paired-end sequence was reverse complemented and merged with the corresponding forward read. A quality cut-off was applied in order to remove the sequences that did not contain the barcode, those with an average base quality value (Q) lower than 30 and those that did not provide a perfect match in the overlapping part between the two paired ends. The barcode was removed and sequences were sorted into Operational Taxonomic Units (OTUs) using the UPARSE-OTU algorithm (Edgar, 2013). The minimum identity between each OTU member sequence and the representative sequence (i.e. the sequence that showed the minimum distance to all other sequences in the OTU) was set to 97%. The taxonomic classification of each OTU was carried on with the stand-alone version of RDP Bayesian Classifier (Wang et al., 2007), using a 50% confidence level (Claesson et al., 2010). Chloroplast sequences were not excluded by further analyses because their abundance can provide information on PM origin.

Three independent extractions, amplifications and sequencing on each sample were performed in order to test the robustness of the proposed experimental approach and the three replicates featured nearly identical OTU distribution profiles (data not shown).

Statistical analyses

Cluster analysis using the Bray-Curtis similarity index was applied to the bacterial communities belonging to the different aerosol samples. Similarity test (ANOSIM) was performed to detect differences in the bacterial community structure followed by the determination of discriminating genera by means of SIMPER routine. This analysis indicates the average contribution of each genus to the similarity and dissimilarity between groups of samples. Non-metric Multidimensional Scaling (NMDS) analysis was performed using the Bray-Curtis dissimilarity matrix and the first NMDS dimension was then plotted with chemical data in order to gain information from the correlation between abiotic and biotic components of dust samples. Additionally, the chemical peculiarities of the samples based on similarities highlighted by the NMDS were interpreted by using principal component analysis (PCA). Statistical analyses and graphical representations were carried out using the R statistical environment (Version 4.0.1 - R Core Team, 2020) and ggplot2 package (Wickham, 2016).

Results and discussion

General characteristics of air masses and PM samples

Nine Saharan dust advections and ten long-range transports from other geographical origin have been considered in this work. The air mass origins were identified on the basis of back-trajectory (BT) analysis. The BTs for the identified provenance groups are summarized in Fig. 1 , for the 500 m endpoint. The other endpoints (50 and 1000 m above the ground) provided similar results and have been included in supplementary material (Fig. SM3). Saharan dust advection samples have been marked with the code SH. As for the other provenances, three main macro-areas have been identified, namely regional (RG), North-western (NW) and North-eastern (NE). RG air masses have been defined as those remaining over the terrestrial and marine sectors of central Italy for at least 48 h before sampling. Table 1 summarizes all the PM samples collected during 2014 and 2015.

Fig. 1

Air masses arriving at the MM site at 500 m above ground level (see text). Back-trajectories were calculated hourly (grey lines) for the sampling time corresponding to each sample and grouped in four cases: Saharan (SH, lefthand upper panel), Northeastern (NE, righthand upper panel), Northwestern (NW, lefthand lower panel), and regional (RG, righthand lower panel). Red lines are plots of the average of each group of BTs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

Sample characteristics in terms of provenance and aerosol mass concentration in the PM10 and coarse (PM10-PM2.5) fractions. Provenance classification is based on BTs analysis (see Fig. S1 in the Supporting Information).

Sample code	Provenance	PM10[μg m−3]	PMcoarse[μg m−3]
SH_20140220	Northern Africa – Algeria	19.3	11.7
SH_20140404	Northern Africa – Tunisia – Libya	27.5	13.4
SH_20140522	Northern Africa – Algeria	18.7	7.2
SH_20140624	Northern Africa – Tunisia – Mediterranean	12.7	3.9
SH_20141015	Northern Africa – Algeria – Tunisia – Libya	28.4	15.7
SH_20141107	Northern Africa – Libya	8.4	0.9
SH_20141130	Northern Africa – Tunisia – Libya	83.9	51.6
SH_20141201	Northern Africa – Tunisia – Algeria	86.9	55.7
SH_20150506	Northern Africa – Morocco – Algeria	30.1	17.5
NE_20140310	North East – Eastern Europe	17.5	6.8
RG_20140424	Regional – NE	7.1	0.9
NW_20141022	North West – France	6.1	1.7
NE_20141029	North East – Eastern Europe	11.5	2.1
NW_20141212	North-North West – France	6.3	1.1
RG_20150315	West-Tyrrhenian – Regional	19.6	4.8
RG_20150526	Regional – NE	5.2	0.3
RG_20150601	Regional – NW	15.3	5.0
RG_20150607	Regional – NE	13.3	2.3
NW_20150615	West – Iberian Peninsula	9.5	4.2

As a general trend, the Saharan dust samples are characterized by higher aerosol mass concentrations with respect to the non-Saharan advections (see Table 2 ), i.e. an average +68.4% for PM10 and +85.3% for PMcoarse, defined as PM10-PM2.5, and lower PM2.5/PM10 ratio, reading 0.52 ± 0.18 for SH and 0.76 ± 0.12 for non-SH samples. The increase in the concentration of the coarse fraction is typical of natural crustal aerosol sources such as desert dust (Formenti et al., 2011). Moreover, Ca and Fe, typical crustal markers resulted higher in Saharan dust on average (Table 2). The insoluble fraction of Ca, defined as Catot-Ca2+, was close to 60% for Saharan dust, slightly lower for RG air masses and much lower for the NW. This is consistent with both the source area mineralogy and the different atmospheric processes during the long-range transport (Avila et al., 2007). Biomass burning markers such as ammonium and organic carbon (OC) were higher in non-Saharan samples, and particularly enriched in NE samples, possibly due to the frequent wildfires recorded in Eastern Europe regions. The latter have been found to exert a distinct impact on the Monte Martano site, as previously reported in Petroselli et al., 2018a, Petroselli et al., 2018b. The average OC and EC values are in agreement with those reported in (Sandrini et al., 2014) at MM for the year 2009. Total PAHs were on average the highest for RG followed by SH and NE and NW air masses. Benzo(a)Pyrene, the reference PAH for health effects, has the same order in abundance.

Table 2

Aerosol mass concentration and chemical variables average values. Comparison between the investigated Saharan dust advections and non-Saharan (NE, NW, RG and total non-SH average) samples. Errors are given as standard deviations. Mean values for the 2014/2015 are also reported. All data in μg m−3 except when indicated.

	SH	NE	NW	RG	2014/2015
PM10	35.1 ± 29.4	14.5 ± 4.2	7.3 ± 1.9	12.1 ± 5.9	10.6
PM2.5	15.4 ± 9.6	10.1 ± 0.9	5.0 ± 0.5	9.4 ± 4.0	7.7
PMcoarse	19.7 ± 20.0	4.5 ± 3.3	2.3 ± 1.6	2.7 ± 2.2	3.9
PM2.5/PM10	0.52 ± 0.2	0.71 ± 0.2	0.70 ± 0.1	0.81 ± 0.1	0.73
Catot	3.0 ± 3.0	0.14 ± 0.2	0.38 ± 0.3	0.66 ± 1.0
Ca2+	1.2 ± 0.3	0.3 ± 0.2	0.3 ± 0.2	0.3 ± 0.2	0.7 ± 1.0
Fe	1.5 ± 2.1	0.18 ± 0.1	0.19 ± 0.17	0.13 ± 0.05	0.7 ± 1.3
Ti	0.06 ± 0.07	0.02 ± 0.05	0.02 ± 0.05	0.1 ± 0.08	0.05 ± 0.01
NH4+/NO3−	0.10 ± 0.09	0.46 ± 0.2	0.39 ± 0.2	0.46 ± 0.8	0.29 ± 0.2
NO3−/SO42−	1.8 ± 0.6	0.64 ± 0.3	2.9 ± 0.2	0.68 ± 0.8	0.9 ± 0.5
K+	0.15 ± 0.2	0.19 ± 0.01	0.05 ± 0.04	0.07 ± 0.1	0.26 ± 0.1
SO42−	1.4 ± 1.4	4.2 ± 0.5	1.2 ± 1.5	2.3 ± 0.9	1.7 ± 1.4
OC	3.3 ± 1.5	3.9 ± 1.7	1.9 ± 0.5	3.2 ± 1.7	2.8 ± 0.9
EC	0.3 ± 0.2	0.3 ± 0.1	0.12 ± 0.03	0.2 ± 0.1	0.20 ± 0.1
ΣPAH (pg m−3)	280 ± 150	200 ± 130	53 ± 14	440 ± 700	230 ± 220
BaP (pg m−3)	11 ± 8	7.9 ± 5	1.3 ± 1.3	16 ± 28

Characterization of bacterial communities

Sequencing of 16S rRNA gene fragments led to the recovery of 1286659 high-quality sequences, which clustered, across all samples, into a total of 10,513 operational taxonomic units (OTUs) calculated at 97% of sequence similarity. The average number of OTUs per sample was 2239. Although a considerable fraction of the total biodiversity (18.6% on average) could not be classified at genus level, a total of 879 different genera were identified across all samples. Among them, a total number of 116 genera were found whose relative abundance was higher than 0.5% in at least one sample. These genera were considered abundant (abundant genera hereafter) and further analyzed. 32 bacterial genera manifested a relative abundance higher than 0.5% on average in all samples; they are shown in Fig. 2 . Overall, the most abundant genera were Sphingomonas (8.47%), followed by Acidovorax (3.89%), Acinetobacter (3.33%), Stenotrophomonas (3%), Hymenobacter (2.68%), Methylobacterium (2.57%), Propionibacterium (2.1%) and Massilia (1.85%). Most of these genera have already been described as members of airborne bacterial assemblages, both from urban and suburban environments and from dust storm intrusions (). Particularly, the beta-proteobacterium Massilia (Oxalobacteriaceae) was found in many extreme environments impacted by desert dust (Chuvochina et al., 2011); however, it was also consistently retrieved in urban and suburban Chinese areas, together with Sphingomonas (Gao et al., 2017; Wei et al., 2015). Methylobacterium is known to be resistant to desiccation and to γ radiation together with Arthrobacter, also abundant in our samples (Favet et al., 2013). Microvirga was already found in desert-coming air-masses and some species of this genus can reduce nitrogen gas to ammonia (Favet et al., 2013; González-Toril et al., 2020). It should be also noted that some of the 32 most abundant genera are human- and animal-associated bacteria and include known pathogens, such as Haemophilus, Staphylococcus, Streptococcus and Propionibacterium (Brock et al., 2012). Moreover, we retrieved some bacterial genera that, despite being ubiquitous in the environment, also contain many opportunistic pathogens and a few pathogens, such as the Pseudomonadales Acinetobacter and Pseudomonas, or Clostridium sensu stricto and Clostridium XI (Brown, 2014). However, analyses based on 16S rRNA sequences do not allow to distinguish pathogenic from non-pathogenic species or strains.

Fig. 2

Cluster dendrogram and associated barplot of the relative abundance of the most abundant genera in each air-mass sample. Cluster analysis was performed on the 116 genera whose abundance was higher than 0.5% in at least one sample (abundant genera). Barplots represent only the 32 bacterial genera that showed a relative abundance higher than 0.5%. Chloroplasts were also included in the analysis.

Dependence of bacterial diversity on air mass origin

The results of the cluster analysis on the database containing only the abundant genera, summarized by the dendrogram in the right panel of Fig. 2, were used for a data-driven visualization of the samples that are reported in the barplot following the dendrogram order. The average 1-D distances reported in the dendrogram revealed a high β-diversity among the bacterial communities of the SH samples, which however generally clustered together. At genus level, the structure of bacterial communities clearly showed differences due to the sample provenance rather than to other factors such as seasonality. Interestingly, among the PM samples, the non-Saharan samples collected during regional movements of air masses (RG) and, to a lesser extent, during long-range intrusions (NE and NW), showed a high richness of genera with low abundance, indicating a highly diverse and even community. Conversely, PM during Saharan intrusions showed a lower richness of genera, indicating that these microbial communities were dominated by fewer typical phylotypes. This is in contrast with previous studies, which generally reported higher diversity during dust intrusion events compared to non-dust events (González-Toril et al., 2020; Griffin, 2007; Mazar et al., 2016; Polymenakou et al., 2008; Sanchez De La Campa et al., 2013). It may be hypothesized that, in the case of central Italy, both regional air masses and long-range intrusions from NE and NW mainly cross more heterogenous areas than Saharan intrusions, thus collecting a wider variety of microorganisms. Particularly, when BTs indicated that regional air masses were prevalent, it is also possible that local sources played a major role in shaping bacterial communities. In fact, a wide variety of potential local sources, such as soil surface, leaf surface, water bodies and even animal faeces, have been reported to influence airborne bacterial and fungal community structure (Bowers et al., 2011; Mu et al., 2020; Qi et al., 2020; Sadyś et al., 2014).

To gain further insights about similarities and differences among all the aerosol samples and hypothesize possible effects of the air masses with different origins, a non-metric multidimensional scaling (NMDS) analysis was performed using the computation of Bray-Curtis distances between bacterial communities. This statistical method has been applied to the database containing only the abundant genera, and results are shown in Fig. 3 . The results of NMDS analysis showed a good clustering (stress value <0.15) of the samples according to their provenance group. In particular, the NMDS1 dimension (x-axis, Fig. 3) seems to separate well the different clusters, and particularly the Saharan dust samples from the others, while NMDS2 describes the variability within each group.

Fig. 3

Representation of NMDS analysis of air-mass samples, based on the 116 abundant genera (see text). Chloroplasts were also included in the analysis.

NMDS is a helpful exploratory analysis but does not allow explaining the similarities or dissimilarities among samples, and additional information from other analysis is needed. For example, samples associated with RG and NW air masses show a partial overlap in the NMDS analysis, which is understandable based on the phenomenology of back trajectories (see Fig. 1). In fact, RG air masses at MM tend to the terrestrial and marine western sectors of peninsular Italy while samples of NE origins, although limited in number in present study, form a distinct cluster in NMDS and show a different pattern in the back-trajectory plot. On the other hand the similarity between NW_20141022, SH_20150506, and RG_20150607 samples as for the NMDS1 can be explained considering the source characteristics. In order to support the NMDS results we run a principal component analysis (PCA, results in Fig. SM4, Supplementary Material) based on chemical variables. The result indicates that similarity of the three samples is also driven by common chemical characteristics. In particular, the Na concentration, combined with the Benzo(a)Pyrene to Benzo(e)Pyrene ratio, suggests a similar residence time in the atmosphere and collective marine influence (Simo et al., 1991), configuring the SH_20150506 sample as an atypical Saharan dust intrusion.

Furthermore, we performed a SIMPER analysis to highlight the genera that typify each community. The results of pairwise comparison between the dust samples individuated by the above analysis is shown in Fig. 4 . In the plot genera are reported along descending triangles based on their average contribution to the Bray-Curtis dissimilarity. SH samples are represented at the center of the plot, for convenience, and share common genera with all the other groups. The genera that mostly contributed to characterize Saharan intrusions were Sphingomonas, Acidovorax and Stenotrophomonas. Although Sphingomonas represents a quite ubiquitous genus, it is often found in desert-originated air samples (Griffin et al., 2006). Family Comamonadaceae, to which the genus Acidovorax belongs, has been reported both in Saharan dust layers deposited in Alpine snow (Meola et al., 2015) and in Asian long-range transport plumes (Smith et al., 2013). By contrast, Stenotrophomonas was more frequently described as a major member of bioaerosols from anthropic activities, such as wastewater treatment (Han et al., 2019) and livestock and poultry production (Chien et al., 2011). Moreover, SH samples were characterized by higher abundances of genera Bacillus, Adhaeribacter and Hymenobacter than NW and RG samples. Bacillus has been often described as one of the most abundant genera in Saharan dust samples (Griffin, 2007; Sanchez De La Campa et al., 2013). Adhaeribacter and Hymenobacter can be particularly abundant in arid soils and desert sands (Favet et al., 2013).

Fig. 4

Representation of SIMPER analysis and pairwise comparisons of aerosol samples. Genera are reported along the descending triangle on the basis of their average contributions to the average overall Bray-Curtis dissimilarity. Chloroplasts were also included in the analysis.

In contrast, chloroplasts were significantly more abundant in non-Saharan than Saharan samples. This is not surprising, since areas crossed by air masses from RG, NW and NE provenances are more vegetated than northern Africa. Interestingly, most of the other genera that were more abundant in the European or in Regional samples are human and/or animal associated bacteria, such as Lactobacillus, or they also include pathogens, such as Streptococcus, Propionibacterium and Haemophilus (Brock et al., 2012; Brown, 2014). An interesting case is represented by Acinetobacter: since it is ubiquitous in soil and water and easily transported by air, it has recently proposed as an indicator of airborne transport from the Sahara (Barberán et al., 2015). However, despite its presence in our Saharan samples, it was significantly more abundant in the NE group. Overall, we could therefore hypothesize that bacterial communities transported by non-Saharan air masses were significantly affected by the crossing of densely populated areas. Therefore, although it has been demonstrated that a number of diseases are linked to desert aerosols (Middleton, 2017), the concern about Saharan intrusions might be reduced from a public health point of view in this context, since air masses from European and regional origin were more enriched in human-associated bacteria than Saharan air masses. A more detailed boxplot representation of the distribution of the 116 more abundant genera for the four air mass origins is reported in Supplementary Material (Fig. SM5).

Correlations between microbiology and chemistry

Chemical and microbiological data were combined to check possible correlations between the variables and the sample provenances. In particular, some typical markers of Saharan dust, biomass burning, and industrial activities, the two latter being particularly enriched in non-Saharan samples, were identified among the chemical analytes. Moreover, the analysis of the β diversity showed that the microbial communities of long-range transported Saharan dust were significantly different from those sampled when other air masses were present, strongly supporting the hypothesis that desert dust can impact the bacterial composition of the aerosol at our latitudes (Gat et al., 2017; Mazar et al., 2016; Rosselli et al., 2015). On the contrary, non-Saharan samples showed similar communities among each other, which in fact clustered together (Fig. 2, Fig. 3). Nevertheless, as observed for the chemical characteristics, even if similarities existed within the PM samples sharing the same origin, the differences were not negligible and suggested that each event was independent of the others. This has been already observed in previous works, at least for dust events. In fact, significant differences in bacterial community structure were reported during different dust events that impacted the same area, even when two events were very close in time (Federici et al., 2018; Yamaguchi et al., 2014).

In order to combine information about microbiology and chemistry, the NMDS1 dimension from the statistical analysis on bacterial communities was correlated with the concentration of chemical variables, normalized against the PM values (w/w). Some of the statistically significant correlations (p = ***) are shown in Fig. 5 .

Fig. 5

Scatterplots of the significant relationship between one-dimensional non-metric multidimensional scaling (NMDS1) and the concentration (m/m in PM10) of some chemical variables. Dashed lines are linear correlations.

Specifically, a significant correlation was found between NMDS1 and PM2.5 (i.e. PM2.5/PM10 ratio). PM ratio was lower for Saharan dust with the exception of the outlier SH_20141107 which corresponded to the weakest Saharan dust event, with a PM10 concentration of 8.4 μg/m3. The correlation was significant also between NMDS1 and PMcoarse., which was higher for SH because Saharan intrusions are constituted by coarser particles. Organic carbon (OC) content correlated significantly with NMDS1. OC in Saharan dust was lower than in non-Saharan samples because the latter can have a higher anthropogenic contribution. The highest OC/PM10 values were found for SH_20140624 (which had also high sulphate concentration) and SH_20140522 of NW provenance. As stated above, also many bacterial genera that were significantly more abundant in non-SH than in SH samples (e.g. Lactobacillus, Streptococcus, Propionibacterium and Haemophilus) were generally related to anthropic and built environments. This confirms the relevance of the impact that densely populated areas may exert on bacterial populations transported by air masses.

Anthracene was the only PAH correlating significantly with NMDS1, being higher for SH samples. The sum of low molecular weight PAHs (LW) was also higher for SH samples. Calcium concentrations showed no correlation with NMDS1, which was interpreted as due to the high local contributions of this element. Iron, on the other side, was richer in SH samples and correlated negatively with NMDS1. Ammonium and sulphate concentrations were generally higher for non-Saharan air masses. Innocente et al. (2017) also reported that high ammonium and sulphate concentrations were associated with long-range transport from north-west in Milan (North Italy), and those air masses presented a high percentage of Propionibacterium. This is in agreement with our SIMPER analysis, which indicated the genus Propionibacterium as significantly more abundant in NW than in SH samples (Fig. 4). However, Innocente and colleagues also reported that this correlation was weak, and ionic composition of air masses was much more clearly related to air mass provenience than to bacterial community structure.

Indeed, also in our case, since NMDS1 was strongly correlated to air mass origin, it was not possible to fully understand whether variations in chemical variables were more correlated to variations in bacterial community structure or to dust provenance.

Conclusions

In this work, we have characterized the bacterial communities of 19 air masses of different origin, sampled as PM10 at the remote site of Monte Martano, in Central Italy. This EMEP station is representative of the Central Mediterranean area. The main results of the present work can be summarized as follows:

Four distinctive air masses were identified: previous similar work on this topic was substantially limited to Saharan (SH) dust air masses while in the present study we extended the characterization to regional (RG), North-Western (NW), and North-Eastern (NE) air masses.

At genus level, the distribution of the bacterial populations in air masses clearly showed differences due to the sample provenance. In fact, PM10 during Saharan intrusions showed a relatively low number of genera, while non-Saharan samples, particularly those collected during regional movements of air masses (RG), showed a high number of different genera with low abundance, indicating a highly diverse and even community.

The abundance of bacterial genera including potential human and animal pathogens is higher in non-Saharan than in Saharan samples and stressed the relevance of anthropic impact on bacterial populations transported by air masses that crossed densely populated areas.

A significant correlation was found between the NMDS1 dimension from the statistical analysis on bacterial communities and specific chemical variables determined on the samples. Particularly, organic carbon and ammonium concentrations, which are markers of anthropogenic contribution, were found to increase along the sequence SH, NW, NE, and RG, and be also correlated to an increasing abundance of human and animal associated bacteria. However, the strong correlation between bacterial community structure and air mass origin hampered the possibility to disentangle the effects of variations in bacterial populations and of dust provenance on variations in chemical variables.

CRediT authorship contribution statement

CP, RS, BS, BM: chemical analysis, GLP: statistical analysis, SC: Atmospheric modelling: AF, IG, EF, EM, EC, CC: microbiological analysis:, DC: Conceptualization, Methodology, Data curation, Original draft preparation. CP, IG and DC, Writing and editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.