Characteristic Metabolites Drive the Self-Assembly of Microeukaryotic Communities during Spontaneous Fermentation of Icewine

Abstract

Wine can be regarded as a nutritional source for the human diet. It contains many nutrients such as vitamins and minerals, organic acids, flavonoids, and terpenoids. The varietal aroma of wines originated from the symbionts of the grapes and epiphytic microbiota, which combinedly grew from the sexual reproduction of the plant through seeds or during clonal reproduction. Nowadays, more and more studies focus on the assembly process of epiphytic microbiota and design a synthetic microbial community based on regional characteristics to improve wine quality and biofunctions. This review synthesizes the current concepts on the construction of synthetic microbiota, analyzes the advantages and difficulties in designing a functional yeast community, and lists the practical tools for data processing and model construction. A well-designed yeast community will possess high robustness against environment interference, higher efficiency of fermentation, and higher yield of targeted bioproducts.

1. Introduction

The application of commercial yeast in an incompatible fermentation scenario potentially causes the risk of stuck and sluggish fermentation and product homogenization. It even has a detrimental effect on the composition of regional microbiota [1], which play important roles in characterizing aroma profile and wine quality [2]. Currently, more and more researchers and wineries have realized the limitations of commercial yeasts and are gradually turning their attention to the study of non-Saccharomyces strains and the interactions between local microbiota and environmental factors [3, 4, 5]. Spontaneous fermentation has returned to the fore as a natural and regionally characterized winemaking mode [6]. Due to microbial succession, spontaneous fermentation usually presents higher sensorial complexity [4]. Although spontaneous fermentation can largely strengthen the ‘terroir’ characteristics of wines, the contamination by undesired microbes is the first issue and easily cause unpleasant flavors, such as a musty taste and mousiness [7, 8]. Besides, the high dimension of ecology of spontaneous fermentation results in many difficulties in control of fermentation and stable quality of wine production. With the development of high-throughput sequencing technology, researchers can conduct parallel analyses of the epiphytic microorganisms from vineyards and wines. It is effective to acquire an overview of microbial structure and position the pathogenic factors [2, 9, 10]. Therefore, the prerequisite for a well-controlled spontaneous fermentation is to identify the characteristic microbial resources in the production area, reduce the ecological dimension of the microeukaryotes community using self-assembled yeast communities, and establish a proper fermentation condition according to regional characteristics. In our previous study, two main environmental factors have been determined, namely, initial abiotic factors around vineyards and initial microbial community on the grape skin. Through the ‘top-down’ approach, we specially designed a self-assembled non-Saccharomyces yeasts group, which was composed of six selected yeast strains and highly expresses carbohydrate-active enzymes [2, 5]. The yeast group was first used to obtain a predictable and controllable fermentation mode during the spontaneous fermentation of icewine.

Icewine is a type of dessert wine produced from grape juice extracted from frozen grapes [11]. It is normally characterized by a fruity and flowery aroma and smooth taste. Icewine making is heavily dependent on the local environmental abiotic factors during late harvest. Fully ripe grapes hang on the vines for several months to suffer natural freeze-thaw cycles and be desiccated for the concentration of sugar (>35°Brix) [12]. This period has been proven to be effective for the accumulation of varietal aroma in grapes [2, 5]. For instance, high-odor-active compounds, such as terpenes and phenylalanine-derivates, can be largely developed during late harvest through a series of interactions between epiphytic microeukaryotes and grapes [2, 5]. Our previous studies have revealed the freeze-thaw cycles are the inevitable abiotic factors that cause cellular degradation and compartmentation of grape skin [2] and improve the aroma complexity of grapes during late harvest [5]. Therefore, selecting an icewine region is the essential prerequisite for producing high-quality icewine, which must present the regional characteristics [13]. Huanren area (Liaoning Province, Northeast China) is a representative icewine-producing area in China. It is characterized by a year-round cold climate, fewer problems of pests and disease [2]. Unique ecology shapes distinctive and stress-tolerant microeukaryotic communities, which remain active in pressed grape juice and become determinant to icewine fermentation [2, 13]. However, the mechanism and driving factors of the microeukaryotic assembly process are still unclear during icewine fermentation.

Metabolomics mainly studies the variety, quantity, and change rule of the metabolites with molecular weight less than 1500 Da caused by external stimulation, pathophysiological changes, and gene mutation. It is an extension of transcriptomics and proteomics, which accurately reflect the physiological state of organisms. Therefore, it effectively reveals the biological processes of biomarkers, the mechanism of the biological activities, and the regulatory pathways [14]. In terms of detection modes, metabolomics is mainly divided into untargeted analysis and targeted analysis. The untargeted metabolomics is usually based on a high-resolution mass spectrometer (triple TOF or QE) [15]. It can perform unbiased, large-scale, and systematic detection of various metabolites in experimental samples, providing an ‘aerial photography’ perspective to reflect the metabolic disturbance in the plants and microorganisms [16].

This study jointly used untargeted metabolomics and high-throughput sequencing technology to verify the effects of the self-assembled non-Saccharomyces yeasts group on the development of metabolites during AF of icewine. Meanwhile, we also studied two types of spontaneous fermentation using the extracted microbiota from the ripe Vidal blanc grape and the frozen grape after late harvest, respectively (see Figure 1). Our work first elucidated the mechanism and driving factors of microeukaryotic assembly in the different fermentation scenarios. It provides a reference to artificially design functional microeukaryotic communities and to improve the stability of spontaneous fermentation.

Figure 1.

Four inoculation methods were carried out for icewine fermentation. CS and HS were inoculated using the epiphytic microorganisms from the fully ripe grape berries and the grape berries in late harvest, respectively. HN was inoculated using the six core epiphytic yeast strains, which were important microeukaryotes in the grape skin. Control was conducted in icewine fermentation using a commercial yeast strain (BV818). Untargeted metabolomics and high-throughput sequencing technology were jointly used to analyze the self-assembly process of micro eukaryotic communities and the driving factors (important metabolites) during icewine fermentation.

2. Materials and methods

2.1 Vinification of V. blanc icewine

Spontaneous fermentation using the epiphytic microorganisms from the fully ripe grape berries (CS): The ripe V. blanc grape berries were manually collected in October 2021 and preserved at −20°C. To collect epidermal microorganisms, 200 g of frozen grape berries were soaked in 500 mL sterilized 10× PBS solution for 20 min and vortexed to elute epiphytic microorganisms [2]. The 50 mL suspension was centrifuged (10,000 rpm, 10 min) to remove the supernatant. Sediment was collected and first inoculated into 100 mL filtered grape juice for cultivation (25°C, 24 h). Subsequently, the cultivation medium was inoculated into 2 L filtered grape juice to carry out AF. The fermentation was supplied with 0.4 g/L nutrients (Nutrient-Vit, Lallemand, Montreal, PQ , Canada) and maintained the temperature at 18°C. At the beginning of AF, pumpover was carried out twice a day to facilitate yeast growth [4].

Spontaneous fermentation using the epiphytic microorganisms from the grape berries in late harvest (HS): based on previous research, V. blanc grape at the end of December can perform the optimal aroma characteristics [5]. Herein, the V. blanc grape berries were manually collected at the end of December 2021, and the epidermal microorganisms were inoculated into the filtrated grape juice. The specific procedures of inoculation and fermentation can refer to that of CS.

Icewine fermentation using the self-assembled non-Saccharomyces yeasts group (HN): based on the previous study, six core epiphytic yeast strains were separated from V. blanc grape skin [2]. These yeast strains play an important role in the self-assembly of epiphytic microeukaryotes and highly contribute to the development of the varietal aroma of V. blanc during late harvest. The percentage of each strain was calculated according to the importance of the random forest. The inoculation radio of each yeast strain is shown in Table 1. Each strain was cultured using a YPD medium to ensure yeast cell concentration higher than 10⁸ CFU/mL before inoculation. Subsequently, each strain was respectively inoculated into 2 L filtered grape juice according to the inoculation radio.

No.	Yeast species	Importance	Initial ratio (%)	Inoculation (mL)
1	Pichia kudriavzevii	88.05	21.7	8.7
2	Kluyveromyces lactis	80.59	19.8	7.9
3	Lachancea thermotolerans	77.00	18.9	7.6
4	Issatchenkia orientalis	63.69	15.7	6.3
5	Torulaspora delbrueckii	57.59	14.2	5.7
6	Candida dubliniensis	39.62	9.7	3.9

Table 1.

The percentage of each non-Saccharomyces yeast in artificial self-assembled yeast community.

Icewine fermentation using commercial yeast as control (BV818): a commercial yeast BV818 (AngelYeast Co., Ltd., Yichang) was inoculated into the filtrated grape juice for starting AF as control. This yeast belongs to Saccharomyces bayanus, which is stress-tolerant and highly active in a hyperosmotic juice. The fermentation was supplied with 0.4 g/L nutrients (Nutrient-Vit, Lallemand, Montreal, PQ , Canada) and maintained the temperature at 18°C. No sulfur dioxide or other antimicrobial agent was used in these four fermentation scenarios. All fermentations were conducted in triplicate.

2.2 Viticulture and grape juice filtration

V. blanc grapes were manually collected in the core icewine-producing area of Northeast China at the end of December 2021. The grape vines were planted at a commercial vineyard in Sidaolingzi village adjacent to Huanlong Lake, Huanren County, Benxi City, Liaoning Province, China (latitude 41°30′N, elevation 125°23′E) on a 25° slope. A vertical trellis system was installed in the vineyard with east-west oriented rows. The vines were 9 years old with a plant density of 2.5 m × 1 m. The vines were pruned to five to six shoots, and the overall crop yield was maintained at approximately 10,000 kg/ha [11]. Grape berries were destemmed and pressed to collect grape juice. Subsequently, the grape juice was carried through centrifugation (Kc > 6000) and cross-flow filtration (0.2 μm) for the removal of particles and indigenous microorganisms.

2.3 Determination of physicochemical parameters of icewine

Basic wine physicochemical parameters, such as hue, color intensity, Brix°, total sugar, total acid, glucose/fructose, lactic acid, malic acid, acetic acid, and glycerol, were determined using Y15 enzymatic autoanalyzer (Biosystems S.A., Barcelona, Spain). These analyses were performed using the appropriate enzymatic reaction kits purchased from Biosystems. Before detection, the Y-15 was calibrated with external standards that were technically supported by the Biosystems enterprise (www.biosystems.es). The pH value was measured using a PB-10 pH meter (Sartorius, Göttingen, Germany). The alcoholic degree of icewine was determined based on the National Standards of P. R. China (GB/T 15038-2006) [2, 4].

2.4 Untargeted metabolomics analysis

2.4.1 Sample collection and preparation

Icewine samples were centrifuged (4°C, 10,000 rpm, 10 min). The supernatant was added into precooled methanol/acetonitrile/water solution (2:2:1, v/v). The mixture was carried out vortex, ultrasonic homogenization (4°C, 30 min), stewing (−20°C, 10 min), and centrifugation (4°C, 12,000 rpm, 20 min). Then, the supernatant was conducted vacuum freeze-drying. Before metabolomics analysis, 100 μL acetonitrile aqueous solution (acetonitrile:water = 1:1, v/v) was added to redissolve the dry sample. Subsequently, the solution was vortexed and centrifuged (4°C, 12,000 rpm, 15 min) for UHPLC-MS analysis. Twenty microliters of each sample were mixed as a QC sample [17].

2.4.2 UHPLC conditions

Analyses were performed using a UHPLC (1290 Infinity LC, Agilent Technologies) coupled to a quadrupole time-of-flight (AB Sciex Triple TOF 6600). For HILIC separation, samples were analyzed using a 2.1 mm × 100 mm ACQUIY UPLC BEH 1.7 μm column (waters, Ireland). In both ESI positive and negative modes, the mobile phase contained A = 25 mM ammonium acetate and 25 mM ammonium hydroxide in water, and B = acetonitrile. The gradient was 85% B for 1 min, linearly reduced to 65% in 11 min, then reduced to 40% in 0.1 min, kept for 4 min, and then increased to 85% in 0.1 min, with a 5 min re-equilibration period employed.

2.4.3 ESI-Q-TOF MS/MS conditions

The ESI source conditions were set as follows: Ion Source Gas1 (Gas1) as 60, Ion Source Gas2 (Gas2) as 60, curtain gas (CUR) as 30, source temperature: 600°C, IonSpray Voltage Floating (ISVF) ± 5500 V. In MS-only acquisition, the instrument was set to acquire over the m/z range 60–1000 Da, and the accumulation time for TOF MS scan was set at 0.20 s/spectra. In auto MS/MS acquisition, the instrument was set to acquire over the m/z range 25–1000 Da, and the retention time for product ion scan was set at 0.05 s/spectra. The product ion scan is acquired using information-dependent acquisition (IDA) with a high sensitivity mode selected. The parameters were set as follows: the collision energy (CE) was fixed at 35 V with ±15 eV; declustering potential (DP), 60 V (+) and −60 V (−); exclude isotopes within 4 Da, candidate ions to monitor per cycle: 10.

2.4.4 Data processing

The raw MS data (wiff.scan files) were converted to MzXML files using ProteoWizard MSConvert before importing them into freely available XCMS software. CAMERA (Collection of Algorithms of Metabolite Profile Annotation) was used for the annotation of isotopes and adducts [18]. In the extracted ion features, only the variables having more than 50% of the nonzero measurement values in at least one group were kept. Compound identification of metabolites was performed by comparing the accuracy m/z value (<10 ppm) and MS/MS spectra with an in-house database established with available authentic standards [19].

2.5 Metabarcoding of internal transcribed spacer (ITS) sequence

The quantity and quality of extracted DNA were measured using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively. To analyze the taxonomic composition of the microeukaryotes communities, ITS1 of microeukaryotes 18S rRNA genes via a two-step amplification procedure using primers ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′), and ITS1R (5′-GCTGCGTTCTTCATCGATGC-3′) were carried out. Specific DNA extraction, PCR, and Illumina MiSeq sequencing (2- by 150-bp reads) were performed using the Illumina MiSeq platform at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China) [4]. Each sample was extracted for use in three replicates, and each extraction was analyzed twice. The Quantitative Insights into Microbial Ecology (QIIME2) pipeline was employed to process the sequencing data. Amplicon sequence variants were obtained using the procedure of denoise, quality control, splicing, and mosaicking, which were carried out using denoise-paired QIIME DADA2.

2.6 Statistical analysis

Statistics for the characteristic metabolites in V. blanc icewine, such as principle component analysis, partial least squares discriminant analysis (PLS-DA), and heatmap cluster with Euclidean distance, were carried out using the MetaboAnalyst 5.0 platform. ANOVA was conducted using Duncan’s multiple-range test (p < 0.05). Additionally, data visualization of α-diversity and β-diversity, spearman correlation (p < 0.05) of characteristic metabolites and microeukaryotes, volcano plots, principal coordinates analysis (PCoA), and random forest algorithm were carried out using relative packages in R version 4.2.2. The icewine metabolomics was mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to obtain the metabolic pathways in AF.

3. Results

3.1 Identification of the characteristic metabolites

The results of PCA and heatmap cluster of metabolomic illustrated different metabolite profiles between the four fermentation scenarios and V. blanc grape juice (Figure 2A). The metabolomics result was divided into two clusters. Organic acids and derivatives, benzenoids, hydrocarbon derivatives, and nucleosides were promoted during AF and mainly grouped in cluster 2. Their concentration increases, resulting from the growth of microorganisms. The mean of normalized scores in spontaneous fermentations (CS and HS) showed the highest positive values. In other words, spontaneous fermentation developed more organic compounds than the external yeast inoculation fermentations (HN and BV818). The weight of each metabolic compound represented its importance in the fermentation [20]. Obviously, fermentations generated abundant key metabolites and showed different metabolite profiles from that of grape juice. The important terpenes and flavonoids were grouped and marked on the right side of the heatmap. Some of them are first identified in icewine and have specific medical functions such as resveratroloside, neobavaisoflavone, narirutin, homoplantaginin, chrysosplenetin, and deoxysappanone b 7,3′-dimethyl ether acetate, etc. Besides, many high-valued compounds enriched in the grape juice were inhibited after AF, such as nobiletin, tectoridin, robinin, quercetin, Icariin, ginkgetin, etc.

Figure 2.

Comparison of characteristic metabolites in the icewine fermentation and grape juice. (A) Heatmap cluster of all identified metabolites in the four icewine fermentation and grape juice, especially the high-valued terpenes and flavonoids, were calculated and highlighted on the left side of the heatmap diagram. The importance of each treatment was shown on the top of the heatmap diagram. Herein, the Euclidean method was used for distance measurement, and the clustering algorithm was Ward. (B) Score plot of PLS-DA for analysis of observation of four icewine fermentations. (C) Ranking of characteristic metabolites by weighting the sum of VIP scores in compound 1 of PLS-DA. Colored boxes on the right side indicated the correlations between each metabolite and fermentation. Note: CS: spontaneous fermentation using the epiphytic microorganisms from the fully ripe grape berries; HS: spontaneous fermentation using the epiphytic microorganisms from the grape berries in late harvest; HN: icewine fermentation using the self-assembled non-Saccharomyces yeasts group icewine fermentation using a commercial yeast strain, which belongs to Saccharomyces bayanus.

Partial least squares discrimination analysis (PLS-DA) was carried out to compare the differences in each fermentation scenario (Figure 2B). Component 1 and component 2 contributed 76.6% of PLS-DA in total. It indicated the metabolite profile of each fermentation can be discriminated with the first two components. The four fermentations showed a far distance from the grape juice. Herein, the spontaneous fermentations (HS and CS) were close to HN, while BV818 showed a different metabolite profile from the others. To study the characteristic metabolites and their correlations to each fermentation. The top 20 characteristic metabolites were ranked according to their VIP scores of PLS-DA (Figure 2C). Most characteristic metabolites in icewine were positively correlated to the spontaneous fermentation and HN except glycerophosphocholine, succinate, and neohesperidose. In terms of spontaneous fermentation, some characteristic metabolites commonly showed positive correlations to CS and HS, such as d-gluconate, trans-3′-hydroxycotinine o-β-d-glucuronide, d-galacturonic acid, 2-isopropylmalic acid, gramine, adenine, and two dipeptides (Leu-His and Thr-Leu).

3.2 Spatial dispersal of microeukaryotic communities and correlations to the characteristic metabolites during AF

The spatial dispersal and co-occurrence network of each fermentation scenario jointly showed the dynamic changes of microeukaryotes during AF (Figure 3). Accordingly, the microeukaryotes in BV818 were highly developed at the end stage of AF, where ASVs of the microeukaryotes showed a higher enrichment but lower diffusivity (Figure 3A). The inner interactions between microeukaryotes in each module were positive. The major microeukaryotic modules in BV818 were modules 1, 2, 9, 7, and 14. Herein, the cross interactions between modules 1, 2, 7, and 9 showed higher negative similarity. It indicated the major activities between these modules were dominated by negative cohesion, such as competition, parasitism, or amensalism [21]. Moreover, module 14 showed fewer interactions with other modules. It indicated the related microeukaryotic communities could independently develop and be neutralized to other species during AF. HN showed a contrary spatial dispersal but a similar co-occurrence network to BV818 (Figure 3C). The strongest diffusivity of microeukaryotes was observed at the beginning of AF, while the weakest diffusivity was presented at the end of AF. The major microeukaryotic modules of HN were 1, 2, 6, and 12. The cross interactions between modules 1, 2, and 6 were also dominated by negative cohesion; meanwhile, these modules showed fewer interactions with module 12. Both CS and HS belonging to spontaneous fermentation and the same microbial origination from a vineyard resulted in similar spatial dispersals and co-occurrence networks of microeukaryotes during AF (Figure 3B and D). ASVs of microeukaryotes were evenly distributed in all the stages of AF. Comparatively, the microeukaryotes of CS mainly developed at the beginning and end stages of AF, while the microeukaryotes’ development of HS was vigorous at the middle stage of AF. The spontaneous fermentation presented more modules and more complex interactions than BV818 and HN. The major microeukaryotic modules in CS were 1, 2, 6, 8, 10, and 15. Herein, module 1 and module 10 showed positive cross interactions, namely, cooperation and mutualism. However, other microeukaryotic modules were mainly controlled by negative cohesion. The major modules in HS were 1, 2, 6, 7, 9, and 10. Herein, module 1 and module 7 showed positive cross interactions, while other modules were controlled by negative cohesion. In general, the cross-interaction between microeukaryotic phyla during AF was dominated by negative cohesion, and the abundance patterns of microeukaryotic communities depend on different fermentation modes [22].

Figure 3.

(A) Spatial dispersal and co-occurrence networks of fungi genera of BV818 on three stages during icewine fermentation. (B) Spatial dispersal and co-occurrence networks of fungi genera of CS on three stages during icewine fermentation. (C) Spatial dispersal and co-occurrence networks of fungi genera of HN on three stages during icewine fermentation. (D) Spatial dispersal and co-occurrence networks of fungi genera of HS on three stages during icewine fermentation. Each node represents one ASV that is labeled by genus. A node was verified by a robust (Spearman’s correction coefficient R > 0.6) and significant (pFDR < 0.05) correlation. The size of each node is relational to the number of connections, while nodes with the same color display the same module. The thickness of each connection between two nodes is relational to the strength of Spearman’s correlation coefficient. Note: CS: spontaneous fermentation using the epiphytic microorganisms from the fully ripe grape berries; HS: spontaneous fermentation using the epiphytic microorganisms from the grape berries in late harvest; HN: icewine fermentation using the self-assembled non-Saccharomyces yeasts group, icewine fermentation using a commercial yeast strain, which belongs to Saccharomyces bayanus.

Mantel test was used to study the correlations between the specific order of microeukaryotes and fermentation mode (Figure 4A). As a result, Saccharomycetes, Microbotryomycetes, and Tremellomycetes showed strong positive correlations (p < 0.05), which potentially constructed a stable microeukaryotic community during AF and increase the efficiency of alcohol conversion. However, significant negative correlations were found between some orders, such as Sgaricomycetes, Lecanoromycetes, Pezizomycetes, Cystobasidiomycetes, Zoopagomycetes, Exobasidiomycetes, Sordariomycetes, Leotiomycetes, Wallemiomycetes, Microbotryomycetes, and Tremellomycetes. It probably results in interspecific competition and internal friction of resources during fermentation. Moreover, Ustilaginomycetes and Zoopagomycetes were highly correlated to BV818 and HN, respectively.

Figure 4.

(A) Correlations of the networked community structures (Bray–Curtis distance) with important fungi orders. The edge width corresponds to Mantel’s r value, and the edge color denotes the statistical significance. Significance levels were labeled as “*” (p-value < 0.1), “**” (p-value < 0.05), and “***” (p-value < 0.01), respectively. Pairwise correlations of these variables were shown with a color gradient denoting Pearson’s correlation coefficients. (B) The correlation matrix for the top 25 high-valued terpenes and flavonoids and the top 25 active fungi genus was established using a random forest algorithm, and the importance of the coefficient was shown with different cycles. Meanwhile, the explained ratio of the top 25 active fungi genera was represented on the right side of the diagram. Note: CS: spontaneous fermentation using the epiphytic microorganisms from the fully ripe grape berries; HS: spontaneous fermentation using the epiphytic microorganisms from the grape berries in late harvest; HN: icewine fermentation using the self-assembled non-Saccharomyces yeasts group, icewine fermentation using a commercial yeast strain, which belongs to Saccharomyces bayanus.

The correlation matrix between the functional compounds and the top 25 active microeukaryotic genera was established using a random forest model and Pearson’s correlation coefficients (Figure 4B). Accordingly, many microeukaryotes showed positive correlations. They potentially contributed to the development of terpenoids and flavonoids in icewines, such as Starmerella, Mrakia, Tausonia, Rhodotorula, Mrakiella, Verticillium, and Curvibasidium. Herein, Starmerella played an important role in promoting the metabolism of bioflavonoids and polyflavonoids, flavans, and flavones. Fatty alcohols were strongly correlated to many microeukaryotes, such as Mrakia, Rhodotorula, Cystofilobasidium, Mrakiella, Filobasidium, Bettsia, Verticillium, and Curvibasidium. Homoisoflavans positively correlated to Mrakia, Rhodotorula, Mrakiella, Verticillium, and Curvibasidium. Moreover, Cystofilobasidium and Golubevia contributed to the development of tryptamines and derivatives, while Pichia was regarded as an important factor in affecting the development of monoterpenoids and tryptamines.

3.3 Characteristic metabolites of V. blanc grape induced stochastic assembly of microeukaryotic communities during AF

One of the reasons why spontaneous fermentation could not be widely applied is that many winemakers treat it as a ‘black-box’, a complex, multi-strains participated system. With the development of systematic biology, the neutral community model (NCM) will be a useful tool to excavate the deep principle of microbial assembly in co-fermentation or spontaneous fermentation [23]. In this study, the NCM predicted a large fraction of the total correlation between the occurrence frequency of ASVs and their relative abundance variations (Figure 5A), with 83%, 73%, 80%, and 76% of the explained microeukaryotic community variations for the BV818, CS, HN, and HS, respectively. The NCM indicates that microeukaryotic communities comply with stochastic assembly instead of deterministic assembly in different fermentation scenarios. Higher R² suggests the sample is closer to the NCM. Compared to spontaneous fermentation, commercial yeast inoculation (BV818) fits the stochastic process. Moreover, a smaller Nm value suggests more restrictive species dispersal. Therefore, the level of species dispersal in BV818 (Nm = 374) was the lowest in all fermentations, while HS had the maximum species dispersal (Nm = 1005). CS and HN showed similar Nm values. This result is also in accord with the spatial dispersal of microeukaryotic communities during AF; namely, the microeukaryotic communities of HS can evenly develop at three stages of AF, while the ASVs of microeukaryotes of BV818 are mainly concentrated at the end of AF.

Figure 5.

(A) Fit of the neutral community models of fungi communities assembly in the four fermentations. The solid blue lines indicate the best fit to the NCM, and the dashed blue lines represent 95% confidence intervals around the model prediction. OTUs that occurred more or less frequently than predicted by the NCM are shown in different colors (gray dots represent the frequency of occurrence within the 95% confidence interval ranging around the model prediction, orange dots represent above prediction, and purple dots represent below prediction). Nm indicates the metacommunity size times immigration, R² indicates the fit to this model. Meanwhile, a modified stochasticity ratio was used to assess the relative importance of stochasticity. (B) Coefficients between |β-NTI| and the relative abundance of the six featured metabolites were calculated to evaluate the correlative fitting curves, p-values were evaluated, and significant correlations were determined at p < 0.05. (C) Contributions of deterministic and stochastic processes on community assembly within four icewine fermentations. β-NTI calculations of phylogenetic turnover between different icewine fermentations indicate that variable selection has greater effects on fungi dynamics and fermentation progress. (D) A modified stochasticity ratio was carried out to assess the relative importance of stochasticity, which was developed with 50% as the boundary line in each icewine fermentation. Note: CS: spontaneous fermentation using the epiphytic microorganisms from the fully ripe grape berries; HS: spontaneous fermentation using the epiphytic microorganisms from the grape berries in late harvest; HN: icewine fermentation using the self-assembled non-Saccharomyces yeasts group, icewine fermentation using a commercial yeast strain, which belongs to Saccharomyces bayanus.

β-Nearest taxon index (β-NTI) between samples was calculated to describe the phylogenetic turnover of microeukaryotic communities in the different fermentation scenarios. Linear regression curves between |β-NTI| and the relative abundance of the six characteristic metabolites were constructed (Figure 5B). Among these six characteristic metabolites, β-sitosterol is a plant steroid, empenthrin, gentiopicroside, zerumbone, and qingyangshengenin belong to terpenoids, and phlorizin is a kind of flavonoids. Empenthrin, β-sitosterol, gentiopicroside, phlorizin, and zerumbone positively correlated to |β-NTI| while gentiopicroside and qingyangshengenin showed negative correlations to |β-NTI|. However, β-NTI of all the fermentations were in the range of −2 to 2 (Figure 5C). Although gentiopicroside and qingyangshengenin could affect the stochasticity of microeukaryotic communities, the stochastic process still dominated the phylogenetic turnover of microeukaryotic communities during AF. To assess the relative importance of determinism and stochasticity, modified stochasticity ratio (MST) was carried out to supply the results of NCM (Figure 5D). In four fermentation scenarios, MST of CS (36.1%) was much lower than the boundary (50%), which indicated determinism affected the microeukaryotic assembly of CS. The MST of HS (54.7%) was slightly higher than the MST of HN (53.0%) without a significant difference (p < 0.05). As expected, BV818 showed the highest MST (63.8%), which indicated stochasticity played an important role in shaping microeukaryotic communities of BV818. This result suggested the assembly process of spontaneous fermentation could be more stable than the inoculation of commercial yeasts as long as the source of grape is consistent and fermentative conditions are standard.

3.4 The specific process of stochastic assembly of microeukaryotic communities and the selection of fermentative specialists

Principal coordinates analysis (PCoA) was carried out to observe the degree of similarity between the microeukaryotic communities of different fermentation scenarios (Figure 6A). The total contribution of the first two components of PCoA can explain 66.11% of variations. The confidence intervals of HN and BV818 have been separated from the spontaneous fermentations (CS and HS). Comparatively, the degree of distribution of BV818 was more concentrated than that of the others while HN and HS showed higher degree of distributions than CS and BV818. The result of PCoA indicated that the microeukaryotic composition of BV818 was different from co-fermentations (HN, CS, and HS), and the concentrated distribution of BV818 also indicated low diversity of microbial participation during AF. Moreover, the top five abundant microeukaryotic genera were selected and located on the PCoA plot. Herein, Pichia (ASV_382) was located in the confidence ellipse of HN (p < 0.05), while Vishniacozyma and Starmerella (ASV_1133, ASV_496, ASV_327, and ASV_446) highly related to CS and HS.

Figure 6.

(A) Principal coordinate analysis (PCoA) based on Sørensen dissimilarly displaying differences in community composition of fungal genera in the four icewine fermentations, the top 10 highest active fungi genera were selected and labeled. (B) The relative influence of each community assembly process among four icewine fermentation was defined by the percentage of site pairs governed by each process. (C) Partitioning the turnover and nestedness components of β-diversity of fungi communities in four icewine fermentations. Duncan’s multiple range test was carried out for the significant differences of the model at p < 0.05. (D) The SPEC-OCCU plots showed the most abundant OTUs in each icewine fermentation; the x-axis represents occupancy, i.e., how well an OTU is distributed across the specific fermentation scenario; and the y-axis represents specificity, i.e., whether they are also found in other fermentation scenarios. The fermentative specialists of fungi species in each fermentation were selected as both occupancy and specificity higher than 0.7.

Except for the NCM, phylogenetic-bin-based null model analysis was carried out further to infer community assembly mechanisms in the stochastic process (Figure 6B) [20]. The important assembly processes of microeukaryotic communities were heterogeneous selection, homogeneous selection, dispersal limitation, homogenizing dispersal, and drift. Specifically, the assembly process was mainly composed of homogeneous selection (43.1%) and drift (56.2%) in BV818. For HN, dispersal limitation (0.05%), homogeneous selection (45.1%), and drift (49.7%) were the important parts of the stochastic assembly. For CS, the specific stochastic assembly included homogeneous selection (25.5%), dispersal limitation (0.05%), homogenizing dispersal (0.04%), and drift (65.4%). The heterogeneous selection was the characteristic assembly process in HS, which possessed 0.05% of the stochastic assembly. Besides, the proportion of homogeneous selection, dispersal limitation, homogenizing dispersal, and drift was 15.7%, 10.5%, 0.03%, and 65.4% of the stochastic assembly in HS, respectively. In general, drift and homogeneous selection were the major assembly processes during AF of icewine, and HS showed a more complex composition of the assembly process than other fermentations. Partitioning β-diversity of richness was carried out to quantify the result of species replacement between different fermentations (turnover) and species gains or losses between different microeukaryotic communities (nestedness) (Figure 6C). As a result, species replacement played an important role in shaping microeukaryotic communities. The greatest difference between the contributions of β-diversity components was observed in the four fermentation scenarios, where the turnover was about six times higher than nestedness. Specifically, the range of turnover in each fermentation was from 50% to 55% while the range of nestedness was from 5% to 8%. The nestedness of HN was the highest level among the four fermentations. For turnover, BV818 and CS were higher than HN and HS.

To inspect the distribution of OTUs from each fermentation and the specificity of these microeukaryotes during AF, specificity and occupancy were calculated for each OUT, which was then projected onto a plot (SPEC-OCCU plot, Figure 6D). As indicated by the spread of OTUs across occupancy, OTUs from BV818 and HN communities showed highly varied occupancy while the majority of OTUs from CS and HS exhibited more homogenous occupancy, which mainly concentrated at the range of 0–0.25. To find specialist species attributable to each fermentation, we selected OTUs with specificity and occupancy greater or equal to 0.7 (dotted boxes). These microeukaryotes are specific and could be highly fermentative in their fermentation scenarios. The number of these specialist OTUs was BV818 (3 OTUs represent), HN (6 OTUs), CS (2 OTUs), and HS (3 OTUs), respectively. Specifically, Mrakia frigida is the fermentative specialist in BV818, the specialists of HN included Lachancea thermotolerans and Pichia kudriavzevii, while Starmerella bacillaris was the fermentative specialist in spontaneous fermentation. Meanwhile, the changes in the relative abundance of microeukaryotic species verified these yeasts are the key factors to influence fermentation process.

3.5 The analysis of metabolic pathways in four icewine fermentation scenarios

A structural equation modeling (SEM) accounting for the yeast strains, characteristic metabolites, and key assembly factors was built to analyze the interactions of variables in four icewine fermentation scenarios (Figure 7A). Based on the previous results, five yeast strains (Vishniacozyma victoriae, Lachancea thermotolerans, Pichia kudriavzevii, Mrakia frigida, and Starmerella bacillaris) were assigned as important components of microeukaryotic community in the SEM. Average variation degree (AVD), alpha diversity, MST, and bNTI were assigned as key factors of the assembly process of microeukaryotic community. Terpenoids, flavonoids, and steroids were assigned as characteristic metabolites of icewine fermentations. As a result, M. frigida showed a strong connection to yeast combination, which also significantly related to microeukaryotes assembly process during icewine fermentation and except for P. kudriavzevii, other four yeast strains, V. victoriae, L. thermotolerans, S. bacillaris, and M. frigida negatively affected the microeukaryotes assembly process. Moreover, the characteristic metabolites showed negative effects on the microeukaryote assembly process. Steroids showed a strong and positive influence on the total metabolites of icewine. AVD showed an indirect effect on the microeukaryotes assembly process, while AD, MST, and bNTI directly and positively related to the microeukaryotes assembly process during alcoholic fermentation.

Figure 7.

(A) Structural equation modeling (SEM) accounting for the yeast strains, characteristics, metabolites, and key assembly factors in four icewine fermentation scenarios. (B) Top 20 enriched KEGG pathways of the fungal community in four icewine fermentation scenarios. Numbers adjacent to measured variables are their coefficients with composite variables. Numbers adjacent to arrows are path coefficients are the directly standardized effect size of the relationship. The thickness of the arrow represents the strength of the relationship. Note: CS: spontaneous fermentation using the epiphytic microorganisms from the fully ripe grape berries; HS: spontaneous fermentation using the epiphytic microorganisms from the grape berries in late harvest; HN: icewine fermentation using the self-assembled non-Saccharomyces yeasts group, icewine fermentation using a commercial yeast strain, which belongs to Saccharomyces bayanus.

The top 20 enriched KEGG pathways of four icewine fermentations were shown in Figure 7B, where the rich factors of each KEGG pathway were the ratio of differential genes in this pathway, while the p-value indicated the importance of this pathway in the icewine fermentation. Accordingly, protein digestion, absorption, and ABC transporters were the most important pathways of microeukaryotic metabolism in all fermentation treatments. Compared to commercial yeast, CS, HN, and HS showed similar categories of important KEGG pathways, which also contained mineral absorption, biosynthesis of amino acids, and aminoacyl-tRNA biosynthesis. Lysine degradation was the characteristic pathway for the spontaneous fermentations (CS and HS). Starch and sucrose metabolism were relatively important to HS and BV818, while proximal tubule bicarbonate reclamation was the common pathway in BV818 and CS. Interestingly, valine, leucin, and isoleucine biosynthesis and nicotinate and nicotinamide metabolism were the specific pathways in CS. It indicated more potential volatile metabolites could be synthesized by applying the indigenous yeast community during alcoholic fermentation.

4. Discussion

In this study, the top-down approach was first used to design a self-assembled non-Saccharomyces yeasts group (HN) [24]. Herein, ‘top’ referred to the special ecosystem in which the local epiphytic yeasts of the V. blanc grape performed a strong ability to hydrolyze the glycosidic bound terpenes and enhanced the varietal aroma of the grape during late harvest [2, 21]. We tested the availability of HN to conduct AF and explored its biological functions of producing high-valued metabolites in icewine. Besides, two types of spontaneous fermentations were carried out to observe the assembly process of indigenous microeukaryotes during AF. The microbial starter of the spontaneous fermentation was extracted from the fully ripe V. blanc grape and the grape berries in the late harvest. Therefore, the spontaneous fermentations simultaneously represented the regional microeukaryotic composition and the “terroir” characteristics of icewine. As a control, BV818 showed a different metabolite profile from the other fermentations according to PLS-DA. It indicated the application of commercial yeasts may affect the indigenous microbial communities and change wine regional characteristics [25, 26]. Terpenoids and flavonoids were identified as the characteristic metabolites in icewine. Some high-valued metabolites of icewine were first reported in our work. For instance, resveratroloside is a monoglucosylated form of stilbene in red wine. It effectively alleviates postprandial hyperglycemia and cardioprotective effects [27, 28]. Neobavaisoflavone has the potential for anti-human glioma cancer and inhibition of inflammatory mediators. Deoxysappanone b 7,3′-dimethyl ether acetate has been reported to have anti-cryptosporidial activities [29]. Based on PLS-DA, the spontaneous fermentation, and self-assembled non-Saccharomyces yeasts group had the advantage of producing more characteristic metabolites.

The spatial dispersal and co-occurrence network of each fermentation scenario performed a dynamic distribution of microeukaryotes and the correlations between microeukaryotic communities during AF. As a result, negative cohesion, such as competition, parasitism, or amensalism, played an important role in the cross-interactions between microeukaryotic communities. The microeukaryotic dispersal of BV818 was mainly concentrated at the end of fermentation, while HN showed the opposite result that the microeukaryotes mainly distributed at the beginning of fermentation. Understandably, most commercial yeasts were Saccharomyces cerevisiae, which usually performs killer activity against other fungi at the beginning of AF, then dominates and completes the fermentation. Although the spontaneous fermentation obtained quite a balanced development of microeukaryotic communities at the three stages of AF, fermentation stress and frequent microbial succession jointly resulted in more complex interactions between microeukaryotic communities in the spontaneous fermentations [22, 30]. Moreover, many microeukaryotic genera were positively related to the development of terpenoids and flavonoids in icewine. It is due to these compounds commonly present in glycosidic bound forms, while the epiphytic microeukaryotes of grapes can generate a series of glycoside hydrolases (GHs) and glycosyl transferases (GTs), which efficiently break glycosidic bonds and release the aglycones [2, 31]. For instance, monoterpenyl disaccharide glycosides originated from plants and microbes, such as α-arabinofuranosidase, β-glycosidase, and β-galactosidase, hydrolyze glycosidic bound terpenes into free forms [32]. The result of NCM indicated the microeukaryotic community assembly in icewine was predominately driven by births, deaths, and immigration (stochastic process) instead of external abiotic factors and niche differentiation (deterministic process) [23]. Therefore, the hydrolytic/fermentative microeukaryotes in AF are inclined to spontaneously conduct stochastic assembly through the response or sharing of certain organic compounds [33]. Based on the result of random forest ranking of metabolites and β-NTI of ASVs, β-sitosterol, four terpenoids (empenthrin, gentiopicroside, zerumbone, and qingyangshengenin), and phlorizin were identified as the characteristic metabolites to affect the stochastic assembly of microeukaryotic communities. These metabolites may induce the common response of microeukaryotes to biotic and abiotic stresses and contribute to the assembly process of microeukaryotic communities. For instance, β-sitosterol could be used as an antibiotic to enhance the antibacterial activity of microeukaryotes [34]. Gentiopicroside is significantly correlated with the abundance of dominant endophyte genera among different tissues of Gentiana officinalis [35]. Moreover, fungal life history strategies determined homogenous selection and dispersal limitation attributed to the large proportion in structuring microeukaryotic communities, e.g., growth, dispersal, or dormancy methods [36].

People had been drinking natural wines for thousands of years before sulfur dioxide was applied in the wine industry. Current studies have admitted the importance of a certain strain of commercial yeast in the winemaking industry, while indigenous fungi have been misunderstood as microbial contaminators during AF. In our study, the dominator of AF was indigenous non-Saccharomyces yeasts instead of commercial yeast. Some fermentative specialists not only complete AF but highly contribute to the stochastic assembly of microeukaryotic communities. In particular, S. bacillaris dominated the spontaneous fermentation and has been reported high glycerol production, moderate nitrogen demand, and more active adaptation to extremely high sugar concentration [37]. L. thermotolerans and P. kudriavzevii are multi-stress-tolerant and widely present in different kinds of fruits [38]. Appropriate use of these indigenous non-Saccharomyces yeasts can enhance the aroma complexity of icewine [4]. S. cerevisiae was difficult to be the advantage yeast strain since the growth was always inhibited throughout AF of icewine. Inoculation of BV818 induced the development of M. frigida, which is highly active in producing β-galactosidase and killer toxins [39]. The result of SEM also confirmed M. frigida could be a key defender to restrict the development of the invasive exotic species and ensure the stability of indigenous microeukaryotic communities. Our results elucidated the importance of indigenous non-Saccharomyces yeasts in stochastic assembly, while inoculation of commercial yeast will induce reconstruction of microeukaryotic communities and potentially affect the “terroir” characteristics of icewine. Nevertheless, steroids, especially β-sitosterol, are not only commonly upregulated in different fermentation scenarios but play an important role in driving the stochastic assembly of microeukaryotic communities during AF. Vitamin digestion, absorption, and lysine degradation are the characteristic enriched pathways in spontaneous fermentation. It mainly relates to the stress-resistant functions of microeukaryotes, such as antioxidation, low-temperature resistance, and salt resistance [40].

5. Conclusion

This study first explored the assembly process of microeukaryotic communities in four typical icewine fermentation scenarios. Herein, top-down design was carried out to construct a non-Saccharomyces yeasts group for icewinemaking. Untargeted metabolomics revealed terpenoids and flavonoids were the characteristic metabolites of icewine. Besides, some microeukaryotic families were strongly correlated to the development of the key classes of terpenoids and flavonoids. Through NCM and the β-diversity of microeukaryotes, a stochastic process dominated microeukaryotic community assembly during AF. Moreover, the spontaneous fermentation showed an evenly spatial distribution of microeukaryotes at the three stages of AF. Co-occurrence networks suggested that negative cohesion played a pivotal role in microeukaryotic interactions between different modules. Five terpenoids and a flavonoid showed their high contributions in driving stochastic assembly. Meanwhile, mechanisms of microeukaryotic community assembly are mainly involved in the homogeneous selection, drift, and dispersal limitation during AF of icewine. Indigenous non-Saccharomyces yeasts, M. frigida, L. thermotolerans, P. kudriavzevii, and S. bacillaris were identified as the fermentative specialists during AF of icewine. Furthermore, the enriched KEGG pathways of microeukaryotes related to their stress-resistant functions, which result in higher production of terpenoids, flavonoids, and steroids in the spontaneous fermentation.

Acknowledgments

We acknowledge Huanren Senpatina Icewine Domaine Co., Ltd. for kindly providing the experimental raw materials.