The grape samples were collected from six prime Marselan producing regions in China, namely, Ningxia’s Helan Mountain’s Eastern Region (NX), Fangshan District in Beijing (FS), Huailai County in Hebei Province (HL), Xiangning County in Shanxi province (SX), Penglai City in Shandong Province (YT), and Changli County in Hebei Province (CL) in 2017 vintage. Among them, the grapevines in the SX and YT could overwinter without soil burial, but the grapevines in FS, HL, CL, and NX had to be buried with soil for protection against the winter chill (Wang et al., 2018). For each region, three of the most representative wineries (vineyards) were considered for sampling, except for SX and CL, in which only one and two wineries, respectively, had planted Marselan grapevines (Supplementary Table 1).

For all regions, the Marselan grape samples were picked on a sunny morning within 3 days of harvest. Private wine producers authorized the sampling, and the field study did not involve any endangered or protected species. For each winery, three bags of healthy and undamaged grapes, about 1 kg each, were respectively, collected using sterilized garden scissors. To ensure the representativeness of the sampling, multiple bunches of grape samples were randomly picked from different positions in the vineyards, containing different ranks and orientations. These samples were placed into sterile plastic bags and transported to the laboratory chilled on ice within 1 day (Pinto et al., 2015).

Here, 45 bags of grape samples from fifteen wineries (three parallel samples per winery) were aseptically destemmed, and hand-squeezed in a clean bench, after which 350 mL of grape must and pomace from each sample were placed into sterile 500 mL jars, sealed with sterile sealing films. The physical and chemical indexes of the initial must is shown in Supplementary Table 2. Stationary spontaneous fermentation was performed at room temperature (controlled at 25 ± 2°C) (Portillo and Mas, 2016).

According to the pre-experiment of 2016 and the fermentation process of 2017, the Marselan must and krausen samples were collected on days 0 (must), 3, 5, and 8 for microbial diversity analysis. It should be noted that the fermented samples meaned all samples including must and krausen samples. During each stage, 9 mL must, or krausen samples were collected from corresponding spontaneous fermentation samples, placed into sterile 10 mL centrifuge tubes, and stored at −80°C for DNA extraction. In addition, the fermented jars were sealed with sterile plastic and sealing films to provide anaerobic fermentation conditions on day 3, and the fermented mash was transferred from the 500 mL jars to 250 mL sterilized triangular flasks to remove the skin residue on day 5. The fermentation process was monitored by measuring the °Brix every 2 days, and fermentation was considered completed when the °Brix displayed no change three consecutive times, the clarified fermented samples on day 20 was taken as the naturally fermented wine samples. Then, respective 50 mL must, and spontaneously fermented wine samples were stored at −20°C for physical and chemical index analysis. The fermentation process is shown in Supplementary Figure 1.

The genomic DNA of the fungal communities in the Marselan must and fermented samples were extracted using a modified CTAB method. The purity and concentration of the DNA were quantified with a nano spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States) and 1% (w/v) agarose gel. The fungal ITS2 gene PCR was performed via the ITS3-2024F (5′-GCATCGATGAAGAACGCAGC −3′) and ITS4-2409R (5′-TCCTCCGCTTATTGATATGC −3′) primers. The primary PCR reactions contained 20 ng DNA template, 2× AmpliTaq TSINGKE^® MasterMix (TsingKe Biotech, Beijing, China), 10 pmol of each primer. A secondary PCR to index the amplicons was performed with TaKaRa Taq DNA Polymerase. Amplification was conducted as follows: 94°C for 4 min, followed by 34 cycles of 94°C for 1 min, 57°C for 45 s, 72°C for 1 min, and a final extension at 72°C for 10 min. All the samples were amplified in triplicate, and no-template controls were included in all steps of the process (Zhang et al., 2019b). The Qiagen Gel Extraction Kit and PCR Clean-up (Qiagen, Germany) were employed to refine the PCR products. Furthermore, fluorometric evaluation using the Qubit 2.0 dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, United States) obtained the purity of the PCR mixture, after which sequencing analysis continued via the IonS5TMXL sequencing platform (IonS5TMXL, Thermo Fisher Scientific, United States) from Novogene, Beijing, China.

The pH and °Brix of the samples were measured with a HANNA 211 pH meter (HANNA, Padua, Italy) and a refractometer (ATAGO, Tokyo, Japan), respectively.

The glucose, fructose, ethanol, and glycerol content were determined with HPLC using a Waters 2414 RI Detector and a BIO-RAD Aminex HPX-87H resin-based column (300 mm × 7.8 mm) (Sun et al., 2015), and eluted with 5 mM H₂SO₄ at 55°C, 0.5 mL/min.

The acetic acid (AA), total acid (TA), polyphenol (PP), and anthocyanins (AC) content were analyzed with corresponding Randox kits on the Randox Monaco Analyzer (Randox, Monaco, United Kingdom).

Raw sequencing reads obtained from the IonS5TMXL platform were paired and pre-merged using FLASH software (Version 1.2.7), as well as filtered with the QIIME software (Version 1.7). All quality filtered sequencing reads were then clustered into operational taxonomic units (OTUs) with a minimum identity of 97%, by applying UPARSE software (Version 7.0) (Liu et al., 2018). Additionally, the OTU table underwent a series of filtering steps, including removing low-quality bases and chimeras, removing the OTUs with <3 counts across all samples, removing possible contaminants (mitochondrial and chloroplast sequences). The remaining OTUs samples were rarefied at a value equal to the median amount of all sequences (80160 sequences) to compensate for the uneven sequencing depth between the samples. For each fungal representative sequence (OTUs), taxonomy was assigned based on the UNITE fungal ITS database.

The Shannon index, Simpson index, abundance-based coverage estimator (ACE), chao1, and phylogenetic diversity (PD) whole tree values were calculated to compare the intra-group diversity (α-diversity). Moreover, Wilcoxon tests were used to evaluate the differences in average α-diversity indices between the groups (Blount et al., 2019). Furthermore, unconstrained principal coordinates analysis (PCoA) was conducted based on weighted UniFrac distance to evaluate the differences between the fungal communities of different regions and wineries, using the vegan package (version 2.5.3) in the R software (version 4.0.0). Based on the PCoA results, orthogonal partial least squares discriminant analysis (OPLS-DA) was performed to characterize the similarities and differences between the fungal microbiota of different regions using R software with the DiscriMiner package (version 6.3-73) (Shankar et al., 2013). Moreover, the multivariate analysis of variance (Adonis) and multi-response permutation procedure (MRPP) was conducted to identify the significant differences between the regions and wineries using the vegan R package (Cai et al., 2015).

Subsequently, the linear discriminant analysis (LDA) effect size (LEfSe) algorithm was performed to identify the representative fungal taxa of each region, by utilizing the Huttenhower Lab Galaxy Server (Liu et al., 2018). After the environmental factors (Gc and Pi) were screened via the variance inflation factor (VIF) and step selection using the vegan and dplyr (version 0.4.3) R package, detrended correspondence analysis (DCA) result determined to apply redundancy analysis (RDA) to explore the correlation between the environmental factors and fungal communities. Moreover, variation partitioning analysis (VPA) was performed to evaluate the contribution rate of each factor (Li et al., 2017). In addition, Spearman’s correlation analysis was applied to estimate the correlation between the environmental factors, the Pi variations of the spontaneously fermented wine (difference value with must), and the major fungal genus (Huang et al., 2019).

Other data were expressed as mean ± standard deviation, and the statistical significance between the groups was analyzed with one-way analysis of variance (ANOVA) using SPSS 19.0 software (SPSS Inc., Chicago, IL, United States).

To characterize the fungal communities throughout the Chinese Marselan wine-producing regions, 180 triplicated must, and spontaneously fermented samples were collected from 15 wineries (or vineyard blocks) in six regions. It is worth noting that eight samples were removed according to the abnormal Shannon index and species annotation, while six samples were eliminated since they did not meet the sequencing requirements (Supplementary Table 3). The remaining valid sequences in 166 samples with average ITS2 reads of 80,055 were clustered into 1611 OTUs with 97% similarity (Supplementary Data Sheet 1). Compared with the fungal reference database, taxonomic assignment revealed at least three fungal phyla (Ascomycota, Basidiomycota, and Mucoromycota), 22 classes, 66 orders, 148 families, and 306 genera after removing the samples assigned to bacterial taxa, in addition to some unknown groups, indicating the relative extent of uncharacterized fungi. The Ascomycota dominated 87.68% of all OTUs, followed by Basidiomycota at 4.28% (Supplementary Data Sheet 2). Aureobasidium, Alternaria, and Cladosporium were the dominant fungal genera in the must samples, accounting for 24.40, 17.29, and 14.82%, respectively. The mycobiota of the must also contained Colletotrichum, Rhodotorula, Metarhizium, Botrytis, Papiliotrema, Hanseniaspora, and other trace fungi with a relative content exceeding 2%, while the unclassified genus accounted for 7.97% (Figure 1 and Supplementary Table 4).

Considering the fermentation process, Aureobasidium, Hanseniaspora, Alternaria, Cladosporium, and Rhodotorula comprised the primary fungal genera, accounting for 24.88, 16.21, 9.16, 8.78, and 6.14%, respectively. Notably, the proportion of Saccharomyces throughout the fermentation process increased to 5.32%, while the unclassified genus accounted for 6.57%. The Marselan fermentation mycobiota contained Colletotrichum, Botrytis, Metarhizium, Aspergillus, Acremonium, and Papiliotrema accounting for more than 1% (Figure 1 and Supplementary Table 4). The rarefaction curves and Good’s coverage data of the Marselan must and fermented samples indicated that the ITS2 gene library was generally well-constructed (Supplementary Figure 2 and Supplementary Table 5).

Moreover, 283 fungal species were detected in the Marselan must samples of YT, which was the highest of all six regions, while SX displayed the least number at only 208 fungal species. The number of observed species in descending order was as follows: YT, HL, FS, CL, NX, and SX (Figure 2B). Although the variations in ACE, Chao1, the Shannon indices, and the Simpson index were generally consistent with the species number trends, FS, the third-highest on the list of fungal species, exhibited the lowest Simpson index value (Supplementary Table 5A). At the winery level, Chateau Rongzi (SX.rz) in SX presented the lowest number of fungal species, while the Winery Baihuagu (HL.bhg) in HL had the highest (Figure 2D and Supplementary Table 6A).

During the fermentation process, the number of fungal species in the six Marselan regions decreased slightly, but YT still had the highest number of species and NX the lowest, while CL, ranking fourth in observed species, demonstrated the minimum Simpson index value (Supplementary Figure 3B and Supplementary Table 5B). At winery level, Winery Yunkou (NX.yk) in NX displayed the least number of species at 208, while HL.bhg in HL exhibited the highest at 299 (Supplementary Figure 3D and Supplementary Table 6B).

The relative abundance diagram of the fungal communities in different regions indicated that whether it was must or fermented samples, the fungal microbiota composition of NX and HL regions were relatively similar, while YT and CL were relatively similar. And the FS and SX regions were obviously different from other regions, while FS contained a higher proportion of Colletotrichum, SX contains a higher proportion of Rhodotorula (Figure 2A and Supplementary Figure 3A). Furthermore, except for individual wineries in FS (FS.rdb) and YT (YT.tll) regions, the fungal structure among different wineries in the same regions were relatively similar for both must and fermented samples (Figure 2C and Supplementary Figure 3C).

As a visual and intuitive display, the PCoA method was conducted to explore the biogeographic distribution patterns of the fungal communities in different regions and wineries. The first two principal coordinates (PC) axis for the must account for 58.57% of the total fungal variance. Although extended distances were evident between the wineries within the FS and YT regions, the fungal communities displayed distinct clustering at the region level. SX was noticeably different from other regions, located at the left bottom of the third phenomenon, indicating the remarkable impact of the high Rhodotorula abundance. While the fungal compositions of NX and HL regions were similar displaying closer distances, mainly locating at the first phenomenon, the sample points of FS, YT, and CL regions spanned the second and third quadrants, and some sample points of YT and CL regions displayed a partial overlap (Figure 3A). Although not as clear as the boundaries between different regions in must, the fungal microbiota of the fermented samples also revealed obvious clustering phenomenon at the region level, except for the FS and YT regions (Figure 3C).

Through PCoA analysis, we could preliminarily confirm the regional microbial clustering landscape of fungal communities in the main Chinese Marselan wine-producing regions, supervised OPLS-DA analysis was further applied to verify the regional distribution pattern. The OPLS-DA plots display that the fungal microbiota exhibit the clustering phenomenon at the region level, and the differences among regions are more distinct. For must fungal communities, HL, NX and SX regions are closer, manifesting their relatively approximate fungal structure. FS and CL regions are relatively close, indicating that their fungal composition is relatively similar. And YT is located in the third quadrant alone, indicating its unique fungal structure (Figure 3B). The OPLS-DA plot of the fermented samples displayed the same regional distribution characteristics, but there were some overlapping sample points between different regions, revealing that the diversity of fungal microbiota between different regions decreased (Figure 3D).

Although the regional distribution pattern of fungal communities in the main Chinese Marselan wine-producing regions has been confirmed by the ordination analysis, the similarities and differences varied from different regions and wineries. Consequently, Adonis and MRPP were applied for difference quantitation to further assess the fungal community homeomorphism at different geographic scales. The Adonis and MRPP tests based on weighted UniFrac indicated that substantially significant differences (P = 0.001) were observed between the fungal compositions of the must from different regions (R²_Adonis = 0.560, A_MRPP = 0.296) and wineries (R²_Adonis = 0.830, A_MRPP = 0.518). During the fermentation process, extremely significant discrepancy were still evident in the fungal consortium of both different regions and wineries. However, the differences exhibited a decline (Region: R²_Adonis = 0.308, A_MRPP = 0.183, Winery: R²_Adonis = 0.478, A_MRPP = 0.287) (Table 1). Differences among the wineries in each region were also investigated, indicating significant discrepancy between the fungal composition of wineries within each region except CL. The Adonis of the wineries in CL displayed no significant differences based on weighted UniFrac distance (P = 0.4). During the development of fermentation, the differences between the fungal communities of the wineries in the different regions also decreased (Supplementary Table 7).

As we have confirmed the regional distribution pattern of fungal communities in the Chinese main Marselan wine-producing regions, LEfSe algorithm with a 4.0 LDA score threshold was applied to identify the discriminative fungal taxa in the Marselan must and fermented samples from the different regions. Among the fungal communities of Marselan must, 47 fungal species were verified as being differentially abundant in six regions. Of these, Filobasidiaceae (family), Phoma (genus), and Aspergillus (genus) were significantly enriched in NX; Mycosphaerellaceae (family), Colletotrichum (genus), and Lasiodiplodia (genus) were prominently abundant in FS; genus Alternaria and Filobasidium were significantly enriched in HL; genus Rhodotorula, Botrytis, and Monilinia were prominently abundant in SX; genus Cladosporium, Metarhizium, and Acremonium were significantly enriched in YT; genus Meyerozyma and Papiliotrema were prominently abundant in CL (Figure 4).

Considering the Marselan fermented samples, 36 fungal taxa were verified as being differentially abundant in the six regions, which were slightly less than the Marselan must numbers. Among them, genus Aspergillus was the only strain exhibiting a significantly higher abundance in NX; genus Colletotrichum and Lasiodiplodia were significantly enriched in FS; genus Alternaria and Filobasidium were prominently abundant in HL; genus Rhodotorula, Botrytis, and Monilinia were significantly enriched in SX; genus Cladosporium, Metarhizium, and Acremonium were prominently abundant in YT; genus Aureobasidium, Papiliotrema, and Phoma were significantly enriched in CL (Supplementary Figure 4). These results indicated that of all the sampling regions used in this study, the Marselan grape regions could be effectively distinguished by the demonstrative microorganisms at different biological classification levels, from order to genus.

Redundancy analysis was applied to explore the relationship between the Marselan must fungal communities, and the environmental factors (Gc and Pi). The geoclimatic data of 2017 for the six Marselan regions were supplied by the Wuhan Donghu Big Data Trading Center Co., Ltd. (Table 2)¹. After removing all the autocorrelation variables, the RDA plot indicated that ten of the 17 tested indexes significantly correlated with the microbial composition. These included six Gc, namely solar radiation (SR), average temperature (AT), average soil temperature (AST), LAT, LNT, and evaporation capacity (ET), as well as four Pi, namely total sugar (TS), PP, TA, and pH. While these environmental factors accounted for 58.17% of the taxonomic information, the first two axis explained 39.60% of the fungal variation (Figure 5A). Further VPA indicated that the Gc was principally responsible for shaping the fungal communities of the Marselan must, accounting for 36.20% of the fungal variation, while the Pi explained 10.20% of the total variability. Furthermore, 41.80% of the community distribution was not constrained by these two groups of determining factors (Figure 5B).

Then, Spearman correlation analysis was used to further investigate the correlation between the variations in the major fungal genus and environmental factors. As such, the relative abundance of genus Saccharomycopsis, Colletotrichum, Acremonium, Metarhizium, Lasiodiplodia, and Penicillium were positively correlated with LNT, AT, and AST, and negatively with TS and PP. Moreover, the relative abundance of genus Rhodotorula, Alternaria, Filobasidium, and Naganishia displayed a significant positive correlation with SR, ET, and TS, and a negative correlation with LNT, AT, and AST, while various some other fungal genus were affected by the different Gc and Pi (Figure 6).

Although Marselan must generated more than 8% v/v alcohol via spontaneous fermentation (except for FS.leb, FS.rdb, HL.gz, and CL.cc), most of the wineries did not meet the residual sugar standard for dry wine of below 4 g/L, including NX.zh and NX.yk in NX, FS.leb, FS.lxb, and FS.rdb in FS, HL.bhg and HL.gz in HL, as well as CL.cc and CL.lgs in CL. Furthermore, the AA content in most of the wineries exceeded the legal standard of 1.2 g/L, included NX.bs, NX.zh, and NX.yk in NX, FS.leb and FS.rdb in FS, HL.bhg in HL, as well as CL.cc and CL.lgs in CL (Table 3, OIV, 2013).

Spearman correlation analysis was performed to investigate the correlation between the fungal composition of the Marselan fermented samples and the Pi variations of the spontaneously fermented wine (difference value with must), the results indicated that AC and PP content correlated positively with genus Filobasidium and Rhodotorula, and negatively with genus Colletotrichum and Acremonium. Notably, the fermentation rate (FR) was confirmed to correlate positively with genus Filobasidium and Rhodotorula, and negatively with genus Paramycosphaerella, Colletotrichum, and Acremonium, which was consistent with AC and PP content (Figure 7).