LTD (Figure 1A), MTD (Figure 1B), and HTD (Figure 1C) samples were obtained from three representative Baijiu enterprises located in the core production regions of Shanxi, Sichuan, and Guizhou provinces, China, respectively (Figure 1D). These enterprises were selected for their long-term stable production practices and representativeness of typical Daqu production processes. In October 2023, mature Daqu bricks that had undergone natural fermentation and were stored for 3 months were randomly collected from the Daqu storage rooms of each enterprise. Three Daqu bricks were obtained from each site, yielding a total of nine samples. After collection, all samples were thoroughly ground, homogenized, and subdivided into sterile polyethylene bags. Each sample was divided into two portions: one was immediately subjected to microbial isolation and cultivation, while the other was stored at −80 °C for subsequent genomic DNA extraction and TGS (PacBio platform).

To obtain a comprehensive representation of fungal diversity, six types of culture media were employed: potato dextrose agar (PDA; 20% potato infusion, 2% glucose and 2% agar), PDA supplemented with 1% Daqu extract, bengal red agar (1% peptone, 2% glucose, 0.1% dipotassium phosphate, 0.005% magnesium sulfate, 0.0033% rose bengal and 2% agar), G25N agar (25% glycerol, 0.37% yeast extract, 0.10% dipotassium phosphate, 0.225% sodium salt, 0.0375% potassium chloride, 0.0375% magnesium sulfate heptahydrate, 0.00075% ferrous sulfate heptahydrate and 2% agar), yeast extract-peptone-dextrose agar (YPD; 1% yeast extract, 2% peptone, 2% glucose and 2% agar), and YPD supplemented with 0.2% (v/v) acetic acid. Microbial isolation and cultivation were performed using the conventional dilution plating method. Samples were initially subjected to a tenfold dilution and vortexed for 5 min to achieve thorough homogenization. Serial tenfold dilutions were then prepared up to 10⁻⁵. Aliquots of 100 μL from four dilution levels (10⁻², 10⁻³, 10⁻⁴, and 10⁻⁵) were spread onto six different media to maximize the recovery of diverse fungal isolates. PDA was incubated at 45 °C for the isolation of thermophilic filamentous fungi (Han et al., 2023). PDA supplemented with 1% Daqu extract and Bengal red agar were incubated at 30 °C for the isolation of mesophilic filamentous fungi (Li et al., 2022). G25N agar was used at 30 °C to isolate osmotolerant fungi (Pitt, 1973). YPD and acidified YPD (YPD + 0.2% acetic acid) were incubated at 30 °C for yeast isolation, with the acidified medium effectively suppressing filamentous fungi growth (Guan et al., 2019). The Daqu extract for PDA enrichment was prepared by boiling 10 g of crushed Daqu in 600 mL of distilled water, followed by filtration and volume adjustment to 1 L with distilled water to provide trace nutrients. All media contained 2% (w/v) agar and were sterilized at 115 °C for 30 min. After cooling to approximately 55 °C, chloramphenicol (0.2%, v/v; 100 mg/mL stock solution) was added to inhibit bacterial growth. The inoculated plates were incubated in an inverted position at the corresponding temperatures for 24–72 h until visible colonies appeared. Representative colonies were selected based on morphological characteristics and subcultured on the corresponding medium to obtain pure isolates. The purified strains were confirmed by streak plating and stored in 25% (v/v) sterile glycerol at −80 °C for long-term preservation.

Yeast genomic DNA was extracted following the method described by Wei et al. (2025a), while genomic DNA from filamentous fungi was extracted using a plant DNA magnetic bead extraction kit (Geneonbio, Changchun, China). The quality and concentration of the extracted DNA were assessed using 1% (w/v) agarose gel electrophoresis and a NanoDrop spectrophotometer (Thermo Scientific, USA). To identify yeast species, the D1/D2 region of 26S rDNA was amplified using the primers NL1 (5′-GCA TAT CGG TAA GCG GAG GAA AAG-3′) and NL4 (5′-GGT CCG TGT TTC AAG ACG G-3′) (O’donnell, 1993), while the ITS region was amplified using the primers ITS1 (5′-TCC GTA GGT GAA CCT GCG G-3′) and ITS4 (5′-TCC TCC GCT TAT TGA TAT GC-3′) for filamentous fungi identification (White et al., 1990). The amplification reaction volume was 25 μL, consisting of 22 μL of 1.1 × T3 Super PCR Mix (Tsingke Biotech Co., Beijing, China), 1 μL of DNA template, and 1 μL of each primer. The PCR cycling conditions were as follows: initial denaturation at 98 °C for 2 min; followed by 35 cycles of denaturation at 98 °C for 10 s, annealing at 55 °C for 15 s, and extension at 72 °C for 15 s; with a final extension at 72 °C for 5 min. To accurately identify certain filamentous fungi, the calmodulin gene (CaM) was amplified using the primers CMD5 (5′-CCG AGT ACA AGG ARG CCT TC-3′) and CMD6 (5′-CCG ATR GAG GTC ATR ACG TGG-3′) (Samson et al., 2014), and the β-tubulin gene was amplified using the primers Bt2a (5′-GGT AAC CAA ATC GGT GCT GCT TTC-3′) and Bt2b (5′-ACC CTC AGT GTA GTG ACC CTT GGC-3′) (Hubka and Kolarik, 2012). These amplifications were performed using the same reaction volume and PCR cycling conditions as described above. PCR products were visualized by electrophoresis on a 1% (w/v) agarose gel under UV light. The single band of the expected size was directly sequenced by Bio Magic Gene Company.¹ Sequences were analyzed using Chromas Lite v.2.1 and compared to the GenBank nucleotide database using BLAST.²

To complement the molecular identification results, representative strains of the dominant fungal species were subjected to systematic morphological characterization. Macroscopic colony features on individual plates were documented using a digital camera (Canon EOS 90D, Japan) or a stereomicroscope (Nikon SMZ25, Japan). For microscopic observations, high-resolution images of cells or hyphae were obtained using a scanning electron microscope (Zeiss Axio Imager A2, USA). All morphological features were compared with standard taxonomic references for yeasts (Kurtzman et al., 2011) and filamentous fungi (Pitt and Hocking, 2022) to validate the molecular identification results.

The total Genomic DNA was extracted from samples following a modified method described by Luo et al. (2023), with the tissue grinding step performed at 60 Hz for 2 min to ensure thorough cell disruption. DNA purity was evaluated by measuring the A260/A280 ratio using a NanoDrop™ One spectrophotometer (Thermo Fisher Scientific, USA), and DNA integrity was assessed by 1% (w/v) agarose gel electrophoresis. DNA samples were normalized to an equal concentration using nuclease-free water prior to PCR amplification. The full-length fungal ITS region was amplified using primers ITS1 and ITS4, with PCR conditions and cycling parameters as described by Wei et al. (2025b). Positive and negative controls were included to ensure amplification reliability, and no amplification or detectable contamination was observed in the negative controls. PCR products were purified with AMPure PB magnetic beads (Beckman Coulter, USA) to remove residual primers and impurities. Library preparation and sequencing were performed by Annoroad Gene Technology Co., Ltd. (Beijing, China) on the PacBio Sequel II platform.

PacBio amplicon sequencing data were processed following the standard SMRT analysis workflow to ensure the generation of high-quality sequences. Initially, barcode sequences of each sample were identified and demultiplexed using Lima v1.7.0 to obtain circular consensus sequences (CCS). Primer sequences were subsequently removed with Cutadapt v2.7 (Kechin et al., 2017), and only ITS fragments within the 300–1,000 bp length range were retained. To further improve sequence accuracy, fungal chimeric reads were identified and removed using the UNITE CHIME reference database (release October 16, 2022) (Nilsson et al., 2019). Subsequently, the unique sequences were denoised using the UNOISE3 algorithm to generate amplicon sequence variants (ASVs) with 100% nucleotide identity, with the minsize parameter set to 10. The resulting ASVs were taxonomically assigned by comparison against the UNITE UTAX reference dataset (released July 25, 2023) (Nilsson et al., 2019), using a 97% sequence similarity threshold as the criterion for species-level annotation. To further support taxonomic assignments, representative ASVs with high relative abundance were cross-validated using BLASTn searches against the NCBI nucleotide database.

The sequences obtained by Sanger sequencing were aligned with those of type strains from related species using MEGA v.7 (Kumar et al., 2016). A phylogenetic tree was subsequently constructed using the Neighbor-Joining (NJ) method in MEGA v.7 with 1,000 bootstrap replicates (Kumar et al., 2016; Lachance, 2022). The resulting phylogenetic tree was visualized using the Interactive Tree of Life (iTOL) online platform (Letunic and Bork, 2016). To ensure comparability among samples, the ASV table was rarefied to the minimum sequencing depth (3,052 reads per sample) using the rrarefy function. Subsequently, the fungal community composition and diversity were analyzed in R (v4.3.2) using the ggplot2, vegan, dplyr, and igraph packages. Alpha diversity indices (ACE, Chao1, Richness, and Shannon) were calculated using the diversity function, and differences among groups were assessed by analysis of variance (ANOVA) followed by Tukey’s HSD test. Furthermore, differential taxa were further identified using the MicrobiotaProcess package following the LEfSe statistical framework, which applies Kruskal–Wallis and Wilcoxon rank-sum tests; taxa with an LDA score > 2.0 were considered significant biomarkers. The results were visualized with ggtree and ggplot2 packages. In addition, shared and unique taxa among groups were visualized using Upset, Venn, and Venn network plots generated by the UpSetR, VennDiagram, and igraph packages, respectively. Macroscopic and microscopic images of the dominant fungal species were arranged and compiled using Adobe Illustrator 2020.

A total of 135 fungal isolates were obtained from the three types of Daqu based on culturomics, including 42 yeasts and 93 filamentous fungi. Phylogenetic analysis based on 26S rDNA sequences revealed that the yeast isolates belonged to six genera and seven species (Figure 2A), among which Saccharomycopsis fibuligera (12.6%, Supplementary Figure S1A), S. cerevisiae (6.7%, Supplementary Figure S1B), and Wickerhamomyces anomalus (5.9%, Supplementary Figure S1C) were the predominant species (Figure 2C). Meanwhile, ITS sequence-based phylogenetic analysis indicated that the filamentous fungi were classified into eight genera and 14 species (Figure 2B). Among these, five species within the genera Aspergillus, Thermomyces, and Paecilomyces could not be reliably resolved at the species level based solely on ITS sequences. Subsequent phylogenetic analyses using calmodulin (Supplementary Figure S2A) and β-tubulin (Supplementary Figure S2B) gene sequences further identified these isolates as Aspergillus chevalieri, A. latus, A. flavus/oryzae, Thermomyces lanuginosus, and Paecilomyces variotii. Overall, the dominant filamentous fungi across the three Daqu types were Lichtheimia ramosa (15.6%, Supplementary Figure S1E), A. flavus/oryzae (11.1%), Rhizomucor pusillus (8.9%, Supplementary Figure S1D) and Aspergillus fumigatus (4.4%, Supplementary Figure S1F) (Figure 2C). Comparative analysis of fungal composition among the three types of Daqu showed that filamentous fungi were most abundant in LTD, dominated by R. pusillus, followed by HTD, dominated by A. chevalieri; and least abundant in MTD, where L. ramosa was the major species. Yeasts were most frequently isolated from LTD and to a lesser extent from MTD, with S. fibuligera as the dominant species in both, whereas no culturable yeast was recovered from HTD.

In this study, fungal DNA from three types of Daqu was analyzed using full-length ITS sequences. After quality filtering, a total of 61,890 high-quality ITS sequence pairs were obtained, with an average of 6,876 reads per sample. These sequences were classified into 109 fungal ASVs, representing to three phyla, 33 genera, and 54 species. Overall, Ascomycota (54%) and Mucoromycota (44%) were the dominant phyla, with Mucoromycota prevailing in MTD, Ascomycota predominating in HTD, and LTD exhibiting an intermediate distribution pattern (Figure 3A). At the species level, L. ramosa and Pichia kudriavzevii were dominant in LTD and MTD, whereas Thermoascus crustaceus was the predominant species in HTD (Figure 3B). Alpha diversity indices were applied to comprehensively assess the richness and evenness of the fungal communities (Herrera et al., 2023). The ACE and Chao1 indices were used to estimate species richness, with ACE emphasizing differences in abundance between rare and common taxa, and Chao1 estimating the number of unobserved rare taxa (Deng et al., 2021). The Richness index reflects the total number of detected species, while the Shannon index incorporates both abundance and evenness, representing overall community diversity (Supriatna, 2018). Comparative analysis among the three types of Daqu revealed that LTD exhibited the highest ACE, Chao1, and Shannon indices, MTD showed the highest Richness index, whereas HTD displayed the lowest values across all four indices, although these differences were not statistically significant (Figures 3C–F). Principal component analysis (PCA) further demonstrated clear separation of fungal community compositions among the three Daqu types, with PC1 and PC2 explaining 52.4% of the total variance. Notably, the LTD samples showed the least within-group variation, whereas HTD exhibited the greatest dispersion (Figure 3G).

To elucidate differences in fungal community characteristics among the three types of Daqu, Linear discriminant analysis Effect Size (LEfSe) analysis was conducted based on third-generation sequencing data to identify potential fungal biomarkers (Figure 4). The results showed that four specific fungal taxa were identified in both LTD and HTD. In LTD, the identified biomarkers included Rhizopus oryzae, P. kudriavzevii, Pichia exigua, and Pichia fermentans, whereas in HTD they were T. crustaceus, A. chevalieri, P. brunneolus, and Aspergillus teporis. By comparison, MTD exhibited the fewest characteristic taxa, with only L. ramosa and W. anomalus detected as biomarkers. These findings highlight distinct fungal community compositions among the three Daqu types, providing insights into the fungal ecological characteristics shaped under different fermentation conditions.

To further explore the complementarity between sequencing- and culture-based approaches, we compared the fungal communities obtained by both methods across each Daqu type. The results revealed pronounced differences in fungal community structures among the different types of Daqu. Overall, TGS detected greater fungal richness, capturing a higher number of uncultured or difficult to culture taxa, whereas culturomics primarily recovered a limited set of cultivable dominant fungi. Specifically, TGS identified the largest number of characteristic species in HTD (29 species), while culturomics detected the most in MTD (seven species). Both approaches shared eight, eight, and four fungal species in LTD, MTD, and HTD, respectively (Figures 5B,D,F). Notably, L. ramosa and T. crustaceus were consistently detected as dominant species across all three types of Daqu by both methods, indicating their high ecological adaptability to varying temperature conditions (Figures 5A,C,E). In contrast, several dominant species were detected exclusively by only one of the two approaches. For instance, S. fibuligera and P. variotii were isolated exclusively through culturomics, whereas P. kudriavzevii was detected solely by TGS (Figures 5A,C,E). These differences suggest that integration both methods provide a more comprehensive understanding of the fungal community composition across different types of Daqu.

The composition and diversity of fungal communities in three types of Daqu were systematically compared using TGS and culturomics. The results revealed significant differences in fungal community structures among the different Daqu types, with both methods exhibiting a degree of consistency and complementarity. Based on the TGS results, 16 fungal species were shared across all three Daqu types (Figures 6A,B). HTD and MTD shared seven species, LTD and MTD shared five species, and no common fungal species were detected between LTD and HTD. The numbers of unique species identified in LTD, MTD, and HTD were seven, five, and 10, respectively, indicating distinct temperature-dependent characteristics of the fungal communities. According to the culturomics results (Figures 6C,D), eight fungal species were shared among the three Daqu types. HTD and MTD shared two species, LTD and MTD shared eight species, while no common species were detected between LTD and HTD. The numbers of unique species detected by culturomics in LTD, MTD, and HTD were two, three, and four, respectively. Overall, both approaches revealed distinct differences in the fungal community composition across the three Daqu types. Notably, TGS detected a greater number of both shared and unique species compared with culturomics, highlighting its superior resolution in uncovering fungal diversity in Daqu.