This study used four soybean varieties: two FM-susceptible varieties, “ND12” and “GX1,” which have a yellow seed coat and are conventional cultivars in Southwestern China, and two highly resistant germplasms, “C103” and “2162,” which have a black seed coat and are grown in Sichuan Province. The soybean plants were grown in pots in the experimental field of Sichuan Agricultural University in Ya’an, China (103°00′E, 30°08′N). Six seeds were sown per pot, and 2 weeks after sowing, the seedlings were thinned to three plants per pot. The soybean plants were transferred from the field to a solar greenhouse approximately 5 days before growth stage R7 (beginning maturity). Plants in the greenhouse were exposed to day/night temperatures of 21/13°C and humidity of 85–100% for 7 days during the remainder of the seed development and maturation period, according to Keigley and Mullen (1986) with some modifications. The soybean fruits (seeds and pods) were harvested at stage R8 (full maturity), when mildew covered the plants. Three biological replicates of each soybean variety were evaluated. The fruits from the middle portion of the plants were collected and immediately stored in sterile plastic tubs reported to avoid aerial contamination. The tubs were sealed and stored at -80°C until further analysis.

The mold levels of the seeds and pods of these four varieties were determined by estimating the percentage of the seed/pod area that was mildewed. The details of this method were described in our previous study (Liu et al., 2016). Briefly, mildew rate= × 100 (f, mold seed/pod levels; n, number of surveys). Mildew index = × 100 (f_i, mold seed/pod levels of various grades; s_i, corresponding mold grades; n, number of surveys; s_mas, top mold grade). These mildew indexes were subjected to regression analyses using all the experimental data to construct mathematical models using Microsoft Excel 2016 and SPSS software.

Samples were collected from the same five surfaces of soybean fruits from three randomly chosen sterile plastic tubs. Total genomic DNA from 5 g of soybean samples was extracted directly using the FastDNA^® spin kit (MP Bio, Santa Ana, CA, USA) according to the manufacturer’s instructions. The DNA was diluted to 1 ng/μL using sterile water and amplified using the primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2 (5′-GCTGCGTTCTTCATCGATGC-3′), which target the ITS1 region of the fungal rRNA gene (Krumins et al., 2009). The PCR amplifications were conducted using TransGen AP221-02 (TransStart Fast Pfu DNA Polymerase) and an ABI GeneAmp^® 9700. Amplification was assessed by analyzing 3 μL of product on 2% agarose gels. Samples with an intense band between 100 and 500 bp were chosen for further experiments. Sequencing libraries were generated using the NEB Next^® Ultra^TM DNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations, and index codes were added. The library quality was assessed using the QuantiFluor^TM-ST Blue fluorescence quantitative system (Promega), and the library was sequenced using the Illumina MiSeq platform at Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China, according to standard protocols.

Operational taxonomic unit (OTU) clustering, species annotation and diversity analysis were conducted according to Zhang et al. (2016). OTUs were assigned using USEARCH (version 7.1¹) with a threshold of 97% pairwise identity (Zhang et al., 2016). The normalized OTU abundance information was utilized for alpha (single sample analysis) and beta (multivariate sample community) diversity analyses. Metrics of community richness included the Chao1 and ACE estimators. Community diversity was assessed using the Shannon and Simpson indices. Coverage, indicating sequencing depth, was evaluated using Mothur version v.1.30.1 (Schloss et al., 2011). For beta diversity, principal component analysis (PCA) was applied for the multivariate examination of OTUs (Wang et al., 2012), and Bray-Curtis dissimilarity (Srinivasan et al., 2012) was calculated using the R package (Jami et al., 2013).

Fungi on FM-damaged soybean fruits were isolated using the suspension plating method (Zhou et al., 2014). In brief, 2 g of each sample was added to 50 mL of sterilized water, and the mixture was rotated for 1 day to produce a 10^-1 slurry (w/v); the suspension was then diluted to final concentrations of 10^-2, 10^-3, 10^-4, and 10^-5. A 0.1-mL aliquot of the 10^-3, 10^-4, and 10^-5 suspensions was placed in a 90-mm-diameter Petri dish containing 20 mL of Martin’s rose bengal agar (henceforth Martin medium) with streptomycin (at approximately 40°C) and spread evenly using sterile coated rods. Sterilized water was used as the control. The plates were kept in the dark at 37°C for 1–2 days. Single fungal colonies with morphological differences were selected from the mixed populations, sub-cultured and purified using the streak plate method (Bondici et al., 2013). Distinct colonies, based on morphological type and color, were selected for isolation and subsequent taxonomic and molecular identification. All isolates were sub-cultured and initially grouped into morphotypes based on morphological characters. The fungi were first identified to the genus or species level based on the morphological characteristics reported in the original taxonomic literature (Chakrabarti and Ghosal, 1989; Amoah et al., 1995) and then further identified based on ITS sequence comparison.

To confirm the pathogenicity of the isolated fungi, preliminary pathogenicity testing of the fungal isolates was performed according to Koch’s postulates. Each isolated fungus was first grown on Martin medium in the dark at 28°C for 7 days until reaching sporulation. Each fungus was mixed with sterilized water to produce a fungal suspension of 2 × 10⁶ spores/mL. Healthy seeds of the FM-susceptible variety “ND12” with no defects were infected with the fungal suspension in 5-cm Petri dishes containing moistened filter paper to encourage spore germination on the seed surface. Sterilized water was used as a blank control. All treatments were performed in triplicate. The Petri dishes were sealed with film to maintain a suitable relative humidity (98∼100%) and incubated at 25–28°C in darkness for 7 days. The samples were inspected daily until the 7th day, at which point the levels of mold on the infected seeds were determined according to the method of Li et al. (2005) with slight modifications. These procedures were described in our previous study (Deng et al., 2015b). The pathogenicity of each fungus was confirmed using a previously reported fungal recovery procedure (Ali-Shtayeh et al., 2003).

The isolated pure pathogenic fungi with different morphotypes were grown on Martin medium in the dark at 28°C for 5 days until reaching sporulation. The mycelia were placed in a sterilized mortar, immediately frozen in liquid nitrogen, and ground into a dry powder for further analysis. Genomic DNA was extracted from pathogenic fungi as described by Zhou et al. (2014). The molecular identities of strains were determined by PCR amplification and sequencing of the nuclear ribosomal DNA ITS region using the following primers: ITS1 (5′-TCCGTAGGTGAACCTGCGG-3′) and ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) (Labbé et al., 2008). Each PCR mixture (50 μL) contained 2.5 μL of genomic DNA, 2.0 μL of each primer (10 μmol/L), 35.5 μL of sterile deionized water, 0.5 μL of Taq PCR Master mix, 5.0 μL of 10× PCR buffer, and 2.5 μL of dNTPs (2.5 mmol/L). The PCR amplification program consisted of 94°C for 5 min, followed by 35 cycles of 94°C for 40 s, 54°C for 1 min and 72°C for 1 min, with a final extension of 72°C for 10 min. The PCR products were sequenced by Sangon Biotech (Shanghai) Co., Ltd., using an ABI-PRISM 3730 automated sequencer (Applied Biosystems, USA). The obtained ITS sequences were compared using BLAST analysis of the NCBI GenBank database.

After denoising and chimera detection, we obtained 62,382 sequences, which were then clustered into 164 OTUs with 97% sequence similarity. The numbers of OTUs detected for the tested soybean varieties “2162,” “C103,” “GX1,” and “ND12” were 146, 85, 71, and 93, respectively (Table 1). As indicated by the rarefaction curves shown in Figure 1, the number of observed OTUs (Y-axis) increased exponentially with the increasing randomness of ITS1 sequences (X-axis) but then tended to plateau. Although the number of observed OTUs was slightly lower than estimated based on Chao1 and ACE, these values were close to the estimates. The values of the Shannon and Simpson community diversity indexes were 1.74–1.89 and 0.29–0.37, respectively. The coverage data were good, with values higher than 0.9996, indicating sufficient sequencing depth. The above results indicate that the sequencing data are reasonable and that the number of sequences analyzed can sufficiently represent the fungal diversity of all samples.

The fungal community structures at the phylum and genus level in the different soybean varieties are shown in Figure 2. As shown in Figure 2, the four FM-damaged soybean samples belonged to four phyla: three formally described fungal phyla, Ascomycota, Basidiomycota, and Chytridiomycota, and an unassigned fungal phylum that includes a total of 69 fungal taxa. More than 85% of the total fungal sequences belonged to Ascomycota and Basidiomycota (Figures 2A–D); and were mainly distributed among twelve fungal taxa, including Mycosphaerellaceae, Sordariomycetes, and Alternaria (Figures 2E–H).

The rDNA ITS sequences of FM infecting four soybean varieties were detected using Illumina high-throughput sequencing. Of these varieties, “2162” and “C103” are black soybeans, whereas “GX1” and “ND12” have yellow seed coats. As shown in Figure 2, the relative abundances of fungal genera differed among the varieties. Of the total fungal sequences found on the black soybean “2162,” 59.71% belonged to Mycosphaerellaceae and 12.14% to Sordariomycetes (Figure 2E). Similarly, Mycosphaerellaceae (13.0%) and Sordariomycetes (42.83%) were the dominant fungal taxa found on “C103” (Figure 2F). In contrast, the dominant fungal communities on the yellow soybean varieties “GX1” and “ND12” differed from those found on the black soybean fruits: Pleosporales (45.69%) and Nectriaceae (27.87%) dominated on “ND12,” and Gibberella (39.46%) and Alternaria (32.01%) were the dominant fungi on “GX1”. These results were confirmed by PCA (Figure 3A) and hierarchical clustering analysis (Figure 3B) based on the beta diversity of each fungal community. As shown in Figure 3A, based on the first principal component, the four tested soybean varieties were clearly separated, with the black soybeans (“2162” and “C103”) being separated from the yellow soybeans (“GX1” and “ND12”). In addition, “GX1” and “ND12” were separated along the second principal component, whereas “2162” and “C103” were not. A similar result was obtained by using hierarchical clustering analysis, with the two black soybean varieties clustering together but the yellow soybean varieties being grouped into two independent clusters (Figure 3B). These results indicate that the fungal diversity among FM-damaged soybean fruits was highly host specific, especially between varieties with different colors of seed coats.

To further investigate the common relationships among fungal communities of different soybean varieties, the OTU distribution was statistically analyzed, and these results are presented as a Venn diagram in Figure 4. For “2162” and “C103,” 77 OTUs (14+13+42+8) were amplified; 54 OTUs (0+42+12+0) were amplified for “GX1” and “ND12,” indicating a more diverse fungal community on black soybean fruits. Additionally, 42 OTUs were amplified from all four different soybean varieties. Further analysis of these 42 OTUs demonstrated that 37 of them belong to Ascomycota, mainly Alternaria, Gibberella, Hannaella, Colletotrichum, Periconia, Cladosporium, Epicoccum, Boeremia, and Fusarium (Supplementary Table S1).

The current study also determined the organ-specific mildew indexes of the seed (MIS) and pod (MIP) in different soybeans with different field mildew resistances. As shown in Table 1, the MIS values for the four tested soybean varieties decreased in the order GX1 > ND12 > 2162 > C103; in contrast, the MIP values decreased in the order 2162 > C103 > GX1 > ND12. To further clarify their relationship, regression analyses were conducted between fungal diversity and organ-specific mildew indexes. The relationship between MIP/MIS values and the indexes of community richness (ACE, Chao1) and OTUs indicated that quadratic models fit the experimental data well (Figures 5A,B). The community richness indexes (ACE, Chao1) and OTUs increased as MIP values increased based on a quadratic relationship with high determination coefficients (R² > 0.99) (Figure 5A). In contrast, these community richness indexes decreased as the MIS values increased based on a quadratic relationship with moderate determination coefficients (R² > 0.40) (Figure 5B).

Regression analyses were also conducted between organ-specific mildew indexes (MIP, MIS) and the community diversity indexes (Shannon, Simpson) (Figures 5C,D). Unlike the community richness indexes, a higher Shannon index means higher community diversity, while a higher Simpson index means lower community diversity (Wang et al., 2012). As shown in Figure 5C, with the increase in MIP values, the Shannon community diversity index decreased based on a quadratic relationship (R² = 0.2973), and the Simpson community diversity index increased based on a quadratic relationship with high determination coefficients (R² = 0.8886). These two regression curves indicated that community diversity decreased as MIP values increased. Furthermore, as shown in Figure 5D, with the increase in MIS values, the Shannon and Simpson community diversity indexes decreased based on a quadratic relationship with low determination coefficients (R² < 0.15).

Four single fungal colonies were isolated from the FM-damaged grains of the susceptible soybean variety “ND12”. These colonies were named F1, F2, F3, and F4; sub-cultured; and purified based on morphological characters. The morphological characters of the fungal colonies under study were recorded over 7 days of growth on Martin medium. Microscopic characteristics of conidia or hyphae were measured using an ocular micrometer (Figure 6). Based on the comparisons of colony and microscopic characteristics in classical studies, we identified these isolated single fungal colonies (F1, F2, F3, F4) as Aspergillus flavus (Figures 6A,B), A. niger (Figures 6C,D), Fusarium moniliforme (Figures 6E,F), and Penicillium (Figures 6G,H), respectively (Ding et al., 2003; Musmade et al., 2013).

To study genetic variation among the above four strains isolated from FM-damaged soybean grains, ITS regions were amplified using rDNA ITS1 and ITS2 primers and then sequenced. The lengths of the sequences among the selected strains ranged from 500 to 650 bp. To search for homologous sequences, the obtained ITS sequences were compared by using BLAST analysis with the NCBI GenBank database. The results of the ITS sequence comparison and molecular identification are shown in Table 2, confirming the results of the morphological identification. The F1, F2, F3, and F4 fungal strains were identified as A. flavus culture-collection MUM:10.220, A. niger isolate PJ-2, Fusarium sp. BAB-4364, and P. chrysogenum strain A096, respectively, with high similarity (Table 2).

The above four isolated candidate pathogenic fungi were re-inoculated onto the surfaces of grains of the FM-susceptible variety “ND12,” and the phenotype was observed after 7 days. As shown in Figure 7, the blank control (Figure 7CK) soybean grains were simply swollen, without obvious mildew characteristics. In contrast, the other re-inoculated grains showed obvious mold, with a high mildew index (Deng et al., 2015b) above 85. Compared with the pre-harvest naturally fungal-infected soybean seeds in the field (Figure 7FM), the grains re-inoculated with F. moniliforme (Figure 7F3) and P. chrysogenum showed a similar mildew phenotype of dark brown color, with a shrunken and cracked surface (Figure 7F4). These results indicate that the two fungi identified as F3 and F4 caused the same symptoms as experienced by the field FM-damaged soybean grains. Combining the above information regarding fungal community composition and diversity, Fusarium was determined to be the dominant strain on FM-damaged soybean grains; this fungus is widely parasitic to several varieties of soybean. Therefore, we conclude that the pathogenic fungus on the FM-damaged pre-harvested soybean grains was F. moniliforme.

Grain fungi can be divided into two large categories, field and storage molds, based on the location of the outbreak (Liu et al., 2016). Storage fungi can grow at low water activity, which occurs readily under a high-temperature (25–30°C) and humid environment. In contrast, FM infection typically occurs in the field when the environment is abnormally cold (temperature 13–21°C) and humid (humidity 85–100%). The ability of fungi to germinate, grow and sporulate on grain depends on the availability of water and temperature as well as on the intergranular gas composition (Zhang et al., 2014), and these factors influence grain spoilage by fungi. Different environmental conditions result in different fungal community compositions. For instance, storage mold on Zea mays is caused by A. flavus and Penicillium spp. (Pinto and Fonseca, 2006). To date, there have been many in-depth studies on post-harvest fungal ecology, especially with regards to the impact of fungal growth on stored grain. F. culmorum, A. ochraceus, Penicillium verrucosum, and Fusarium section Liseola, as well as their released fungal toxins, are sources of contaminants in stored cereal grains (Magan et al., 2003). Abnormally cold and wet climatic conditions occasionally result in FM outbreaks in the field before soybean harvest in Southwest China. FM reduces soybean yield and quality and has become the main obstacle to soybean production in Southwest China (Deng et al., 2015b). However, there is scant information to date regarding FM on soybean grains, and the pathogenic fungi have not been elucidated. The current study investigated the fungal diversity among FM-infected soybean grains as well as identifying the pathogenic fungal, filling an important gap in this research.

In the current research, we surveyed the organ-specific mildew indexes of four soybean varieties with yellow and black seed coat colors. Our results indicated that the MIP values were higher than the MIS values, which implies that FM primarily occurred on the pod. To confirm this hypothesis, we analyzed the relationship between organ-specific mildew survey data and fungal diversity indexes for different varieties of soybean. The results revealed a close relationship between the FM fungal diversity and the organ-specific mildew indexes of soybean fruit. In particular, there were high determination coefficients (R² > 0.99) in the quadratic relationships between the community richness indexes and MIP values. With increases in the MIP, the fungal community richness increased sharply. In stark contrast, the community richness decreased as MIS increased.

Moreover, similar to our previous research, soybean seeds with yellow coats mildewed more easily than did black soybean seeds, indicating that soybean coat color has some type of special relationship with resistance to seed mildew in the field (Deng et al., 2015b). However, this was a speculation without direct evidence to confirm the function of the black seed coat. Our current research provides another possible mechanism for FM resistance. Although black soybean seeds were not more easily infected by FM, their corresponding pods were seriously damaged. In contrast, the opposite results were observed for yellow soybeans, with higher MIS but lower MIP values. It is plausible to suggest that pod mildew on soybeans may mitigate FM damage to the seed in a trade-off between pod and seed mildew.

Our recent research demonstrated a potentially important function of soybean pods in FM resistance at the phytochemical level (Liu et al., 2016). With increases in the contents of infection components, such as proteins, carbohydrates, and fatty acids, the MIP increased and the MIS decreased, demonstrating indirectly that pod mildew can help mitigate the corresponding seed mildew to a certain extent. Seed pods are the outermost barrier between legume seeds and pests and pathogens; our current research provides further evidence of the protective effect of soybean pods in FM resistance.

Additionally, the results of the present study also suggest that the richness and diversity of the FM fungal community were highly host specific among different soybean varieties. Of particular interest was another trade-off between richness and diversity of the FM fungal community. For example, although the FM community richness increased with MIP values in black soybean pods, the FM community diversity in black soybean pods was lower than that in yellow soybean pods. Our current study indicated that the fungal community was more diverse on black soybean seeds, suggesting the existence of antagonistic fungi, although more evidence is needed to confirm this possibility.

In our study, four fungal strains were isolated from the FM-damaged grains of the susceptible soybean variety. The integrated use of multiple methods, including morphological and molecular identification and re-inoculation, confirmed that F. moniliforme was the dominant FM pathogenic fungus among our samples, and we identified the strain as Fusarium sp. BAB-4364. Various Fusarium species not only cause seed rot and pre-emergence and post-emergence seedling damage but also result in wilting during the growth of numerous species of plants, including various cereals, vegetables, pulses, and oilseed crops (Muzahid et al., 2010). F. moniliforme is one of the most important pathogens causing several damaging crop diseases and producing enormous losses; of particular note are stalk rot and ear rot disease of maize, which have now achieved great importance in maize-growing areas (Musmade et al., 2013). Indeed, F. moniliforme produces mycotoxins, which are potential health hazards to humans and animals that frequently consume maize and maize products (Abeni et al., 2014). Species of Fusarium known to be seed-borne pathogens cause enormous losses to crops worldwide, and these fungi are frequently present on plant seeds; some are serious pathogens, whereas others are weakly pathogenic and still others merely saprophytes (Muzahid et al., 2010). Some species of Fusarium are seed-borne pathogens of soybean, especially F. moniliforme, which is responsible for wilt, root rot and foot rot (Muzahid et al., 2010); other species include F. anguioides, F. equiseti, F. fusarioides, F. moniliforme, and F. poae (Muzahid et al., 2010). However, reports of soybean fruit rot resulting from Fusarium are rare.

Additionally, Fusarium has been described as a latent pathogen, one that can function in both asymptomatic and symptomatic infections (Impullitti and Malvick, 2014). The slightly FM-damaged soybean fruit clearly displayed mildew symptoms on the pod surface, while the seed appeared healthy. However, during asymptomatic infection, the hyphae of Fusarium may already colonize intercellular spaces. When soybean plants suffer serious FM damage, which is a symptomatic infection, the hyphae of Fusarium may colonize plant tissues both inter- and intracellularly, and the biomass of the infected site will be reduced. In this situation, the seriously FM-damaged soybean seed exhibited a mildew phenotype with a dark brown color and a shrunken and cracked surface. The same symptoms were observed in soybean seeds re-inoculated with F. moniliforme (Deng et al., 2015a). In conclusion, the seriously FM-damaged soybean pod was more likely to be a symptomatic infection, while the slightly FM-damaged soybean seed was more likely to be symptomless endophytic colonization.