Alterations in the airborne bacterial community during Asian dust events occurring between February and March 2015 in South Korea

During Asian dust events, a relatively high concentration of particulate matter is transported by wind from arid and semi-arid regions, such as the Gobi and Taklamakan deserts, to nearby countries, including China, Korea, and Japan. The dust particles contain various microorganisms, which can affect human health as well as the environmental microbe population. In the current study, we investigated the characteristics of the airborne bacterial community during Asian dust events between February and March 2015 in South Korea. Bacterial diversity indexes such as operational taxonomic units, Chao1 and Inverse Simpson index were increased, along with total 16S rRNA gene copy number during Asian dust events. The bacterial community structure during Asian dust events was clearly distinguishable from that during non-Asian dust days. The genera Bacillus and Modestobacter were increased 3.9- and 2.7-fold, respectively, while Escherichia-Shigella was decreased by 89.8%. A non-metric multidimensional scaling plot with metadata analysis revealed association of particulate matter concentration, but not temperature, humidity or wind speed, with bacterial community structure, suggesting that the newly transported dust particles contain various microorganisms that influence the airborne bacterial environment.

environmental and technical problems. (i) Since AD events occur over only a few days every year, it is difficult to collect sufficient samples to perform effective statistical analysis. (ii) There are no arrangements in place to collect airborne samples for elucidation of biological features on ordinary days. (iii) The existing culture and PCR-based methods for analyzing bacterial structures, including the plate method, Sanger sequencing, and DGGE, are not suitable for handling unculturable bacteria and multiple data on a large scale 18,20,21 .
In the current study, we monitored and collected airborne dust particles from January to March 2015 in South Korea, and obtained daily samples, including 9 days of AD events and 11 non-AD days. Fortunately, AD events were relatively frequent during this period, allowing us to secure sufficient samples for comparison with normal dust-free days. The whole metagenomics approach using next-generation sequencing with the Roche 454 system was applied for analysis of more than 200,000 sequences.
The main aim of this study was to determine the bacterial community structures of the airborne environment during AD events in Korea and analyze the information obtained using a metagenomic approach. Establishment of bacterial characteristics through continuous monitoring of the airborne environment during AD events may aid in understanding the biological features of dust particles as well as management of associated environmental impacts and human health problems.

Results
Environmental conditions and bacterial diversity analysis. To determine the airborne bacterial communities during AD events, we collected dust samples using a high-volume sampler from January to March 2015 from the rooftop of the university building in Goyang-si, South Korea. During the sampling period, no precipitation occurred to limit the collection of dust particles. The humidity of the sampling period varied between 26.6-71.3%, temperature from − 9.9 to 10.4 °C, and wind speed from 1 to 4.5 m/s. From February to March 2015, nine days of AD events samples were collected. We examined samples from a total of 20 days, including 11 non-AD days. As shown in Supplementary Table S1. Using the classification data, we confirmed operational taxonomic units (OTUs) at the species level during both AD event and non-AD days. As shown in Fig. 1a and Supplementary Table S1, OTU results varied from around 68 (March 10, non-AD) to 733 (February 8, AD; March 19, non-AD), with the average value for AD samples being significantly higher than that for non-AD samples (P-value = 0.019, independent t-test). In addition, the Chao1 index, but not inverse Simpson index (representing richness and evenness diversity of the bacterial community, respectively) was significantly different between AD and non-AD samples (Fig. 1b and Supplementary Table S1, P-value = 0.029, independent t-test). These results suggest that airborne bacterial diversity (at least the richness index) is increased during AD events, and transported dust particles contain various microorganisms that can alter an airborne bacterial community.
16S rRNA gene copy numbers in the collected samples. Since newly transported dust particles during AD events were shown to contain various bacteria that could influence airborne bacterial diversity, we measured the 16S rRNA gene copy number in samples to estimate the relative amounts of different bacterial communities in the airborne environment. For quantification, E. coli ATCC 11775 was used as the standard strain containing one copy of the 16S rRNA gene in the genome. Using the standard curve of real-time PCR results based on the E. coli strain, the 16S rRNA gene copy number was determined in samples during both AD events and non-AD days. During AD events, a relatively higher concentration of particulate matter 10 (PM10), presenting dust particles under 10 μ m in size, was observed in the airborne environment. As expected, the 16S rRNA gene copy number was correlated with PM10 concentration (Spearman correlation coefficient, 0.6496; P-value = 0.00247), and both factors were significantly higher in AD samples, compared to non-AD samples (independent t-test, P-value = 0.047 and 0.027, respectively) (Fig. 2). In particular, the average 16S rRNA gene copy number was 20 times higher in AD event samples, with the highest 16S rRNA gene copy number (AD sample, March 21) being ~8,000 times greater than the lowest 16S rRNA gene copy number in a non-AD sample (February 26). Together with diversity index data, these results support an increase in the number of microbes transported with dust particles during AD events, which can influence the local airborne bacterial community. Non-metric Multidimensional Scaling (NMDS) analysis. To elucidate whether characteristics of the bacterial community are clearly distinguishable based on transported microorganisms during AD events, we analysed the differences in bacterial community structures between AD and non-AD days using non-metric multidimensional scaling (NMDS) at the species level. Among the nine days of AD events, samples from six days are presented along the left direction of axis1, and samples from most days of non-AD days presented along the right direction of axis1 (Fig. 3). We observed no differences between AD and non-AD samples along axis2. The bacterial community structures of dust collected from several AD events (March 1, 16, and 17) were presented in the right direction of axis1 with non-AD days. The results obtained for three days of AD events were correlated with the Chao1 diversity index, which presented a relatively low level of diversity, compared to other AD samples. The differential distribution of bacterial structures between AD and non-AD samples on the NMDS plot was statistically confirmed using analysis of molecular variance (AMOVA) (P-value = 0.026). To further determine the specific factors affecting distribution on the NMDS plot, we calculated the correlation coefficient between the direction of axis and experimental or environmental factors, including 16S rRNA gene copy number, PM10 concentration, temperature, wind speed, and humidity. As shown in Table 1, only two factors (16S rRNA gene copy number and PM10 concentration) were significantly correlated with the negative direction along axis1 of NMDS plots (Correlation coefficient = 0.55 and 0.46; P-value = 0.0069 and 0.029). These results suggest that the airborne bacterial community structure was altered through influx of newly transported microorganisms, but not by general meteorological conditions such as temperature, wind speed, and humidity, during the AD events under investigation. Construction of a phylogenetic tree with the ThetaYC algorithm at the species level additionally revealed considerably different bacterial community structures of AD and non-AD samples   Figure S1, both weighted and unweighted significance values between AD and non-AD days were less than 0.05). Our data indicate that the airborne bacterial community is clearly distinguishable between AD and non-AD episodes owing to newly transported bacteria, with no contribution from meteorological elements. Although the airborne microbial environment was gradually restored through atmospheric flow within a few days, continuously occurring AD events and the consequent transient increase in microorganisms may contribute to several human diseases and exert widespread effects on the ecological system.

Analysis of the bacterial community during AD and non-AD periods at the phylum level.
To identify the bacterial community in the airborne environment between February and March 2015, pyrosequencing results were classified using Mothur software with the SILVA database 22,23 . During the non-AD period, the predominant phylum in the airborne environment was Proteobacteria, followed by Firmicutes, Actinobacteria, Deinococcus-Thermus, Acidobacteria, Gemmatimonadetes, and Chloroflexi ( Fig. 4). Proteobacteria remained the dominant phylum during AD events, but the proportion decreased from 69% to 53.5%. Notably, the relative amounts of several other phyla, including Actinobacteria, Firmicutes, and Acidobacteria, were increased during AD events. Next, we performed a Metastats analysis to identify the phyla displaying significant differences between AD and non-AD periods ( Table 2) 24 . Our data revealed a significant decrease in proportion (by 22.42%) of the predominant phylum, Proteobacteria, in AD samples. Simultaneously, Actinobacteria, Acidobacteria, Gemmatimonadetes, and Chloroflexi were increased 1.9, 2.2, 2.8, and 3.2 fold, respectively, in AD, compared to non-AD samples. In particular, the phylum Actinobacteria displayed the highest proportion (21.53%) except Proteobacteria during the AD period, consistent with previous data obtained in December 2014 in Seoul, Korea (manuscript under submission).

Analysis of the bacterial community during AD and non-AD periods at the genus and species levels.
To elucidate additional alterations in the airborne bacterial community during AD events, we analyzed the classification data at the genus and species levels (Tables 3 and 4). To this end, we focused on selected genera and species representing > 0.5% of the total community. Analysis at the genus level revealed significant differences in the proportions of 21 genera, compared to non-AD samples. Moreover, all three genera belonging to the phylum Proteobacteria were decreased in AD samples. In particular, the proportions of Escherichia-Shigella and Pseudomonas were decreased by 89.8% and 52%, respectively. Although the proportion of genera showing increased abundance were relatively low within the entire classified bacterial community, most genera displayed a two to four fold increase in AD events, compared to non-AD periods. In particular, the proportion of the genus Geodermatophilus belonging to the phylum Actinobacteria was 6.42 times higher in AD than non-AD samples. Additionally, genera Microvirga, Methylobacterium, Rubellimicrobium, Gemmatimonas and Modestobacter  Table 1. Calculation of the correlation coefficient between the nonmetric multidimensional scaling plot and experimental or environmental factors. P-values < 0.05 mean that the indicated factor is correlated with the direction of plot axis. 16S rRNA gene copy number and PM10 concentration, but not other factors, were significantly correlated with the negative direction along axis1 of the plots. presented relatively high proportions during AD events. In contrast to the Metastats analysis at the phylum level, the genus Bacillus belonging to the phylum Firmicutes was increased 3.93-fold in AD samples, constituting the highest proportion of significantly altered genera during AD events (from 1.9% to 7.6%). For analysis of the bacterial community at the species level, closely related species after classification of pyrosequencing data were selected based on the SILVA database. As shown in Table 4, among the species showing significant differences, Bacillus circulans displayed the highest increase in quantity (8.84-fold) and proportion (3.14%) in AD samples, followed by Methylobacterium iners, Sphingomonas starnbergensis, Micrococcus terreus, Rubellimicrobium roseum, and Rubellimicrobium aerolatum. In contrast, the proportions of Escherichia coli, Escherichia fergusonii, Sphingomonas echinoides, and Sphingomonas oligophenolica were decreased by ~90%. Despite slightly inaccurate classification at the species level (since the results were obtained using a limited length of the complete 16S rRNA gene sequence, taken as a whole), the genera Bacillus belonging to the phylum Firmicutes, Escherichia belonging to the phylum Proteobacteria, and Modestobacter belonging to the phylum Actinobacteria were clearly distinguished among the airborne bacterial communities during AD and non-AD periods.

Discussion
In this study, we confirmed an increase in the total bacterial levels and community diversity in the airborne environment during the AD events in Korea. Moreover, the bacterial population was clearly distinguishable between AD and non-AD periods. Since AD events represent a natural phenomenon of transportation of dust  Table 3. Metastats analysis at the genus level. The table contains data on significantly different genera (P-value < 0.05) with a proportion greater than 0.5% of the total bacterial community. a Mean values signify the average proportion of the genera. b P-values < 0.05 indicate significant differences in proportion between the two groups (AD events and non-AD days).
particles from their region of origin and dust contains various microorganisms, we expected the samples collected during AD events to include newly transported, together with pre-existing microorganisms. As predicted, we obtained significantly higher values of diversity indices, including OTUs and Chao1, in AD samples, compared with non-AD samples, consistent with previous findings on dust particles [25][26][27] . Further examination of inverse Simpson index, which presents a richness feature with evenness diversity, confirmed that the evenness diversity index does not distinguish between AD and non-AD samples, in contrast to OTU results and Chao1 index, indicative of richness diversity. A previous report by Jeon et al. 27 also showed that the Chao1 value of AD samples was 2.6 times higher than that of non-AD samples, while evenness values of each sample were determined as 0.851 and 0.887, respectively. Since evenness diversity reflects the number as well as proportion of specific OTUs, it is not surprising that this index is not altered by AD events. An et al. 28 reported that richness values for control and AD samples are not correlated, while Mazar and co-workers proposed that the results are attributable to different sampling methods 26 . With the method of sampling, various conditions, including region and period of sampling, an origin of AD event, and meteorological conditions, may affect the biological characteristics of AD and non-AD samples and therefore consider for comparison. Therefore, it is important to clarify the features of AD events from various viewpoints and accumulate biological information consistently to obtain an overall picture of microbial alterations. The airborne bacterial populations between AD events and non-AD periods were clearly distinguished using NMDS plot and phylogenetic tree analysis. In our experiments, PM10 concentration and 16S rRNA gene copy number, which were significantly correlated with AD events (Fig. 2), but not temperature, wind speed, and humidity, contributed to different bacterial community structures, suggesting that the differences between AD and non-AD samples are attributable to newly transported dust particles containing microorganisms. Several geochemical elements, including SO 2 and NO 2 , have been correlated with PM10 concentration 29 . However, our data, in conjunction with other previous studies 30,31 indicate that geochemical elements, including SO 2 , NO 2 , CO and O 3 , are not correlated with PM10 concentration, and do not differ significantly between AD and non-AD conditions during the period of sampling (independent t-test, P-values for SO 2 , NO 2 , CO, and O 3 were 0.1918, 0.4348, 0.2143, and 0.9675, respectively., Supplementary Table S2). Information on airborne bacterial environments, together with geochemical and geophysical elements, may be valuable in elucidating the characteristics of dust particles between AD events and non-AD periods.
Finally, based on analysis of airborne bacterial communities from both AD and non-AD samples, we established that among the significantly altered genera, Bacillus and Modestobacter belonging to the phyla Firmicutes   Table 4. Metastats analysis at the species level. The table contains data on significantly different species (P-value < 0.05) with a proportion greater than 0.5% of the total bacterial community. a Mean values signify the average proportion of the species. b P-values < 0.05 indicate significant differences in proportion between the two groups (AD events and non-AD days).
and Actinobacteria, respectively, had the highest proportions during AD events between February and March 2015. Consistent with our findings, increased proportions of Firmicutes and Actinobacteria in AD samples have been reported in previous studies 19,32,33 . The genus Bacillus has been isolated from soil, sediment, or sewage, and produces endospores, which may be present in the newly transported dust particles. Bacillus circulans is reported to cause sepsis, bacteremia, abscesses, and meningitis in humans 34,35 , suggesting that newly transported dust particles during AD events containing this microorganism could contribute to several human diseases. The species belonging to the genus Modestobacter is also found in soil or sediment, but does not produce endospores. However, common features of the genus Modestobacter, such as resistance to desiccation and irradiation 36 , facilitate survival in the desert, making it one of the major microorganism constituents of dust particles during AD events. At the time of writing the manuscript, four species of the genus Modestobacter had been identified. At present, there is no evidence to suggest a relationship between this bacterium and human disease. Interestingly, the proportion of the genus Escherichia-Shigella was largely decreased in AD samples (11.3% to 1.2%), consistent with results from a previous study 28 . Yamaguchi et al. also reported a relatively low proportion of the class Gammaproteobacteria containing the genus Escherichia-Shigella in Asian dust samples collected from Japan as well as its regions of origin, such as the Taklamakan and Gobi deserts 32 . Based on the 16S rRNA gene copy number results, which revealed that several hundred to thousand times more bacteria are present in AD samples (Fig. 2), we propose that the genus is rarely transported from its origin and present at a relatively low proportion in AD samples. The reasons underlying the low proportion of Escherichia-Shigella in dry or semi-dried regions remain to be established.
In conclusion, an in-depth investigation of the airborne bacterial community during AD events in this study confirmed that total amounts of bacteria in the airborne environment as well as their richness diversity are significantly altered, compared to the non-AD environment. Continuous monitoring of the characteristics of airborne bacterial community structures should be performed to obtain an overview of results that may provide important insights into the biological features of transported dust particles during AD events.

Materials and Methods
Collection of air samples. During  DNA extraction and pyrosequencing. For genomic DNA extraction, 10 ml PBS was used to suspend collected dust particles on a 20 × 100 mm membrane. After resuspension, dust particles containing microorganisms were precipitated via centrifugation at 10,000 g for 1 min. Bacterial cells were lysed with DNA extraction buffer (100 mM Tris-Cl, pH 8.0; 100 mM sodium-EDTA, pH 8.0; 100 mM sodium phosphate, pH 8.0; 1.5 M NaCl; 1% CTAB) and 1 mg of proteinase K, followed by incubation at 37 °C for 30 min and subsequent addition of 2% SDS. After incubation at 65 °C for 30 min, cell debris and dust particles were precipitated via centrifugation. The supernatant fractions were mixed with an equal volume of chloroform and centrifuged. The upper layer was mixed with 70% v/v isopropanol and genomic DNA precipitated. After washing with 70% ethanol, the genomic DNA pellet was dried and resuspended with distilled water. For preparation of pyrosequencing samples, the 16S rRNA gene region was amplified with a set of primers containing the fourth and fifth hypervariable regions, according to a previous report 37 . The forward primer sequence was 517 F of the 16S rRNA gene (GCCAGCAGCCGCGGTAAT) and the reverse primer was 896 R (CCGTACTCCCCAGGCGG). For pyrosequencing, the forward and reverse primers were conjugated with adaptor sequences and multiplex identifier tag sequences provided by the manufacturer (Macrogen, Korea). Amplification of 16S rRNA gene was performed as follows: one cycle of pre-denaturation at 95 °C for 3 min, 35 cycles of denaturation at 95 °C for 3 s, and annealing and extension at 64 °C for 15 s. After amplification and purification (GeneAll Biotechnology, Korea), equal amounts (30 ng) of each PCR product were collected in a tube, and pyrosequencing performed with the GS FLX system (Roche, Switzerland).

Quantification of 16S rRNA gene copy number.
To determine the 16S rRNA gene copy number of each sample, E. coli strain ATCC 11775 was used to prepare a standard curve. Briefly, serially diluted E.coli strains were inoculated onto tryptone soy agar plates followed by counting growth colonies, and real-time PCR performed using an equal amount of E. coli as the template with the following primer set and amplification conditions: forward primer 517 F (GCCAGCAGCCGCGGTAAT), reverse primer 805 R (GACTACCAGGGTATCTAATCC), one cycle of pre-denaturation at 95 °C for 10 min, 30 cycles of denaturation at 95 °C for 3 sec, and annealing and extension at 60 °C for 15 sec. With the result of counted colonies, a standard curve was constructed using the threshold cycles of real-time PCR results. Based on the standard curve, the 16S rRNA gene copy number was assessed from real-time PCR data obtained using the samples collected as template. Real-time PCR was performed using the same conditions and primer sets as above.
Scientific RepoRts | 6:37271 | DOI: 10.1038/srep37271 Data analysis. Pyrosequencing results were analyzed using standard Mothur software 22 . Briefly, all low-quality sequences, including small fragments (< 200 bp) and chimeric sequences, were removed using Mother and UCHIME software 38 , respectively. Various sequence numbers (2,445 to 14,771) were obtained after trimming low-quality data of pyrosequencing results and 2,445 sequences of subsamples randomly selected for normalization of read numbers. After trimming the sequences, the SILVA database was used for classification of the bacterial community 23 , and statistical analysis of the classified results performed using Metastats software 24 . Following classification, mitochondria, chloroplast, archaea and eukaryote-associated sequences were removed. For diversity analysis, we used a phylotype analysis method based on classified taxonomic information rather than cutoff values (using Mothur command, phylotype), since several sequences classified into the same taxonomy presented different OTUs using a distance-based cutoff method. The airborne bacterial community was analyzed at the phylum, genus, and species levels using taxonomic information, and bacterial diversity indices and NMDS plots computed using phylotype results at the species level. The correlation coefficient between NMDS plot and meteorological elements was calculated via metadata analysis with Mothur software. Statistical values, including Pearson and Spearman correlation coefficients and independent t-tests, were calculated using R software with the Rcmd package 39 .