Analysis of the Fungal Community in Ziziphi Spinosae Semen through High-Throughput Sequencing

Ziziphi Spinosae Semen (ZSS) has been widely used in traditional Chinese medicine system for decades. Under proper humidity and temperature, ZSS is easily contaminated by fungi and mycotoxins during harvest, storage, and transport, thereby posing a considerable threat to consumer health. In this study, we first used the Illumina MiSeq PE250 platform and targeted the internal transcribed spacer 2 sequences to investigate the presence of fungi in moldy and normal ZSS samples collected from five producing areas in China. Results showed that all 14 samples tested were contaminated by fungi. Ascomycota was the dominant fungus at the phylum level, accounting for 64.36–99.74% of the fungal reads. At the genus level, Aspergillus, Candida, and Wallemia were the most predominant genera, with the relative abundances of 13.52–87.87%, 0.42–64.56%, and 0.06–34.31%, respectively. Meanwhile, 70 fungal taxa were identified at the species level. Among these taxa, three potential mycotoxin-producing fungi, namely, Aspergillus flavus, A. fumigatus, and Penicillium citrinum that account for 0.30–36.29%, 0.04–7.37%, and 0.01–0.80% of the fungal reads, respectively, were detected in all ZSS samples. Moreover, significant differences in fungal communities were observed in the moldy and normal ZSS samples. In conclusion, our results indicated that amplicon sequencing is feasible for the detection and analysis of the fungal community in the ZSS samples. This study used a new approach to survey the fungal contamination in herbal materials. This new approach can provide early warning for mycotoxin contamination in herbal materials, thereby ensuring drug efficacy and safety.


Introduction
Herbal medicines, which are commonly used to prevent, diagnose, and cure diseases, have played an important role in health care since ancient times [1]. However, many studies have described the occurrence of mycotoxins in medicinal plants and herbal medicines from various countries [2][3][4], thereby attracting considerable attention worldwide due to drug efficacy and safety. For instance, 58 (8.29%) and 17 (2.43%) of the 700 herbal medicine samples in South Korea are aflatoxin B 1 (AFB 1 )and total aflatoxin (AF)-positive, respectively, and the AFB 1 (up to 73.27 mg/kg) and total aflatoxin used in fungal ecology studies [30][31][32], thereby providing a new prospect for the detection of fungal diversity in herbal medicines.
Here, we first investigated the presence of fungi in ZSS by using the Illumina MiSeq PE250 platform and demonstrated the diversity of fungal contamination in the ZSS samples to provide early warning for mycotoxin contamination in ZSS.

Analyses of the Diversity of Fungal Communities in the ZSS Samples
All 14 ZSS samples were successfully amplified by PCR for the internal transcribed spacer 2 (ITS2) sequences. A total of 912,941 ITS2 sequences that were longer than 200 bp were produced after excluding the chimeric sequences. The optimized sequences were divided into 210 operational taxonomic units (OTUs) after cluster analysis. Table S1 shows the calculated Chao 1, Shannon, and Good's coverage to measure the α-diversity. High Chao 1 value demonstrated a large variation of species in each sample. High Shannon value indicated high community diversity in each sample. Meanwhile, the results of Good's coverage, which is an estimator of sampling completeness, indicated good overall sampling with levels of >99.8%. Rarefaction curve analysis showed that all samples were almost parallel to the x-axis, thereby indicating that the obtained reads were sufficient to represent the overall fungal diversity ( Figure S1). The fungal community in the WM group also had higher Chao 1 estimate and Shannon index than those of the FM group (Figure 1a). The number of unique and common OTUs for the two groups are shown in a Venn diagram (Figure 1b). The results showed that the WM group possessed more unique OTUs than the FM group.
For β-diversity, principal coordinate analysis (PCoA) showed that samples were clustered according to the presence of mold. This result was consistent with the results of hierarchical clustering analysis (Figure 1c-d). Significant differences were reported between the FM and WM groups (analysis of similarity, ANOSIM, R = 0.524 and p = 0.002; Figure S2a). By contrast, the samples from the five producing areas showed an insignificant difference ( Figure S2b).   For β-diversity, principal coordinate analysis (PCoA) showed that samples were clustered according to the presence of mold. This result was consistent with the results of hierarchical clustering analysis (Figure 1c-d). Significant differences were reported between the FM and WM groups (analysis of similarity, ANOSIM, R = 0.524 and p = 0.002; Figure S2.a). By contrast, the samples from the five producing areas showed an insignificant difference ( Figure S2.b).

Fungal Community Composition in the ZSS Samples
Three fungal phyla, namely, Ascomycota, Basidiomycota, and Mucoromycota, were identified in all ZSS samples. Among the three phyla, Ascomycota was the dominant fungus, accounting for 64.36%-99.74% of the fungal reads. Meanwhile, other fungal phyla were detected with extremely low relative abundances (0%-0.33%, Figure 2a).
Further taxonomical classification detected 61 genera, and the genera with high relative abundance are shown in Figure 2b. Among these genera, Aspergillus, Candida, and Wallemia were the most common, with the relative abundances of 13.52%-87.87%, 0.42%-64.56%, and 0.06%-34.31%, respectively. The differences in the fungal community compositions between the FM and WM groups were significant. The relative abundances of the Aspergillus species in the FM group (42.66%-87.87%)

Fungal Community Composition in the ZSS Samples
Three fungal phyla, namely, Ascomycota, Basidiomycota, and Mucoromycota, were identified in all ZSS samples. Among the three phyla, Ascomycota was the dominant fungus, accounting for 64.36-99.74% of the fungal reads. Meanwhile, other fungal phyla were detected with extremely low relative abundances (0-0.33%, Figure 2a).
Further taxonomical classification detected 61 genera, and the genera with high relative abundance are shown in Figure 2b. Among these genera, Aspergillus, Candida, and Wallemia were the most common, with the relative abundances of 13.52-87.87%, 0.42-64.56%, and 0.06-34.31%, respectively. The differences in the fungal community compositions between the FM and WM groups were significant. The relative abundances of the Aspergillus species in the FM group (42.66-87.87%) were higher than those of the Aspergillus species in the WM group (13.52 to 49.67%). Meanwhile, the relative abundances of the Candida species in the FM group (0.42-23.95%) were relatively lower than those of the Candida species in the WM group (3.13-64.56%). The linear discriminant analysis effect size (LEfSe) algorithm, which was used to identify the different relative abundances of the fungal taxa in the two groups, showed the same results ( Figure 2c). Compared with the WM group, Aspergillus was the most dominant genus in the FM group. Meanwhile, the WM group possessed a much higher proportion of Candida, Meyerozyma, and Botryosphaeria species than those in the FM group.
A total of 162 fungal taxa were identified in the ZSS samples. Among these fungal taxa, 70 can be identified at the species level, while the remaining 92 can be resolved at the genus level or higher. Among the 70 accurately differentiated species, potential mycotoxin-producing fungi, namely, A. flavus, A. fumigatus, and Penicillium citrinum that account for 0.30-36.29%, 0.04-7.37%, and 0.01-0.80% of the fungal reads, respectively, were detected in all ZSS samples (Figure 2d). A higher percentage of potential mycotoxin-producing fungi was also detected in the WM group (1.14-36.37%) than that in the FM group (0.61-5.03%).
than those of the Candida species in the WM group (3.13%-64.56%). The linear discriminant analysis effect size (LEfSe) algorithm, which was used to identify the different relative abundances of the fungal taxa in the two groups, showed the same results (Figure 2c). Compared with the WM group, Aspergillus was the most dominant genus in the FM group. Meanwhile, the WM group possessed a much higher proportion of Candida, Meyerozyma, and Botryosphaeria species than those in the FM group. A total of 162 fungal taxa were identified in the ZSS samples. Among these fungal taxa, 70 can be identified at the species level, while the remaining 92 can be resolved at the genus level or higher. Among the 70 accurately differentiated species, potential mycotoxin-producing fungi, namely, A. flavus, A. fumigatus, and Penicillium citrinum that account for 0.30%-36.29%, 0.04%-7.37%, and 0.01%-0.80% of the fungal reads, respectively, were detected in all ZSS samples (Figure 2d). A higher percentage of potential mycotoxin-producing fungi was also detected in the WM group (1.14%-36.37%) than that in the FM group (0.61%-5.03%).

Prevalence of Fungal Contamination in the Commercial ZSS Samples
Fungi are common natural contaminants of herbal medicines due to their widespread distribution in the air. Fungal contamination in medicinal herbs has received considerable concern from the public. An investigation [33] on the occurrence of fungi in medicinal herbs and spices from

Prevalence of Fungal Contamination in the Commercial ZSS Samples
Fungi are common natural contaminants of herbal medicines due to their widespread distribution in the air. Fungal contamination in medicinal herbs has received considerable concern from the public. An investigation [33] on the occurrence of fungi in medicinal herbs and spices from India showed that 92% samples are contaminated by fungi, and 47% of the contaminated samples exceed the permissible limits set by the World Health Organization. Kong et al. [4] also found that 14 functional foods and 10 spices from Chinese markets are infected by molds. Zheng et al. [34] also showed the widespread fungal contamination in 15 medicinal herbs collected from China. Meanwhile, Su et al. [35] found that 48 samples of eight root herbs from Chinese markets are all contaminated by fungi. In all these studies, Aspergillus and Penicillium are identified as the most common contaminants. In the present study, the results indicated that Aspergillus was the most abundant genus in the ZSS samples, followed by Candida, accounting for 13.52-87.87% and 0.42-64.56% of the total reads, respectively. In contrast to previous studies, Penicillium species were detected at relatively low abundances (<1%) in all ZSS samples, except for the LY01 sample (7.77%). The relationship between fungal species and herbal materials cannot be fully explained owing to the complicated contamination reasons [34,36]. Aspergillus and Penicillium species were the principal storage fungi [36]. The substrate composition and storage conditions, such as moisture content, aeration, and temperature, might be the factors to influence the fungal community [34,36].
Meanwhile, Aspergillus, Penicillium, Fusarium, and Alternaria are the most common contaminants in food and herbal medicines, which contain the majority of mycotoxin producers [3,16]. Potential toxigenic fungi, such as A. fumigatus, A. flavus, A. tubingensis, and A. aculeatus, were previously isolated in different herbal medicines [3,34,37]. In the present study, three potential mycotoxin-producing fungi, namely, A. flavus, A. fumigatus, and P. citrinum that are potential producers of AFB 1 , OTA, and citrinin, respectively, were identified in the ZSS samples. The high relative abundances of the potential mycotoxin-producing fungi ranging from 0.61-36.37% in ZSS poses a considerable threat to public health. The relative abundances of the potential toxigenic fungi in normal samples were much higher than those in the moldy samples. Therefore, evaluating the safety of herbal medicines according to the presence or absence of molds is unadvisable.
The differences in the fungal community structure between the FM and WM groups were significant. Candida species were dominant in the normal ZSS samples, whereas Aspergillus was the predominant genus in the presence of molds. Meanwhile, all samples were clustered according to the initial grouping in PCoA, thereby indicating the meaningfulness of grouping samples according to the presence of microscopic molds. However, the differences in the fungal community structures among the SX, HN, LN, HB, and SD groups were insignificant, thereby demonstrating that the presence or absence of molds affected the fungal community structure in ZSS samples to a greater extent than the producing areas.

DNA Marker Selection to Analyze the Fungal Community in ZSS
An effective molecular marker is the precondition for the application of DNA barcoding technique. Schoch et al. [38] evaluated six candidate DNA regions, namely, ITS, LSU, SSU, RPB1, RPB2, and MCM7, for fungal identification and proposed ITS as the primary barcode marker. The length of the ITS region differs substantially among the different fungal genera and species, which distort the community description due to the preferential amplification of shorter sequences, thereby resulting in a biased quantification of the taxon relative abundances in the HTS research [39,40]. ITS1 and ITS2, which were shorter than the complete ITS, were of appropriate length for HTS and have become the most commonly used markers for fungal diversity analyses. Considering the lower phylogenetic richness and fewer OTUs generated by ITS1 region than those of the ITS2 region, ITS2 is recommended for metabarcoding [41]. In the present study, the results showed that the ITS2 region exhibited good identification ability for fungi in the ZSS samples. A total of 162 fungal taxa were detected in the ZSS samples. Among these fungal taxa, 92 cannot be resolved at the species level mainly due to the following reasons. First, some information on species is lacking on the UNITE database. Second, the ITS2 region cannot provide an adequate variation information in distinguishing some sibling species in Aspergillus, Fusarium, Penicillium, and Candida [42][43][44][45].

Prospects of Applying Amplicon Sequencing for the Analysis of Fungal Diversity in Herbal Materials
The application of HTS platform in analyzing the abundance and richness of the microbial species with low abundances in environmental samples overcomes the isolation limitations of the culture-based approach. Most microorganisms cannot be cultivated using traditional cultivation techniques [46]. Xia et al. [47] identified 55 fungal genera from Chinese Cordyceps through Illumina Miseq sequencing that were not observed using culture-dependent methods. The amplicon sequencing platform is widely used for the analysis of microbial diversity in soil, outdoor air, food, and other environmental samples [32,[48][49][50][51]. The application of HTS platform in analyzing the fungal diversity in herbal materials has not been reported to date. In the present study, we utilized the amplicon sequencing platform and targeted the ITS2 sequences to investigate the fungal contamination in herbal materials for the first time. Fungal contamination in the ZSS samples is extremely common. Hence, all samples were infected. HTS platform is a good tool for investigating the occurrence of fungi, especially the potential toxigenic fungi, in herbal materials. The application of amplicon sequencing provides a new approach to analyze the fungal diversity in herbal materials, thereby providing an early warning for mycotoxin contamination in herbal materials.

Sampling
A total of 14 commercial ZSS samples were collected from Henan, Shanxi, Hebei, Liaoning, and Shandong, which are the five main ZSS producing provinces in China. The samples were divided into five groups, namely, SX, HN, LN, HB, and SD, by production areas. Among these commercial samples, six samples were affected with mildew due to inappropriate storage. Samples were also divided into two groups, namely, FM and WM, according to the presence or absence of macroscopic molds. The detailed information of the samples is listed in Table 1.

DNA Extraction
Approximately 3 g ZSS samples were transferred into a 15 mL sterilized centrifuge tube with 10 mL of sterilized water and shaken with a vortex mixer for 20 min. Then, the mixture was filtered through a single layer of sterile gauze, and the microorganisms were collected from the filtrate by centrifugation at 7830 rpm for 15 min (Centrifuge 5430 R, Eppendorf AG, Hamburg, Germany) [52]. The total DNA was extracted using the EZNA ® soil DNA kit (Omega Bio-tek., Inc., Norcross, GA, USA) according to the manufacturer's instructions.

PCR Amplification and HTS
ITS2 sequences were amplified with the primer pairs ITS3 (5 -GCATCGATGAAGAACGCAGC-3 ) and ITS4 (5 -TCCTCCGCTTATTGATATGC-3 ) [53]. PCR was performed with the following conditions: Initial denaturation at 95 • C for 5 min; 35 cycles of denaturation at 95 • C for 45 s, annealing at 58 • C for 50 s, and elongation at 72 • C for 45 s; and a final extension at 72 • C for 10 min. Amplifications were conducted for each sample in triplicate, and the PCR products were pooled to minimize the PCR bias. The PCR products were detected using electrophoresis on a 2% agarose gel and purified using the DNA gel extraction kit (Axygen, Union City, CA, USA). Subsequently, the PCR products were further identified on a 2% agarose gel and quantified using QuantiFluorTM-ST (Promega, Madison, WI, USA). Then, amplicons were sequenced using the Illumina MiSeq PE250 platform (Illumina, San Diego, CA, USA) by AuwiGene Technology Co., Ltd. (Beijing, China). Raw sequences are available in the Sequence Read Archive of the NCBI under the accession numbers SAMN10275058-SAMN10275071.

Sequence Analysis
Raw FASTQ files were demultiplexed, and their quality was filtered using the Quantitative Insights into Microbial Ecology (QIIME, version 1.8, http://qiime.org) software [54]. The reads were truncated at any site, thereby receiving an average quality score of <20 over a 50 bp sliding window. The primers were exactly matched, thereby allowing two nucleotide mismatches, and the reads containing ambiguous bases were removed. Only sequences that overlap longer than 10 bp were merged according to their overlap sequences. The sequences were clustered into OTUs at 97% sequence similarity by using UPARSE (version 7.1, http://drive5.com/uparse/) [55], and chimeras were removed using USEARCH (version 8.1.1861, Http://Www.Drive5.Com/Usearch/) [56]. OTUs were denominated at the kingdom, phylum, class, order, family, genus, and species levels according to the UNITE database [57]. To estimate the α-diversity, we calculated the three metrics, including Chao 1, Shannon, and Good's coverage, after the reads were normalized to the minimum reads (35,924 reads) in each sample. β-Diversity was measured using the weighted UniFrac and unweighted UniFrac distance matrices. PCoA was performed based on the weighted UniFrac distance matrices. Samples were also hierarchically clustered based on the unweighted Unifrac distances by using Unweighted Pair Group Method with Arithmetic Mean (UPGMA). LEfSe was applied to identify the differentially abundant taxa between the FM and WM groups [58]. Statistical differences among the sample groups were tested using ANOSIM (available through QIIME). Venn diagram and heat map were drawn using R tools.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-6651/10/12/ 494/s1, Figure S1: Rarefaction curves of OTUs for the ZSS samples, Figure S2: ANOSIM analysis of the ZSS samples. (a) The differences in the fungal community structures between the FM and WM groups were significant. (b) The differences in the fungal community structures among samples from the five producing areas were insignificant. Table S1: α-Diversity of the fungal community in the ZSS samples.

Conflicts of Interest:
The authors declare no conflict of interest.