The airborne mycobiome and associations with mycotoxins and inflammatory markers in the Norwegian grain industry

Grain dust exposure is associated with respiratory symptoms among grain industry workers. However, the fungal assemblage that contribute to airborne grain dust has been poorly studied. We characterized the airborne fungal diversity at industrial grain- and animal feed mills, and identified differences in diversity, taxonomic compositions and community structural patterns between seasons and climatic zones. The fungal communities displayed strong variation between seasons and climatic zones, with 46% and 21% of OTUs shared between different seasons and climatic zones, respectively. The highest species richness was observed in the humid continental climate of the southeastern Norway, followed by the continental subarctic climate of the eastern inland with dryer, short summers and snowy winters, and the central coastal Norway with short growth season and lower temperature. The richness did not vary between seasons. The fungal diversity correlated with some specific mycotoxins in settled dust and with fibrinogen in the blood of exposed workers, but not with the personal exposure measurements of dust, glucans or spore counts. The study contributes to a better understanding of fungal exposures in the grain and animal feed industry. The differences in diversity suggest that the potential health effects of fungal inhalation may also be different.

Seasonal and climatic effect on fungal diversity patterns. The average OTU richness per sample was 239 ± 73 (average ± SD). The accumulation curve showed that the total number of OTUs detected from winter samples were relatively higher compared to autumn (Fig. 1), but the difference was non-significant in the ANOVA analysis (Table 1; Fig. S2). The sample with the highest species richness was from grain mills located in the continental subarctic climate of the eastern inland of Norway, with short summers and snowy winters (hereafter called climatic zone two or CZ2), and lowest in samples from grain mills located in the continental subarctic climate of central coastal Norway, with short growth season and lower temperature (hereafter called climatic zone one or CZ1) (Fig. 1). However, on average, significantly higher richness (ANOVA, p = 0.014) was observed in the samples from grain mills located in the humid continental climate of the southeastern Norway (hereafter called climatic zone three or CZ3) (281 ± 58) followed by CZ2 (251 ± 78) and CZ1 (171 ± 32) ( Table 1; Fig. S3). In contrast, total abundances of the top 10 most common OTUs varied significantly between both seasons (being high during autumn) and climatic zone (being more abundant in zone one followed by two and three) ( Table 1; Figs. S2, S3). Shannon diversity and evenness index values were significantly higher in winter compared to autumn, but differences among climatic zones were not significant for both indexes (Table 1; Figs. S2, S3).
Seasonal and climatic effect on fungal community structural patterns. The variation partitioning analysis revealed that different climatic zones accounted for highest amount of the variation (23%) in fungal www.nature.com/scientificreports/ community composition, whereas season contributed 4% to the total variation (Fig. S4). Variations explained by season and climatic zones were not shared and 74% remained unexplained. In line with the above results, we also found that 35% of the OTUs were unique to winter season and 19% were unique to autumn, while 47% were shared between both seasons (Fig. 2a). Among the climatic zones, 21% OTUs were shared, while the highest proportion of unique OTUs were detected in zone two (34%) and lowest in zone one (4%). We also observed that sharing of OTUs between climatic zone one and three was lowest, whereas between zone two and three was highest (Fig. 2b). Multivariate PERMANOVA and NMDS analysis resulted in the very similar structural patterns, and confirmed strongest impact of climatic zone on fungal community (R 2 = 0.35; p < 0.001; Table 2). Samples from different climatic zones were well separated from each other and those from the same climatic zone were grouped together in ordination space (Fig. 3a). Marginally significant seasonal effect on fungal com- Table 1. Effect of different seasons (autumn and winter) and three climatic zones on airborne fungal richness, Shannon diversity index, evenness and abundance of the 10 most common operational taxonomic OTUs (OTUs) in industrial grain mills. Statistically significant differences among studies factor variables were analyzed using ANOVA and Tukey's HSD post-hoc test. Bold values represent differences significant at a P-value < 0.05. Seasonal and climatic effect on fungal composition. The shift in fungal composition at different taxonomic levels (from phylum to species level) also existed among both season and climatic zones. In autumn and winter seasons, the Ascomycota dominated in terms of read abundances ( Fig. 4a; Table S1), but occurred in higher abundance during winter (62.7%), compared to autumn (55.1%) season ( Fig. 4b; Table S1). The fungal composition at the order taxonomic level was highly season specific, with Capnodiales and Microascales (both ascomycetes) being most abundant during autumn whereas Pleosporales, Hypocreales Helotiales (all ascomycetes) and Tremellales (basidiomycetes) were more common in winter ( Fig. 4c; Table S1). Order level composition in terms of occurrences was very similar between season ( Fig. 4d; Table S1). Fungal genera also showed differential abundance distribution between seasons, where both Bipolaris and Botrytis were more abundant during winter compared to summer, and significant differences were recorded for Sclerotinia (F-value = 5.07; p = 0.035; Fig. 5a). The most common fungal OTUs also responded to the seasonal variations. In NMDS ordination space OTUs matched to species Cladosporium herbarum (OTU 1), Cladosporium exasperatum (OTU 9004), Microascus sp. (OTU 8) and Fusarium spp. (OTU 2) showed their species optima towards the samples from autumn whereas OTU 30 matching to species Gibberella tricincta was more abundant during winter (Table 3; Fig. 3b). Stationary air samples from companies in CZ2 and CZ3 were dominated by the fungal phyla Ascomycota both in terms of abundance and OTU frequency, whereas in CZ1 basidiomycetes were more common (Table 3; Fig. 4a,b). At an order level, fungal composition among climatic zones varied largely in terms of their abundances but not by their occurrence (Table 3; Fig. 4c,d). The Capnodiales, Pleosporales, and Hypocreales were more abundant in CZ3 followed by CZ2 and CZ1, whereas exactly the opposite pattern was observed for the Tremellales and Sporiodiobolales. Climatic zone also affected the distribution of different fungal genera (Fig. 5b). Abundances of Microdochium (F-value = 4.91; p = 0.019), Gibberella (F-value = 3.93; p = 0.037) and Bipolaris (F-value = 5.44; p = 0.013) were significantly higher in CZ2 and lower in CZ1; whereas Sporobolomyces (F-value = 3.84; p = 0.039) and Botrytis (F-value = 4.95; p = 0.018) were more common in CZ2 and less common in CZ3. Sequences placed in the genus Itersonilia were significantly more abundant in CZ1 as well as CZ2, whereas the pattern was opposite in CZ3 (F-value = 5.91; p = 0.037). OTUs belonging to Sclerotinia were only detected from CZ2 (F-value = 3.92; p = 0.037). Only one Ascomyceta OTU (identified as Gibberella tricincta) and the majority of the most common Basidiomyceta OTUs (matched to Sporobolomyces ruberrimus, Vishniacozyma victoriae, Dioszegia crocea, Itersonilia pannonica and Cryptococcus sp.) showed strong preferences in CZ1 and low preferences in CZ2 and CZ3. Abundance patterns was exactly opposite for most of the common Ascomyceta OTUs (matching to species Cladosporium herbarum, Cladosporium exasperatum, Gibberella tricincta and Alternaria metachromatica) restricted to CZ3.
Correlation analyses between fungal compositional structure and occupational exposure and inflammatory markers. We investigated the relationship between fungal community structure and previously reported mycotoxins in settled dust (Table S2) 2 , as well as personal exposure to dust, fungal spores and 1,3-β-glucans (Table S3) 11 , and biomarkers of exposure in the blood of workers (Table S4) 27 . Using both procrustes (Correlation in a symmetric rotation = 0.37; p = 0.186) and mantel (r = 0.19; p = 0.137) analyses, we found that the overall mycotoxin level in grain dust was not significantly correlated with fungal community composition. However, individual mycotoxin analysis using envfit function with NMDS, showed that some mycotoxins were significantly correlated with the fungal community composition ( Fig. 3c; Table S5). These mycotoxins were variants of enniatin (A1, B, B1 and B3), asterric acid, and rubellin D. Personal exposure levels of dust, fungal spores and 1,3-β-glucans were not significantly correlated with fungal composition (Table S5). Among all tested inflammatory markers in the workers blood (CC-16, SP-D, sP-selectin, IL-6, TNF-α, fibrinogen, sCD40L and CRP), only fibrinogen (R 2 = 0.35; p = 0.04) was significantly correlated with the fungal compositional structure ( Fig. 3c; Table S5).

Discussion
In this study, the diversity and composition of fungal bioaerosols in stationary air samples derived from Norwegian grain elevators and animal feed mills was characterized. A broad spectrum of fungi was identified and the taxonomic groups characteristic of the grain industry in general belonged to the Ascomycota orders Capnodiales, Xylariales, Hypocreales, Pleosporeales and Eurotiales, as well as the Basidiomycetes orders Sporidiobolales, Tremellales, and Cystofilobasidiales. The airborne mycobiome identified in this study differed in fungal diversity and Table 2. Relative importance of different seasons (autumn and winter) and three climatic zones on the airborne fungal community composition in industrial grain mills, as revealed from the PERMANOVA analysis. Bold values represent differences significant at a P-value < 0.05.    www.nature.com/scientificreports/ composition depending on season and geographically and climatically different zones. The fungal composition was related to some mycotoxins in certain zones, and to fibrinogen from the blood of exposed workers, but was not correlated with personal exposure measurements of dust, glucans or spore counts. Differences related to climatic zones accounted for most of the variation in airborne fungal community structure, explaining 23% of the variance, and shown by the distinct fungal community structure and composition of each of the zones. The zones represent differences in geographic location and climatic conditions that likely influence the occurrence and fungal growth in the surrounding environment, and subsequently in the locally produced grain that enters the grain elevators, of which the airborne dust is generated. Relatively homogenous compositions of fungal communities were found in wheat grain dust from 96 fields across 560 km 2 in a Swiss study 21 . This shows that the geographic distance may not necessary be the strongest determinant for fungal diversity differences in general. However, an important role of geographical distances in structuring the fungal communities of other crop types, such as wineyards, have been indicated 29 . Whereas Pellissier and colleagues analysed the mycobiome of wheat dust, the present study included dust from wheat, rhy, oat and barley of various cultivars, which are likely to influence the diversity pattern differently 21,30 . Although all sampling zones are located north of 59 degrees latitude, and therefore have harsh winters, the climate is mildest in CZ3; the Southeast of Norway, followed by CZ2; Eastern inland of Norway and CZ1; the central coastal Norway (Fig. 6). It is therefore plausible that the airborne grain dust mycobiome of grain mills of the Southeast of Norway had higher fungal diversity in terms of species richness and Shannon diversity index than those of the Eastern inland of Norway or the central coastal Norway. CZ2 and CZ3 are geographically closer to each other and shared 19% of OTUs, whereas CZ1 shared only 7% OTUs with CZ2 and 2% with CZ3. CZ1 was also dominated by sequences placed in the phylum Basidiomycota, compared to the other districts that were dominated by the Ascomycota. This is also in contrast to Ascomycota species that are considered more robust, as CZ1 is considered a harsher environment in the northernmost district in the study. However, although short grain growth period, the coastal climate around the Trondheim fjord characterized by a mild climate, and prime agricultural district, the CZ1 may give better growth conditions for fungi placed in the Basidiomycota.
The distinct fungal compositional profiles between each climatic zone indicates that employees in different zones could be differently exposed. Fungal exposure assessments representative of all zones should therefore include measurements among the 21% of OTUs shared among the zones. However, it will be important to investigate what the fungal profile difference may imply for occupational health of the workers in the different zones.
We observed higher fungal diversity and evenness index values, and a tendency to higher richness, in winter compared to autumn. This result may be somewhat surprising, considering that autumn may represent more favorable growth conditions and higher rate of spore dispersal due to higher temperature and nutrient sources. It is also in contrast to other studies that have shown seasonal differences in fungal communities associated with sawmills 31 , plant leaves 32 , roots 33 and soil 34 . However, the grain dust generated in winter results from handling Table 3. Taxonomic affinity and occurrences of the 20 most abundant operational taxonomic units (OTUs) detected industrial grain mill air during different seasons (autumn and winter) and in three different climatic zones. For overall communities both occurrences and % total reads are shown, but for season and climatic zone only % total reads data is provided here. Ide % identification against reference sequence from UNITE database, Cov % sequence coverage. www.nature.com/scientificreports/ more heterogeneous grain types received during winter than grain types that are received during autumn. As the grain at this time of year has already dried, and therefore will not have the need to go through drying processes at the elevators, the grain turnover at the mills may be higher, leading to more heterogeneous fungal dust. The grain also likely originates from several diverse Norwegian farms that has dried the grain themselves. Further, as the need for grain is larger than the local production, a considerable amount of grain imported from abroad is additionally received and handled, generating dust containing fungi from broad variety of geographic environments. Grain handlers would thus be exposed to different fungal species in autumn and winter, although the seasonal difference in species composition accounted only 4% of the variation. A greater relative abundance of sequences placed in the phylum Ascomycota than Basidiomycota was detected in winter and autumn samples, which is in line with the view of Ascomycota species are more tolerant to environmental stress. In this case, this would not mean the harsh winter conditions, as previously indicated for the mycobiome in Norwegian sawmills 31 , but more likely the production process-related strain that fungi are exposed to during drying and moving the grain at the elevators and feed mills. However, we observed seasonal differences at the order and genus levels that may be explained by different nutrient preferences or tolerances of the identified fungi. The abundance of the 10 most common OTUs were higher in autumn than in winter. This relationship might be related to the evenness of the samples, i.e. if there are fewer species during winter, these will appear as more abundant relatively. Hence, the winter season at grain elevators and animal feed mills most likely represent the most diverse fungal airborne exposure potential.
The stability of mycotoxins make them remain in the dust although their fungal producers die and disintegrate. Species not producing mycotoxins may also be present. This can result in low correlations observed between mycotoxins and potential producers detected in the dust. Despite this, we observed a significant correlation between the fungal community structure and enniatins, asterric acid and rubellin D. Enniatins were www.nature.com/scientificreports/ present in high concentrations in the dust, whereas asterric acid and rubellin D were present in lower levels (Table S3). Rubellin D is a phytotoxin produced by Mycosphaerella rubella 35 and Ramularia collo-cygni 36 , and has also antimicrobial, antiproliferal and cytotoxic properties. These fungi belong to the Capnodiales, the most abundant order of the study, and Mycosphaerella sp. occurred in 45% of the samples whereas Ramularia were found in 82% of the samples, thereof 64% Ramularia collo-cygni (32,508 total reads and 1478 average reads). Enniatins are depsipeptides that are naturally produced by the filamentous fungi in the genera Fusarium 37 , Halosarpheia and Verticillium 38 . The main producer, Fusarium, occurred in 91% of the samples (560,726 total reads and 6372 average reads), Verticillium only in 9% of the samples (n = 2; 42 total reads), whereas no Halosarpheia were observed. Enniatins have a range of biological activities, including anti-insectian, antifungal, antibiotic, and cytotoxic 39 . Asterric acid is an inhibitor of vascular endothelial growth factor (VEGF) 40,41 , and is produced by Aspergillus and Penicillium species, which were found in 82% and 55% of the samples, respectively. Although the mean fungal spore exposure and 40% of the samples were above the proposed occupational exposure limit for fungal spores of 10 5 spores m −3 (Table S4), the fungal species composition did not correlate with the spore counts, neither 1, 3-β-glucans nor dust exposure. As PCR-amplified and MiSeq-sequenced fungi will include hypha and fragments in addition to spores, a correlation with SEM counts of spores was neither expected. It is worth noting that results from metabarcoding represent a completely different picture than the traditional occupational exposure measurements, and the details obtained should be taken into consideration for refining the traditional methods into more targeted measurements. The cause-and-effect relationship between complex microbial occupational exposures and health effects is not clear, but the characterization of the microbial complexity of occupational exposures is a step on the way. The investigation of the relationship between the exposure diversity and inflammatory and anti-inflammatory exposure responses and potential health effects of complex exposure is an important next step towards increased understanding. Barrera and colleagues did indeed find differences in the exposure responses between grain handlers and cattle raisers that seemed to be related to differences in the airborne mycobiome 25 .
Several predominant fungal OTUs identified in this study belong to the genera Aspergillus, Penicillium, Alternaria, Cladosporium and Wallemia, which include plant and animal pathogenic and allergenic species that may cause respiratory effects [42][43][44] . Other common crop-plant associated and plant pathogenic fungi, such as species of the genus Fusarium, can have large impact on human health through their immunosuppressive, toxigenic, teratogenic, and carcinogenic mycotoxins 45 . OTUs from the Basidiomycota genus group Cryptococcus was found in all samples, but only 9% matched to the animal-pathogenic infectious species, Cryptococcus neoformans, that can cause cryptococcosis, a systemic fungal disease that affects internal organs and skin of immune deficient individuals 46 . The species Vishniacozyma victoriae, which has been suggested to be inversely associated with asthma 47 , was detected in the majority of the samples.
Although blood concentrations of both CC-16, IL-6 and fibrinogen were significantly different between exposed and controls (Table S5) 27 , none of the markers were associated with dust exposure or any other exposure component (fungal spores, endotoxins or 1,3-β-glucans) measured during work 27 . This was in spite of exposures to elevated bioaerosol concentrations 11 . However, when compared with the fungal species composition characterized in the present study, the fibrinogen concentration was associated with fungal diversity at company level. Apart from playing a role in the acute phase response to systemic inflammation, fibrinogen has been shown to be hydrolyzed by proteinases from Aspergillus melleus and other allergenic organisms to yield cleavage products that drive distinct antifungal mechanisms and initiates allergy (airway hypersensitivity, innate allergic inflammation and antifungal immunity) through Toll-like receptor 4 and CD11c 48 . Indeed, exposed workers had lower levels of fibrinogen compared with controls 27 , but although this may be a possible ongoing mechanism in the exposed workers, most employees in both study groups were within the reference area for fibrinogen of 1.5-4 mg mL −1 . The restriction of the larger individual dataset of the inflammatory markers to the pooled dataset at company level represents a necessary limitation for the NMDS/PERMANOVA analyses. The results as such cannot be used to identify risk factors or quantify their significance. This has to be taken into consideration when interpreting the results. Further analyses beyond this point, such as the identification of specific fungal species or species profiles that represent a particular exposure-response should be conducted with quantitative epidemiological methods in the future.
In conclusion, this DNA metabarcoding study of the fungal ITS2 region in airborne grain dust provides insight into the complete spectrum of fungal bioaerosols across seasons and climatic zones in the Norwegian grain industry. This contributes to a better understanding of fungal exposures in the occupational workforce of the grain and animal feed industry. The differences in diversity suggest that the potential health effects of fungal inhalation may also be different. Fungal diversity differences should therefore be included in future risk assessments. This inclusion can lead to the development of improved exposure assessment methodologies to clarify exposure-response relationships.

Materials and methods
Sampling and sampling sites. Sixty-eight grain dust exposed workers at twenty companies of industrial grain elevators and compound feed mills participated in this study through the winter and autumn season of 2008. Four sets of samples were collected (Table 4): (1) Personal samples of the workers exposure to inhalable dust and bioaerosol components during a full work shift. The workers carried the sampling equipment during work, and the collected dust represented the workers inhalable exposure during all tasks performed that day; (2) Stationary full-shift sampling of airborne dust subjected to mycobiome analyses. Samplers were placed 1.5 m above ground at a fixed position assumed to be representative for the activity at each company; (3) Freshly settled grain dust from the production facilities subjected to mycotoxin analyses; (4) Blood samples from exposed www.nature.com/scientificreports/ workers. Blood were collected after work from the same persons that carried the equipment for exposure measurements.
The sampling spun over autumn season and winter season, and three geographically and climatically different grain districts; (1) central coastal Norway, characterized by continental subarctic climate with short growth season, rain, low temperature and 850 day-degrees, (2) eastern inland of Norway, also with a continental subarctic climate, but with less rain, but more snow and about 1000 day-degrees, and (3) south-eastern Norway, with humid continental climate with humid and warm summers, and 1200 day-degrees (Fig. 6). The districts are referred to as climatic zone 1, 2 and 3, respectively, in the text. All districts are located north of 59 degrees northern latitude. The industrial process of the study plants has been described previously 19 . In grain elevators, various grain was loaded into the elevator, sorted, winnowed, dried, rotated, moved, stored and unloaded on a continuous basis. Grain moving was either air-driven (suck-and-blow), elevator-driven, or done by passive emptying by gravity. In compound feed mills, the grain was milled and mixed with other nutrients such as maize, calcium carbonate, fat, vitamins, and amino acids, pressed into feed pellets and filled into sacks or tanks. In autumn, a large part of the grain loaded into the grain elevator came from local producers and could include humid batches that needed drying. In winter season, the delivered grain had been dried on the farms before delivery, and the amount of imported grain was higher than in autumn.
Stationary sampling of airborne dust for mycobiome analyses. Stationary sampling of airborne dust (n = 25 samples in total) was performed in 20 grain elevators and compound feed mills in Norway using NIOSH 1 stage samplers (NIOSH, Morgantown, USA) connected to SKC high flow portable battery-powered pumps (SKC Inc., Valley View, PA, USA) with an airflow of 4 L min −1 . The samplers included a 1.5 mL Eppendorf tubes, collecting particles above 1.5 µm in aerodynamic diameter (theoretical particle size cut-off), and a 25 mm polycarbonate filter with pore size 0.8 µm placed in a 25 mm black IOM cassette (SKC Inc.), as a back filter, collecting particles smaller than 1.5 µm. The time of sampling was 235 to 475 min (mean 368 min). The samplers were placed in 1-2 m height at work sites representative for the activity of the companies. DNA extraction, amplification and sequencing. DNA was extracted from dust collected in the Eppendorf tube part of the sampling equipment. 150 µL cetyl trimethylammonium bromide (CTAB) buffer were initially added, and the DNA was extracted as previously reported 49 . Fungal ITS2 rRNA gene region was amplified using the fITS7a and ITS4 50,51 primer combination. We followed similar PCR reactions setup and conditions as well as amplicon library preparation method as mentioned previously 49 . In addition to the 24 samples, individual cultured species, 4 parallels of a mock community, 2 field blanks, and a PCR negative control were also included in the sequencing run as control samples, and were used while curating the final OTU tables as shown before 31 . The mock community comprised fungal cultures of Rhizopus microsporous, Cladosporium allicinum, Penicillium variotii, Fusarium graminearum, Aspergillus fumigatus and Wallemia sebi, with known genome size, and was used to estimate index jumping in the sequencing run. Molecular data were generated using Illumina Miseq Paired-End (PE: 2 × 300 bp) sequencing.

Mycotoxin analyses.
Multiple mycotoxins and other metabolites of fungi were extracted from freshly settled grain dust (1.5-15 g) collected from the 20 grain elevators and animal feed mills, and analysed with LC-MS/ MS. The analyses and results have been described previously 2 , and in the present study, the arithmetic mean concentrations from each company (Table S3) was used for correlations analyses with fungal compositional data.
Aerosol exposure measurements and analyses. Fifty-six full-shift personal inhalable samples were collected with PAS-6 personal inhalable samplers 52 and portable pumps (PS101; National Institute of Occupational Health, Oslo, Norway) for 6-8 h with a flow rate of 2 L/min. Samples were analyzed for grain dust, fungal spores and β-1,3-glucans as previously described 11 . In brief, dust was weighed, fungal spores were counted by scanning electron microscopy, and β-1,3-glucans were analysed by enzyme immunoassay. The results of the per- Table 4. Overview of the acquisition and number of samples for each sample set. a For bioaerosols and biomarkers, the number of samples equals the number of exposed workers. b The number of bioaerosol samples represent parallel sampling for analyses of dust, fungal spores and β1 → 3-glucans. CZ climatic zone. The number of companies grouped in CZ1, CZ 2 and CZ3 was 6, 7 and 7, respectively. All samples were collected at the same day during winter and autumn 2008. www.nature.com/scientificreports/ sonal exposure have been published previously 11 . In the present study, the mean concentration of each aerosol component from each company (Table S4) was used for correlation analyses with fungal compositional data.
Blood sampling and analyses of biomarkers. Blood samples were collected after work between 1 and 3PM. Club Cell protein 16 (CC-16), surfactant protein D (SP-D), SP-A, interleukin 6 (IL-6), serum P-selectin (sP-selectin), tumor necrosis factor-alpha (TNF-α), fibrinogen, C-reactive protein (CRP) and serum CD40 ligand (sCD40L) were analyzed by ELISA and reported as previously described 27 . In the present study, the mean concentration of each biomarker measured in the blood of workers at each company (Table S5) was used for correlation analyses with fungal compositional data.
Bioinformatics processing. The Miseq sequence data were processed through the same bioinformatics pipeline as previously reported 49 . Taxonomic assignment were made by comparing the representative sequence against curated reference database UNITE v6 53 . The data treatment strategy generated an OTU table where all the control samples (extraction, negative and field blanks) were perfectly clean with no reads. Furthermore, mock composition represented by mock communities was as per our expectations. This indicated the appropriateness of the data treatment procedure to account for index hopping. Of total 25 samples, one sample failed during sequencing run, leaving 24 samples in the final dataset. The dataset revealed non-normal spreading of reads among the samples (Fig. S5a) and distribution patterns of reads per OTU were also skewed, with only a few very abundant OTUs and many with very low abundances (Fig. S5b). Per sample rarefaction curves of OTU richness indicated that complete diversity was not captured for majority of the samples (Fig. S6). Additionally, a positive relationship was also observed between OTU richness and sequencing depth (R 2 = 0.34; p = 0.004). To overcome this possible sequencing depth bias, we decided to normalize (43,837 reads per sample) the dataset for richness and diversity analysis, and were subsampled randomly. To account for inconsistencies in abundance measures and further to improve variance homogeneity dataset was arcsine-transformed prior to analyses. Two samples that had extremely low reads (35 and 2314 reads respectively) were removed from both community and diversity analysis, and 22 samples remained in final dataset, used for overall community analysis. Due to lack of mycotoxin data from four samples, only 18 samples were used while analyzing fungal and mycotoxin dataset. To achieve homogeneity of variances the values of the studied variables (mycotoxins in settled dust, airborne exposure levels of dust, fungal spores, and 1,3-beta glucan, and blood concentrations of CC-16, SP-D, SP-A, sP-selectin, IL-6, TNF-α, Fibrinogen, sCD40L and CRP) were transformed to zero skewness and expressed on a 0-1 scale prior to analyses. Statistical analysis. Unless stated otherwise, statistical analyses were performed in R v3.5.0 54 . Richness, Shannon and evenness diversity measures and total abundance the most common top 10 OTUs were calculated using R package vegan, and their differences between two seasons (autumn and winter) and three climatic zones were examined using ANOVA followed by a Tukey's HSD post-hoc test from R package agricolae and visualized using violins plots. Seasonal and climatic zone differences in proportional abundances of different genera were tested using ANOVA (Benjamini-Hochberg FDR correction was applied) followed by a Tukey's HSD post-hoc test and illustrated by heat plots based on hierarchical clustering using R package vegan. Species accumulation curves for number of accumulative OTUs for different seasons and climatic zones were calculated using function specaccum from R package vegan with 10,000 permutations. Overall taxonomic composition of airborne fungi (up to family level) was visualized with heat tree generated using R package Metacoder 55 .
The Bray-Curtis dissimilarity index was used to generate community distance matrices and used further in all community structure analysis. To address the relative importance of seasons and climatic zones on fungal community structure, multivariate permutational analysis of variance (PERMANOVA) was used, as implemented in the Adonis function of the package vegan. PERMANOVA analysis was executed using a forward selection practice to improve the final model 56 . First, we examined single variable models and thereafter, including significant factors in the final model in order of their R 2 values. Furthermore, Nonmetric Multidimensional Scaling (NMDS) ordinations analyses were used to visualize the effects of studied factor variables on fungal community composition using the metaMDS function of the package vegan 57 . Centroids of the factors were fitted into NMDS plots using the function envfit. Ordiellipse function were used to plot the 95% confidence intervals (CI) of the different seasons.
To assess correspondence between overall mycotoxin and fungal community dataset, we used both mantel and procrustes analysis with 999 permutation, where distances matrixes of both datasets are compared together. In order to explore fungal community structural relationship between different mycotoxins in settled dust, airborne exposure levels of endotoxin, bacteria, fungal spores, and 1,3-beta glucans, blood concentrations of CC-16, SP-D, P-selectin, IL-6, TNFα, fibrinogen, CD40L and CRP, we further used NMDS ordination with the envfit function. In this test each explanatory variable is separately regressed on ordination axes 1 and 2 by linear regression analysis. The strength of the relationship between each variable and the ordination axes was assessed by the multiple coefficient of determination (R 2 ) of this regression. A significance level of p < 0.05 was used throughout the study.
To further partition the variation of community dissimilarity from different seasons and climatic zone, we followed redundancy analysis (RDA) from the R package vegan (varpart function). The R 2 values were adjusted according to each predictor (factor) level.
Ethic statement. The Regional Ethical Committee of South-East Norway and the Norwegian Data Inspectorate approved the study. All participants gave their written informed consent upon participation in the study. All but one of the workers that received the written information agreed to participate with blood samples. The www.nature.com/scientificreports/ one that refused blood sampling gave written informed consent to participate in the exposure measurements only. All experiments were performed in accordance with relevant guidelines and regulations. The study was approved by the Regional Ethical Committee of South-East Norway, and received support from the Confederation of Norwegian Enterprise (S-2585).