Precipitation samples were collected at two sites in the Greater Tokyo Area, Japan: Kashiwa (35°54′00″N, 139°55′59″E, 50 m above sea level) and Hongo (35°42′55”N, 139°45'56″E, 30 m above sea level) (Figure 1). The Kashiwa site was on the roof of a seven-story building on the Kashiwa campus, the University of Tokyo, Chiba, Japan, which is surrounded by residences, farms, and woods in a suburb of Tokyo. The Hongo site was on the roof of a five-story building on the Hongo campus, the University of Tokyo, Tokyo, Japan, which is located in downtown Tokyo. The sites are 25.5 km apart and neither geologically nor meteorologically separated in the Kanto plain. The upper areas of both sites are wide open and lack any obstructing buildings or structures that would contaminate the precipitation samples. At the Kashiwa site, precipitation was aseptically collected using a US-330 automatic precipitation sampler (Ogasawara Keiki, Tokyo, Japan) following the method of Kaushik et al. (2012). This device consists of a sterile and disposable bottle inside a 4°C refrigerator and automatically collects precipitation by opening the lid only when a sensor detects precipitation. At the Hongo site, precipitation samples were manually collected into a sterile and disposable bottle on ice and immediately stored in a 4°C refrigerator. At both sites, every part of the collection equipment that potentially directly contacted precipitation samples (e.g., disposable collection bottles and channel tubes) was sterilized by gamma rays in advance of each sample collection. The precipitation samples were pre-filtered through 5-μm membrane filters, and microbial cells were collected using 0.22-μm Sterivex filters (Millipore, USA). The Sterivex filters were promptly moved to a −20°C freezer and stored until DNA extraction. Precipitation sampling required no special permission. To prepare negative control samples, we poured 1 L of Milli-Q purified water into the collection equipment and carried out filtration in the same manner.

We collected 25 and 5 precipitation samples containing sufficient amounts of microbial DNA at the Kashiwa and Hongo sites, respectively. The sampling dates spanned more than 1 year from May 2014 to October 2015, encompassing the rainy and typhoon seasons in Japan (Table 1; the six digits, letter, and suffix number for each sample name represent the sampling date (YYMMDD), the site (K for Kashiwa and H for Hongo), and the volume (if multiple samples were collected during the same precipitation event). A precipitation event was defined if there was no precipitation 6 h before and after the event. The volumes of collected and filtered precipitation ranged from 50 to 1,000 mL. For correlation analysis with meteorological data, we excluded the data obtained from samples 140630K_50, 140630K _100, 140810K_50, and 140810K_100, which were retrieved as replicate samples with different volumes. Eight negative control samples were also collected at different dates at the Kashiwa and Hongo sites.

Microbial DNA on the Sterivex filters was retrieved using a ChargeSwitch Forensic DNA Purification Kit (Invitrogen) according to the supplier's protocol with one exception: the filters were directly suspended in the extraction solution from the kit during the cell lysis process. The V5-V6 region of the prokaryotic 16S rRNA gene was amplified using a standard PCR protocol with TaKaRa Ex Taq (TaKaRa) and the following high-performance liquid chromatography-purified primers: 784F (5′- RGGATTAGATACCC -3′) and 1064R (5′- CGACRRCCATGCANCACCT -3′) (Wang and Qian, 2009; Claesson et al., 2010). Amplified DNA was concatenated to multiplex identifier tags that were unique to each sample, and a mixture of 10 samples on average was sequenced in one run on a 454 GS Junior System (Roche) after size selection (350 ± 50 bp). Pre-packaged sterile water for injection (in lieu of water from a laboratory water purification system) was used throughout the DNA extraction, PCR amplification, and DNA library preparation steps to avoid water-mediated contamination.

For raw sequence data from both precipitation and negative control samples, sequence regions at both ends that contained low-quality bases (quality score < 20) were trimmed using DynamicTrim (Cox et al., 2010), chimeric sequences were filtered out using UCHIME with default settings (Edgar et al., 2011), and sequences whose lengths were shorter than 150 bp were discarded. All remaining high-quality sequences were clustered with a 97% identity threshold using CD-HIT (Fu et al., 2012). After discarding clusters that contained negative control sequences (Cho and Jang, 2014), each cluster was designated as an operational taxonomic unit (OTU). For hierarchical cluster analysis of the precipitation samples, the Ward method was used based on Bray-Curtis dissimilarities between their OTU compositions. Non-metric multidimensional scaling (NMDS) analysis was conducted using Bray-Curtis dissimilarities. The taxonomic assignment of each OTU was performed by conducting a BLASTN search (Camacho et al., 2009) against the SILVA database (Quast et al., 2013) and retrieving the top hit sequence that showed e-values ≤ 1E-15. To estimate ordinary habitats for each 16S rRNA sequence, a BLASTN search was performed against MetaMetaDB (Yang and Iwasaki, 2014), and the top hit sequence with an e-value ≤ 1E-10 and an identity ≥ 90% was retrieved. Microbial habitability index (MHI) scores were calculated as previously described (Yang and Iwasaki, 2014).

Amplicon-sequencing data of aerosol and cloud water samples were downloaded from NCBI SRA database (the accession numbers are shown in Supplementary Table S1). Their ordinary habitat analyses were conducted as described above after quality filtering.

The data on the amount of precipitation, temperature, wind speed, and atmospheric pressure were retrieved from the website of the Japan Meteorological Agency (http://www.jma.go.jp/jma/menu/menureport.html), Ministry of Land, Infrastructure, and Transport of Japan. Precipitation, wind speed, and temperature data from the Abiko (35°51′48″N, 140°06′36″E; 16.4 km from Kashiwa) and Tokyo (35°41′30″N, 139°45′00″E; 2.9 km from Hongo) observatories were used for the analyses of the Kashiwa and Hongo sites, respectively (Figure 1). Atmospheric pressure data from the Tokyo observatory were used (this observatory is the closest to both sites that records atmospheric pressure data). The wind speed, temperature, and atmospheric pressure data were averaged over the period of each precipitation event. To analyze long-range transport paths of air masses that caused precipitation by providing water vapor, we estimated backward trajectories of an air mass at 2,000 m altitude for 240 h prior to all precipitation events for each sampling site. The trajectories were calculated based on the hybrid single-particle Lagrangian integrated trajectory (HYSPLIT) model (http://ready.arl.noaa.gov/HYSPLIT.php) provided by the Global Data Assimilation System of National Oceanic and Atmospheric Administration, USA (Stein et al., 2015). The HYSPLIT model uses gridded meteorological data and considers advection and diffusion of air parcels in calculation of their trajectories. This model has been used in a variety of atmospheric simulations focusing on the atmospheric transport, dispersion, and deposition of pollutants and hazardous materials (Stein et al., 2015), while it has also been adopted for estimation of sources of airborne microbes (e.g., Smith et al., 2013; Cho and Jang, 2014; Kobayashi et al., 2015; Xia et al., 2015; Xu et al., 2017).

The amplicon sequence data were deposited in the DDBJ/ENA/GenBank database under BioSample IDs SAMD00059586-SAMD00059614 and SAMD00060461-SAMD00060468. All data were registered under BioProject ID PRJDB5087.

A total of 64,100 high-quality sequences 231 ± 45 bp in length were generated from 30 precipitation and eight negative control samples. The precipitation samples included typhoon rain, rainy season rain, and snow. After removing sequences exhibiting >97% similarity to the negative control samples, 12,089 “effective” sequences comprising 1,297 OTUs remained. To make our analyses based on reads that were not likely from contamination as much as possible, we took a conservative and strict filtering approach, whose extent of read number reduction was similar to that in a previous study (Cho and Jang, 2014). The number of OTUs per sample ranged from 4 to 226 (Table 1). Based on rarefaction curves, the obtained OTUs represented their microbial communities well for some samples, although several samples required additional sequences (Supplementary Figure S1).

Hierarchical cluster analysis of OTU composition in the precipitation samples indicated samples collected during the same precipitation event with different volumes (50, 100, and 200 mL) that were highly similar to each other (Figure 2, open symbols), suggesting that differences in volume have little effect on analysis in the 50–200 mL range. Moreover, microbial communities in samples that were collected on the same day at different sampling sites (Kashiwa and Hongo) were closely positioned in the dendrogram (Figure 2, closed symbols), indicating that the observed OTU compositions reflect the microbial populations in precipitation rather than those in the atmosphere near the ground surface or equipment- or reagent-mediated contamination at each site. NMDS analysis did not show any clear trend, although samples of close dates tended to be clustered together (Supplementary Figure S2).

Among the 12,089 effective sequences, 11,994 (99.2%) were taxonomically assigned at the phylum level. Almost all sequences were assigned to 24 phyla in the domain Bacteria with the exception of 4 (0.03%) and 219 (1.7%) sequences assigned to Archaea and mitochondria, respectively. This strong bias toward bacterial sequences may reflect the actual composition but may also be attributable to amplification bias introduced by primer specificity. The top three and six most abundant bacterial phyla accounted for >80 and >95%, respectively, of the sequence pool of all precipitation samples (Figure 3A). Proteobacteria was the most abundant phylum (23–88%) across all precipitation samples with the exception of the 140630, 140926, and 150116K samples (Firmicutes (89–94%), Actinobacteria (50%), and Firmicutes (49%) were the most abundant phyla, respectively). A particularly exceptional microbial community dominated by Firmicutes was observed in the 140630K sample. Firmicutes, Bacteroidetes, and Actinobacteria were the other dominant phyla in the total sequence pool. In principle, these results were consistent with those of a previous study in which Proteobacteria, Firmicutes, and Bacteroidetes were the dominant phyla in precipitation samples captured in Seoul, Korea (Cho and Jang, 2014), whereas comparatively greater numbers of sequences were assigned to Actinobacteria, Planctomycetes, and Cyanobacteria in this study. At the class level, the abundant groups were Gammaproteobacteria, Betaproteobacteria, and Alphaproteobacteria, followed by Bacilli, Flavobacteriia, Clostridia, Actinobacteria, and Sphingobacteriia (Figure 3B). Notably, the enrichment of these phyla and classes was also reported in previous studies investigating aerosolized (Bowers et al., 2009, 2011a; DeLeon-Rodriguez et al., 2013) and cloud water microbial communities (Kourtev et al., 2011; DeLeon-Rodriguez et al., 2013).

Several OTUs were assigned to genera that potentially contain INA bacteria, i.e., Acinetobacter, Bacillus, Erwinia, Flavobacterium, Luteimonas, Microbacterium, Pseudomonas, Psychrobacter, Sphingomonas, and Stenotrophomonas (Després et al., 2012). We also detected several genera containing known pathogens, including typical human pathogens such as Legionella, Streptococcus, Arcobacter, Rickettsia, and Clostridium, and plant pathogens such as Erwinia, although their abundance was low. We did not detect season-specific microbial groups in the typhoon rain, rainy season, and snow samples with statistical significance, probably partly due to small sample sizes.

Taxonomic distribution exhibited seasonal variability (Figure 3). Notably, the abundance of Proteobacteria decreased from summer to winter (p < 0.01, Mann-Whitney U-test), and a similar trend has consistently been observed in aerosolized microbial communities (Bowers et al., 2011b). To more closely investigate the factors underlying changes in the precipitation microbial communities, we performed a correlation analysis between meteorological characteristics and microbial composition (Figure 4). The relative abundance of the order Bacteroidales negatively correlated with temperature (Spearman correlation ρ = −0.70, p < 0.01 after the Bonferroni correction). Although other correlations were not statistically significant after multiple testing correction, the amount of precipitation, wind speed, and atmospheric pressure showed tendencies of positive correlations with the abundance of the orders Cellvibrionales (ρ = 0.59), Cellvibrionales (ρ = 0.58), and Pseudomonadales (ρ = 0.57), respectively. Notably, the abundance of the order Legionellales, which contains several known pathogens, showed a tendency of a positive correlation with temperature (ρ = 0.47), where aerosolized water is known to facilitate the dispersion of Legionella (Nguyen et al., 2006) and a warm and wet climate is associated with the incidence of Legionnaires' disease (Fisman et al., 2005; Fisman, 2007). Although cell numbers were not measured except for one sample in this study, we note that seasonal variability in cell numbers would also be important, especially because that of atmospheric samples was reported (Kaushik et al., 2014; Dong et al., 2016). Similarly, analyses with particulate matter density and O₃ and NO₃ concentrations are also envisioned, because they would substantially affect aerial microbes (Kaushik et al., 2012; DeLeon-Rodriguez et al., 2013; Wei et al., 2016; Xu et al., 2017).

To estimate the environments from which microbes in precipitation originated, we performed a microbial habitat index analysis using MetaMetaDB (Yang and Iwasaki, 2014), which is a database to estimate the ordinary habitats of microbes based on similarity searches for 16S rRNA gene sequences against amplicon-sequencing and shotgun metagenomic data in public databases. In most samples, animal-associated environments, such as gut microbiota, were estimated to be the most dominant ordinary habitats (52% on average) (Figure 5, Supplementary Figure S4), which is consistent with a previous study in which animal feces were the dominant source of airborne bacteria (Bowers et al., 2011b). Notably, marine-related environments, such as marine and marine sediment, were estimated to be relatively major ordinary habitats for several samples (e.g., 65.1 and 63.1% in the 140810 and 141014K samples, respectively). Soil-related environments, such as soil and rhizosphere, were also estimated to be major ordinary habitats (11.0% on average). For comparison, we also conducted ordinary habitat analyses using amplicon-sequencing data from aerosol (Xia et al., 2015) and cloud water (DeLeon-Rodriguez et al., 2013) samples. The soil-related and animal-associated environments were generally major ordinary habitats as consistent to the present results, whereas marine-related environments were not major possibly because the origins of the microbes or the sampling methods were different from those in our study (Supplementary Figure S5).

The estimated backward trajectories of air masses that led to the precipitation events at the Kashiwa and Hongo sites were classified as terrestrial, oceanic, and hybrid routes. The terrestrial route typically originated from the middle of the Eurasian continent and passed through the East China Sea, the Yellow Sea, and the Sea of Japan; the oceanic route typically originated from the Pacific Ocean and passed through the East China Sea or the Sea of Okhotsk; and the hybrid route comprised both the terrestrial and oceanic areas. Consistent with the typical pattern of the seasonal winds in Asia, the terrestrial and oceanic routes dominated in winter and summer, respectively (Figure 5, Supplementary Figure S3). The estimated ordinary habitats of the precipitation microbes showed agreement with the estimated air mass backward trajectories. For example, Planctomycetes, which contains several aquatic microbes (Fuerst, 1995), was frequently found when the backward trajectories followed oceanic routes (Figures 3A, 5). PERMANOVA analysis showed a significant relationship between the routes and the estimated composition of ordinary microbial habitats (p < 0.05). Notably, the ratios of marine-related environments dominated when the air masses originated from the oceanic route, and animal-related environments dominated when they originated from the terrestrial route. Shannon's diversity indices of microbes became larger when the air masses originated from the terrestrial route (Shannon's diversity indices were 3.74 ± 0.68, 3.05 ± 1.00, and 3.15 ± 1.36 for the terrestrial, oceanic, and hybrid routes, respectively. The index of each sample is shown in Table 1); however, it should be noted that some samples required additional sequences to reach plateaus of rarefaction curves as mentioned already.

Soil, oceanic, and animal-associated microbes are spread in the atmosphere and transported for long distances (Morris et al., 2014; Smets et al., 2016), and precipitation may facilitate this microbial cycle. Sea-living microbes are emitted into the atmosphere via the bursting of bubbles on waves (Fahlgren et al., 2015), whereas soil-living and animal-associated microbes are transported on soil dust (Echigo et al., 2005; Prospero et al., 2005; Maki et al., 2011; Yamaguchi et al., 2014). In high-altitude atmospheric environments, microbes may be under substantial selection pressure due to harsh chemical, physical, and nutrient conditions (Delort et al., 2010; Morris et al., 2013; Smith, 2013). INA microbes play roles in cloud formation (Morris et al., 2013) and may facilitate the return of aerial microbes to diverse environments. The dispersal of pathogenic microbes causes disease epidemics that threaten public health and agricultural plant and animal health (Brown and Hovmøller, 2002; Rodó et al., 2011, 2014; Cao et al., 2014). Continuous long-term monitoring and large-scale analysis of precipitation microbes is thus envisioned to reveal the full impact of atmospheric microbial transport on microbial ecology, microbial evolution, public health, and climate.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.