Fungal microbiome of barley grain revealed by NGS and mycological analysis

Introduction. Barley can be infected with a broad variety of fungi, which can cause considerable loss of crop yield and reduce the quality of grain. Modern vision on the geographical and ecological distribution and biodiversity of micromycetes has been established by traditional, cultivation-based methods. However, more recently, molecular methods have shifted microbiological research to a new level, making it possible to investigate hidden taxonomical biodiversity. Study objects and methods. For this study, we determined the fungal biome on the surface and inside of barley grains using the traditional mycological method and the contemporary molecular method, which employed DNA metabarcoding based on NGS (nextgeneration sequencing) of the ITS2 region. We analyzed five cultivars that were collected in two subsequent crop seasons (2014, 2015). Results and discussion. DNA metabarcoding revealed 43 operational taxonomic units, while 17 taxa of genus or species level were recovered by the traditional method. DNA metabarcoding revealed several minor species and one predominant, presumably plantpathogenic Phaeosphaeria sp., which were not detected in the agar plate-based assay. Traditionally, Fusarium fungi were identified by mycological assay. However, the resolution of DNA metabarcoding was sufficient to determine main Fusarium groups divided by ability to produce toxic secondary metabolites. The combined list of Ascomycetes consisted of 15 genera, including 14 fungi identified to species level. The list of Basidiomycota derived from DNA metabarcoding data alone included 8 genera. Conclusion. It was found that crop season predetermines the fungal community structure; mycobiota on the surface and inside of grain was significantly different.


INTRODUCTION
Barley (Hordeum vulgare L.) is one of the major cereal crops. It occupies fourth place among cereals in the world and second place in Russia by production quantity and cultivation area [1]. The importance of barley has been accepted since ancient times and used in the food, feed and brewing industries due to its versatility, excellent adaptation capabilities and superior properties [2].
The increased interest in barley as a source of food and fodder has resulted in a huge number of studies of associated microorganisms. It is known that barley can be infected with a broad variety of plant-pathogenic and toxigenic fungi, many of which may persist in grains.
The Bipolaris, Pyrenophora, Phaeosphaeria, Alternaria, and Fusarium genera are considered to be prevailing fungi in barley grain worldwide [3,4]. Species of the last two genera are well known as mycotoxin producers, with Fusarium spp. being the most dangerous food and feed contaminants.
Cultivation-based methods have traditionally established modern vision on geographical and ecological distribution and biodiversity of micromycetes. These methods cannot provide accurate data on taxon composition because some microorganisms do not have specific characteristics to be identified, and some appear to be noncultivable. Thus, data on the mycobiome of many substrates, including barley grain, is likely to be incomplete.
In recent decades, molecular methods have shifted microbiological research to a new level, making it possible to investigate hidden taxonomical biodiversity. Next-generation sequencing (NGS), implemented on various independent technical platforms, became the most promising method for conducting research projects aimed at revealing fungal or bacterial composition [5][6][7]. Several studies have focused on a variety of agricultural subjects [8][9][10][11][12]. Advances in this field led to consideration that NGS-based methods are suitable as incipient techniques for seed testing [13].
In Denmark, 454 pyrosequencing of the internal transcribed spacer 1 (ITS1) of the nuclear ribosomal DNA (rDNA) has been chosen to recover the composition of fungal communities associated with wheat grain [14]. NGS revealed a significantly higher level of biodiversity than it was observed in previous culturing studies. Another appropriate 454 pyrosequencing of both ITS regions was done to study the mycobiome of barley grain in western Canada [3]. It demonstrated that geographic location and agronomical practices were the determining factors explaining the observed differences in the fungal communities associated with barley. Such studies may contribute to a better understanding of fungal species compositions in cereals. They may also lead to more accurate food-quality testing and the precise design of crop protection strategies that would reduce the level of fungal contamination of agricultural products.
The objective of this study was to revise the taxonomical variety of fungi contaminated the surface and infected barley grains harvested in the northwestern region of Russia. We hypothesized that grain mycobiome could be significantly differ on the surface and inside the grain, and the difference may depend on the crop year. In our research we used the traditional agar platebased method and the contemporary method based on 454 pyrosequencing of ITS2.

STUDY OBJECTS AND METHODS
Sampling. Grain samples of five spring barley cultivars (Suzdalets, Krinichnyy, Moskovskiy 86, Tatum, and Belgorodskiy 110) were received in 2014 and 2015 from the State Experimental Station (Volosovo, Leningradskaya oblast, Russia, 59°31'N, 29°28'E). Small grain cereals on this station were cultivated with no fungicide treatments. The grain samples intended for fungi isolation and DNA extraction were stored separately at 4°C and -20°C respectively.
Extraction of DNA. Two representative subsamples consisting of 500 grains were picked from each sample. One subsample was placed in a 50 mL polypropylene tube and subjected to superficial sterilization. The grains were consistently washed with 20 mL of 2% sodium hypochlorite solution containing 0.1% of sodium dodecyl sulphate (SDS). They were then washed once with 5% sodium hypochlorite, two times with deionized water (diH 2 O), and finally rinsed with 98% ethanol. With each step the mixtures were actively stirred up within 2 min, and the flushing solution was then decanted without the grain. The ethanol was removed by burning at regular stirring during 10 s. After this, the grains were homogenized in sterile disposable chambers on a Tube Mill Control (IKA, Germany) grinder. The other subsample was similarly homogenized, but the step of superficial sterilization was skipped.
Further, 240 mg of the ground subsamples were transferred to a 2 mL Eppendorf tube, where DNA was extracted with an AxyPrep Multisource Genomic DNA Miniprep Kit (Axygen, USA), according to the centrifugal protocol for plant tissues and fungal mycelium. The DNA concentrations were measured with a Qubit 2.0 (Thermo Fisher Scientific, USA) using a dsDNA HS Assay Kit. Extracted DNA was used for library preparation and subsequent universal tailed amplicon sequencing, as described for the 454 Sequencing System.
Pyrosequencing and primary data analyses. For amplicon library preparation we chose the taxonomically significant ITS2 region, which is commonly used, as well as ITS1, in DNA metabarcoding studies of fungal diversity. To a large extent, ITS1 and ITS2 have similar results when used as DNA metabarcodes for fungi [15][16][17]. However, the ITS2 region lacks the insertions commonly found in ITS1 and thus reduces length variation [18]. This is important, as length variation can bias community pyrosequencing toward shorter amplicons. Also, ITS2 is the best-represented fungal genomic element in the public databases [19,20]. Therefore, in studies similar to our project, use of ITS1 obtained with fungal specific primer (ITS1F) [21] can be helpful in eliminating plant ITS amplification, and may turn out to be the method of choice in cases of mixed plant and fungal genomic DNA. However, it is necessary to take into account that ITS1F (with constricted, specific range toward exclusion of all eukarya except fungal taxa), may not be able to amplify several fungal taxa because it is hampered with a high degree of mismatches relative to the target sequences [22,23].
The amplicon library pool was sequenced with 454 pyrosequencing on the GS Junior sequencer (Roche, USA) according to the recommendations of the manufacturer [25]. The ITS2 locus reads were processed by QIIME Version 1.6.0 (Quantitative Insights into Microbial Ecology) [26]. To reduce the amount of erroneous sequences and thus increase the accuracy of the whole pipeline, the denoising procedure was employed [27].
Next steps included assigning multiplexed reads to samples based on their specific MIDs (demultiplexing), removing the low-quality or ambiguous reads, truncating primers, and other accessorial sequences. Chimeric sequences were detected using the UCHIME algorithm with the Unite database [28][29][30]. All of the reads were clustered into operational taxonomic units (OTUs) at 97% sequence similarity using the UCLUST method [31]. Representative sequences were chosen according to their abundance between similar reads. Low-abundance OTUs, which have less than four copies (singletons, doubletons and tripletons), were deleted over all of the analysis [32]. All 454 pyrosequencing data of the present investigation are available through the Sequence Read Archive (SRA) under BioProject PRJNA353503, with run accession numbers from SRR5022991 to SRR5023010 [33].
Phylogenetic and statistical analyses. Taxonomical identification of representative sequences was carried out by the BLAST method using Genbank databases [34,35]. Query coverage ≥ 99% was recognized as significant. Query identity of ≥ 99% was considered identification at the species level; identity of ≥ 98-95% was considered reliable identification at the genus level. The smaller similarity Ribosomal Database Project classifier, along with the Unite database (minimum confidence at 0.9), were implemented to assign OTUs to a higher taxonomic rank [29,30,36,37].
Alignment of representative sequences was carried out using MAFFT algorithm G-INS-1 [38]. A phylogenetic tree was conducted with MEGA5 using the Maximum Likelihood method, based on the Tamura-Nei model with 1000 bootstrap replicates [39][40][41].
Vegdist and hclust R functions were used for computing Bray-Curtis dissimilarity indices and UPGMA hierarchical clustering of OTUs, showing their coexistence in the samples [42]. Heatmaps were generated with QIIME 1.8.0 with log-transformed abundance data. OTUs were sorted by phylogenetic or hierarchical trees. Beta diversity between samples was calculated by beta_diversity.py script in QIIME with unweighted UniFrac metric [26,43]. To check the robustness of estimated beta diversity, jackknifed analysis, with 96 reads per sample depth and 100 replicates, was performed. The results were visualized with Principal coordinate analysis (PCoA) in a 2D scale plot. Bray-Curtis, weighted and unweighted UniFrac dissimilarity indices were used for measuring the strength and significance of sample groupings with Permutational Multivariate Analysis of Variance (PERMANOVA) and Analysis of Similarity (ANOSIM) with script compare_categories.py [26].
Agar plate-based method of isolation and identification of fungi. Representative subsamples (200 grains) of each cultivar were surface-sterilized by being shaken for 2 min in 5% sodium hypochlorite. Then they were rinsed twice in sterilized water. The surfacesterilized grains were placed on 90 mm Petri dishes (10 grains per dish) with potato-sucrose agar (PSA) and incubated at 24°C for 10-14 days. The isolated fungal colonies from every grain were identified by visual and microscopic observations according to Ellis, Gerlach and Nirenberg, Lawrence, Rotondo, and Gannibal, and Samson et al. [44][45][46][47]. To present data comparable to those obtained with NGS, relative abundance was calculated as the number of all isolates of a certain taxon divided by the total number of fungal isolates (%). A more conventional index of seed expertise, infection frequency, was calculated as the number of grains infected by the fungus (%).
All of the OTUs were assigned to Basidiomycota and Ascomycota phyla (Fig. 1).
One of the most abundant OTU, otu166, assigned as Dothideomycetes, failed to be identified to a more precise taxonomical level. It has similar characteristics (BLASTn 100% query coverage and identity) with several GenBank sequences (e.g., EU552134, AJ279448, HG935454) which can be joined together only at the rank of class. More likely, it coincides with Epicoccum nigrum, the one abundant Dothideomycete identified during mycological analysis.
Eleven minor OTUs (Alternaria related otu44, otu53, otu61, otu180, and otu190; Cryptococcus related otu29, otu60, otu218, and otu264; Davidiella related otu123; and Fusarium related otu89) had no close relation to any known species but appeared in the same samples where a major OTU of a certain genus was abundant. However, potentially such satellite OTUs represent rare and/or poorly studied species, but most likely they are technical errors or sequence variances, which can occur despite all filtering and trimming procedures.
From ITS2 sequences combined into six OTUs and designated as Fusarium, several OTUs can be readily assigned as two synapomorphic (Fig. 4) clades, similar to that described by Watanabe et al. [49]. Two OTUs (F. poae -otu224 and Fusarium sp. -otu259) were abundant, but four OTUs appeared as solitary sequences.
Distribution of Fusarium-related OTUs among clusters corresponded to the prevalent toxic secondary metabolite production. The first cluster consisted of Fusarium fungi that are able to produce the trichothecene group metabolites. The subcluster 1a included F. sporotrichioides and F. langsethiae, which are the main producers of type A trichothecenes (like T-2 and HT-2 toxins); the subcluster 1b included F. culmorum and F. graminearum the producers of type B trichothecenes (like DON, or NIV). Species F. poae (subcluster 1c) produces trichothecenes of types A and B, and enniatins (ENNs). Fungus F. equiseti (subcluster 1d) is able to produce ENNs, but according to some authors, it can also produce a small amount of type A trichothecenes [50,51]. The subcluster 2 brought together Fusarium fungi that are able to synthesize ENNs: F. tricinctum, F. avenaceum, and F. lateritium [52].
Fourteen In general, the mycobiome of nonsterilized barley grains was characterized by a greater abundance of Basidiomycetes in comparison with surface-sterilized grains. The most abundant fungi in nonsterilized grains were Davidiella (Cladosporium) spp. and Cryptoccocus spp. After surface sterilization, the average abundance of Fusarium, Alternaria, Pyrenophora, and Phaeosphaeria, as well as fungi from Dothideomycetes, increased, but the ratio of those taxa depended on the year.
The fungal species composition of non-sterilized grains differed from the mycobiome of surfacesterilized grains primary due to a higher abundance of Basidiomycetes (Cryptococcus spp. and other Tremellales, and Cystofilobasidium macerans and other Cystofilobasidiaceae) and Davidiella (Cladosporium spp.) in the non-sterilized grains. All Basidiomycetes disappeared or became sparse after surface sterilization. The most abundant of them, Cryptococcus tephrensis (otu204) and C. victoriae (otu124), were also revealed inside grains but in fewer samples and in lesser amounts. Several OTUs, e.g., Cryptococcus wieringae (otu183), Mrakiella sp. (otu254), and Dioszegia sp. (otu303), tended to present on seed surfaces during only one year. Mycobiomes observed in two different growing seasons differed by abundance of Pyrenophora sp. in 2014 and Fusarium spp. and Phaeosphaeria sp. in 2015. The Alternaria, Bipolaris, and Epicoccum genera were relatively abundant in both sample sets. More detailed results of fungal coexistence in the samples are introduced in Fig. 6.
Agar plate-based method of identification of fungi. From 87 to 117 fungal isolates per 100 grains were obtained from each sample. As a result of two growing seasons, a total of 18 taxa of seed-borne fungi were identified (Table 1). The taxonomic position of some fungi was vague due to the lack of sporulation (Mycelia sterilia). In both years, the appearance of Alternaria, Bipolaris, Epicoccum, and Fusarium species  In both years, fungi of the genus Alternaria predominated in barley grain samples. The members of two sections, Alternaria and Infectoriae, were determined. More precise identification was not performed, since species concept is debatable for Alternaria and Infectoriae [53][54][55][56] sections.
Contamination of the barley grains by Fusarium spp. varied significantly in 2014 and 2015. In 2014, the Fusarium infection frequency was low (0-4%) and represented by five species, of which F. avenaceum was the most frequent (infection frequency up to 2.5%, relative abundance of isolates up to 2.2%). In 2015, the infection of barley grains with Fusarium spp. was considerably higher (infection frequency 14-19%, isolate abundance 12-17%). Eight Fusarium species were identified; four of them were common for both years.
Comparison of methods. In total, 43 OTUs assigned as Ascomycota (26) and Basidiomycota (17) were revealed by DNA metabarcoding. Only 14 OTUs were assigned to species level. From those species, only two were reoccurred in traditional mycological analysis. The other 12 species either were not detected among isolates grown up from grains on agar medium or were Mycelia sterilia. At the same time, the conventional mycological seed test revealed 17 Ascomycetes, including 11 species, apart from some Zygomycetes and sterile Ascomycetes. Basidiomycetes were not recovered by conventional assay. Two species (Fusarium poae and Bipolaris sorokiniana), one section (Alternaria sect. Alternaria), and three genera (Davidiella [Cladosporium], Fusarium, and Pyrenophora) were formally common for both assays. In general, the list of undoubtedly identified dominant taxa coincides with the results of the NGS mycobiome study of Canadian barley grains [3]. The predominant OTUs from Alternaria were identified as Alternaria and Pseudoalternaria sections when Alternaria and Infectoriae sections were fixed during mycological analysis. In both cases, taxa were identified to the section level. Such precision is sufficient for the majority of practical purposes, e.g., for tests of seed, food, or feed-grain quality. The big section Infectoriae and lately described section Pseudoalternaria are morphologically similar and phylogenetically close groups [46]. This obviously can be the cause of errors, if identification is based on morphological features.
Both methods similarly reflected a very low abundance of Fusarium spp. in 2014 and a higher quantity in 2015. Traditional mycological analysis revealed nine Fusarium species. DNA metabarcoding results were more limited; only one OTU was identified as a certain species, F. poae, but the others were assigned to a clade level. Phylogenetic resolution derived from ITS2 is not useful in defining Fusarium species. Recently, Fusarium-specific primers targeting translation elongation factor 1 (TEF1) were evaluated and successfully applied to analyze Fusarium communities in soil and plant material [57].
The taxonomy of Fusarium fungi is confusing and various classification systems have been proposed [58]. For Fusarium, chemotaxonomy is considered a supplement to traditional morphology-based taxonomy. Several fungal genes involved in trichothecene and enniatins biosynthesis have been defined and used for development of molecular assays aimed at identification. In spite of the ITS sequences used in our analysis, the results strongly suggested the division of fungi based on their ability to produce metabolites. In the future, this will provide an opportunity to predict the severity of grain contamination by some mycotoxins according to the number of certain identified OTUs.
Fusarium avenaceum, F. poae, F. tricinctum, and F. sporotrichioides were the most abundant representatives of the genus. They are the typical pathogens of barley in northwestern Russia [59,60]. Most likely, multi-copy otu259 discovered by DNA metabarcoding is associated with F. avenaceum, which occurred frequently on the barley grain.
Both methods revealed pathogenic fungi from Pleosporaceae: Bipolaris and Pyrenophora. Those fungi have different patterns of appearance through the cropping seasons. DNA metabarcoding demonstrated higher sensitivity. Pyrenophora sp. colonies were not recovered in 2015 at all, but several respective reads were obtained for 7 out of 10 samples.
Davidiella (Cladosporium) associated reads were abundant in DNA metabarcoding assay in non-sterilized samples but only single colonies were detected in the agar plate-based test. Underestimation of relative abundance of the fungus in the latter case can be result of two reasons: rapidly spreading colonies suppress or mask slowly growing fungi, and infected individual grains contain not uniform quantity of fungal biomass that appear as insufficient correlation between the  number of infected grains and the amount of fungal DNA in the whole sample. Four fungal genera revealed by only DNA metabarcoding contained agents of cereal diseases (Neoascochyta, Botrytis, Microdochium, and Phaeosphaeria). The first three taxa were represented by solitary reads. Phaeosphaeria (otu106 and otu215) were found in 14 of 20 samples. In 2015, in surface-sterilized samples, the relative abundance of Phaeosphaeria reads varied between 11 and 36%. Sequences of otu106 had the closest similarity (99%) with representative sequences of Parastagonospora avenae (Septoria avenae or Stagonospora avenae), widespread fungus causing leaf blotch of barley and some other cereals, and Parastagonospora poagena, a recently described fungus from Poa sp. [61,62]. Less abundant OTU, otu215, had a similarity of 98%, with several Phaeosphaeria species and with some unidentified endophytes.

CONCLUSION
DNA metabarcoding, based on high-throughput sequencing, is a sensitive and powerful method of grain mycobiome analysis that provides large amounts of data. However, at this time, not all fungi can be identified to species level by molecular markers, especially by rDNA sequences. In spite of universality, rDNA has a limitation as a taxonomic marker. The resolution of the ITS sequence-based method is not enough to differentiate many fungal species. For instance, many Fusarium species have nonorthologous copies of ITS2. Many other important plant pathogenic and toxigenic fungi also can be identified up to genus level, but that is not always informative in the framework of mycological seed expertise. Erroneous and chimerical sequences, as well as the lack of reference sequences of many species, still limits wide application of NGS-based technologies in biodiversity studies.
The most complete and credible results can be obtained when several approaches are implemented simultaneously. Combining the results of DNA metabarcoding and traditional culture-plating assay allowed us to revise the diversity of fungi colonizing on the surface of and inside barley grains in Leningradskaya oblast (northwest Russia).
Fungal species diversity of barley grain revealed by DNA metabarcoding formally exceeded the traditional microbiological culture-based agar plating: 43 operational taxonomic units (OTUs) vs. 17 taxa of genus or species level. DNA metabarcoding assay allowed seven ascomycete taxa to be added to the total list. Of those additional taxa, only Phaeosphaeria was abundant internal fungus. Seventeen OTUs belonging mainly to surface-seed-borne, yeastlike Basidiomycetes were completely outside the scope of traditional analysis. Meanwhile, routine mycological analysis, in contrast to DNA metabarcoding, resulted in precise identification of practically important Fusarium species. On the other side, due to DNA analysis, one Alternaria taxon was reidentified as Alternaria section Pseudoalternaria instead of section Infectoriae.

CONTRIBUTION
Authors are equally related to the writing of the manuscript and are equally responsible for plagiarism.

CONFLICT OF INTEREST
The authors declare that there is no conflict of interest regarding the publication of this article.