A first set of 47 concentrated maple sap from 11 sites harvested during the flow period in 2015, 2016, and 2017, and previously characterized in Filteau et al. (2017), N’guyen et al. (2018) were used in this study to explore sap bacterial and fungal communities. Those samples were collected on a voluntary basis from maple sap producers in Québec and New Brunswick. As is typical for maple syrup production, sap was collected from hundreds to thousands of trees, mainly Acer saccharum. Sap collection was achieved using a network of plastic tubing lines under vacuum, then concentrated by membrane processes. Sap concentrate samples were then collected right before the evaporation process and the corresponding syrup was sampled after the batch was processed. Samples were frozen on site at −20°C. Supplementary Data S1 summarizes the samples used in this study and their characteristics. The concentration (°Brix) was measured in laboratory with a digital refractometer and the average density of sap concentrate samples was 13.2 °Brix. To obtain the concentrated microbial biomass 50 mL aliquots of sap concentrate were rapidly thawed in a cold-water bath and then centrifugated for 5 min × 915 g at 4°C. The supernatant was partially removed to conserve 5 mL of the final solution and stored at −80°C until analysis. For these samples, we used the reported data for ammonium, organic sulfur, sulfate and phosphate quantification in sap (N’guyen et al., 2018). When available, we also used the corresponding syrup quality classification which was performed by trained analysts following the standard organoleptic assessment procedures officially used by the industry as described in N’guyen et al. (2018). Since the classification nomenclature is not informative for our purposes, we grouped samples by organoleptic conformity: standard (standard and allowing for light off-flavors (v), N = 11) and non-standard (1, 4, 5 and 6 defect types leading to product downgrading or rejection, N = 21).

For follow up experiments, a second set of 42 samples were collected in 2018 and 2019 from five other production sites following similar handling methods. Maple sap concentrated by membrane processing at the farm was frozen on site at −20°C until analysis. After thawing, saps were sterilized by 0.2 μM filtration before analysis.

Colorimetric analysis of available corresponding syrups (N = 37) using the CIELAB color space (L*a*b*) were performed in laboratory with a Spectrophotometer UV-2700 (Shimatzu, Kyoto, Japan) and the UVPC Color Analysis Software, version 3.10.

We previously introduced a dormancy release index representing the remaining Sum of cumulative temperature necessary to reach Bud Break (S_bb) as a temporal and meteorological reference to compare sap samples (N’guyen et al., 2018). It is important to note that this index is independent of the temperature occurring after the sap sample collection date, and thus allows for comparison between samples of different years. Meteorological records from the nearest stations were obtained from the Canadian government public records (http://climate.weather.gc.ca) and missing data were estimated using data from the nearest station available.

The DNA extraction was performed as in (Bleuven et al., 2020). A volume of 1 mL of sap concentrate was used for DNA extractions. Bacterial and fungal contamination levels were quantified by the number of 16S and ITS amplicon copy numbers by droplets digital PCR (ddPCR) with the QX200 droplet generator (Bio-Rad). The procedure was done according to the manufacturer recommendation (Taylor et al., 2015). Briefly, to adjust the DNA concentration of samples into the detection range of the ddPCR, a reference sample was made by pooling 2 µl of each sample. Then, 2 µl pooled sample was then amplified by RT-PCR in 7500 Real-Time PCR system (Life Technologies) with a mix made with 17 µl of PowerUp SYBR Green Master mix (Life Technologies) and 1 µl of 10 µM primers targeting the 16S V3-V4 region (Klindworth et al., 2013) or ITS2 region (Tedersoo et al. (2015). Each sample were diluted by a factor determined by adjusting the reference Ct value to 25Ct. After dilution, the ddPCR protocol (Taylor et al., 2015) was followed for each target, 16S or ITS. The PCR cycle was a first step at 95°C for 5 min followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 56°C for 30 s, extension at 60°C for 30 s and a final stabilisation step at 4°C for 5 min and at 90°C for 5 min. The total DNA copy number of 16S and ITS were then used to estimate the absolute abundance of bacteria and fungi in each sample, respectively.

To explore the microbial community diversity, we amplified and sequenced the 16S subunit gene with primers covering the V3-V4 region (Klindworth et al., 2013) and ITS2 sequence specific regions described in Tedersoo et al. (2015) for bacteria and fungi, respectively. Amplification was performed in a two-step dual-index PCR approach specifically designed for Illumina instruments by the Plateforme d’Analyses Génomiques (IBIS, Université Laval, Quebec City, Canada) as in Bleuven et al. (2020). Briefly, primer sequences were designed by the fusion of the gene specific sequence and the Illumina TruSeq sequencing primers and PCR was carried out in a total volume of 25 µl that contains 1X Q5 buffer (NEB), 1X high CG content buffer, 0.25 µM of each primer, 200 µM of each dNTPs, 1 U of Q5 High-Fidelity DNA polymerase (NEB) and 1 µl of template DNA. The PCR started with an initial denaturation at 98°C for 30 s followed by 35 cycles of denaturation at 98°C for 10 s, annealing at 55°C for 10 s, extension at 72°C for 30 s and a final extension at 72°C for 2 min. Quality of the purified PCR product were checked on a 1% agarose gel. Fifty-fold dilution of this purified product was used as a template for a second PCR step with the goal of adding barcodes (dual-indexed) and the specific sequence required for Illumina sequencing. Cycling for the second PCR were identical to the first PCR with 12 cycles. PCR reactions were checked for quality on a DNA7500 Bioanalyzer chip (Agilent). Barcoded Amplicons were then pooled in equimolar concentration based on agarose gel band intensity, purified with the Axygen PCR cleanup kit (Axygen) and sequenced on the illumina Miseq (paired–end 300 bases with two index reads) at the Plateforme d’analyses génomiques (IBIS, Université Laval, Quebec City, Canada). Please note that primers used in this work contain Illumina specific sequences protected by intellectual property (Oligonucleotide sequences ^© 2007–2013 Illumina, Inc. All rights reserved. Derivative works created by Illumina customers are authorized for use with Illumina instruments and products only. All other uses are strictly prohibited.) The 16S and ITS amplicons were sequenced in two independent Illumina runs.

The bioinformatics workflow applied for this study was similar to (Bleuven et al., 2020). Raw sequences obtained from the Illumina sequencing were first demultiplexed and sequences with low length were removed following the Illumina default parameters. For amplicon assignation to amplicon sequence variant (ASV), several steps were followed based on Dada2 package version 1.5.0 developed by Callahan et al. (2016) and workflow (3) in R version 3.1.2 (http://www.R-project.org). The resulting ASVs were obtained with the following parameters. For the filtration step and quality check, sequences were trimmed at 255 pb for forward reads and 225 pb for reverse reads. The first 10 base pairs were truncated for all reads. Ambiguous bases and more than two expected errors were filtered out for all reads. After filtration, the average of reads count for library was 31,538 (min: 4020; max: 75,505) for 16S and 75,641 (min: 33,291; max: 112,897) for ITS. An average of reads filtered out was 27% +/− 4 for 16S and 57% +/− 14 for ITS mostly due to low length sequences. Then, the amplicon error learning steps for the forward and reverse reads were performed on 100M nucleotides using the DADA2 algorithm with default parameters. Finally, the denoised output reads were generated by dereplication using the derepFastq function from the DADA2 package and all reads with any mismatches were removed. Then forward and reverse reads were merged, and chimera sequences were removed with the uchime_denovo function from the vsearch package with a minimum difference set at 2. Purged ASV sequences with more than 4 read counts were kept. For 16S corresponding ASV, 3928 unique ASV sequences were found with 2100 (53.5%) considered as chimeras representing 7.6%) in 1175149 total sequences. For ITS, 463 unique ASV sequences were found with 103 (40.6%) considered as chimeras representing 1.4% in 1965903 total sequences.

ASVs sequences were then annotated by BLAST matching the NCBI database, excluding the Uncultured/environmental sample sequences. From the top 50 hits, the best score was kept. If several species had a similar score, the sequence was renamed as “Genus.sp”. To exclude potential mislabeling, if several Genus has a similar score, genus represented in more that 1/5 of the best scores were concatenated and the ASV renamed as Genus.species_1-Genus.species_2. Finally, ASVs associated to the same species names were aligned and sequences presenting less that 1% difference were grouped into Operational Taxonomic Units (OTUs). This last step was performed to reduce false positives and to decrease sparsity in the abundance matrix. The relative abundance of each OTU based on sequence counts was converted to absolute abundance estimates using the 16S and ITS total amplicon copy number evaluated by ddPCR in each library. Abundance is presented as the log10 (amplicon copy numbers +1), so that null values can remain informative.

We used the web platform microbiome analyst to compare alpha and beta diversity between samples (Dhariwal et al., 2017). Absolute abundances were used as the input and Shannon diversity index was calculated with unfiltered data. Principal coordinates analysis (PCoA) ordination was performed on OTUs with at least 10% occurrence, using the Bray-Curtis distance. Homogeneity of group dispersions (PERMDISP) and analysis of group similarities (ANOSIM) were performed to compare S_bb and organoleptic conformity groups.

To model the relationships in our high-dimensional dataset, a partial least square approach to model fitting was selected. This approach is particularly useful when there are more explanatory variables than observations or when the explanatory variables are highly correlated, which is the case with our OTU data. We applied a nonlinear iterative partial least squares analysis (JMP14, SAS institute) to OTUs as predictor of sap composition (ammonium, organic sulfur, sulfate and phosphate) and syrup L*a*b* values. Because of missing values, separate models were fitted for each variable and the number of factors was cross validated using the Root Mean PRESS (predicted residual sum of squares) approach. OTUs that were detected in a minimum of five samples were included in the analysis and ranked by Variable Importance for Projection (VIP) to identify the top predictors for each parameter tested.

Further, we used bootstrap forest partitioning (JMP14) to evaluate the contribution of OTUs as predictors of maple syrup conformity (standard or non-standard). The bootstrap forest models (100 decision trees each) were built on multiple predictors; therefore, the method can identify OTUs that might be weak predictors alone but strong when used in combination with other OTUs. The analysis was run ten times on OTUs occurring in at least five samples and the contribution of each OTUs was ranked to identify the consistently best predictors.

The estimated abundance adjusted by ddPCR was used to generate a networks for the Early and Late S_bb groups using the CoNet algorithm (Faust and Raes, 2016) and SparCC (Sparse Correlations for Compositional data) (Friedman and Alm, 2012). OTUs with less than three occurrences were not considered and the Spearman correlation threshold was set at 0.6. Associations observed with both methods were used to construct the final networks. The network characteristics and figures were generated with Cytoscape v3.8.2.

We used the web-based application Piphillin (Narayan et al., 2020) to predict the metagenomic content of maple sap bacterial communities. For greater precision, we used as input the ddPCR-adjusted abundances of each ASV and selected a conservative 99% identity cutoff for genome annotation with the Biocyc 24.0 database. The annotations of 686 (out of 1598) ASV were inferred from 181 genomes yielding 5988 Biocyc inferred features. Spearman’s non-parametric correlations of each feature abundance with the S_bb index and FDR adjusted Pvalues were computed with JMP14.

Allantoin and allantoate were quantified as described in Filteau et al. (2017). Briefly, the maple sap concentration of solid content measured with a digital refractometer was adjusted to 10 °Brix. Then samples were diluted 4 × in 90:10 (acetonitrile:water) containing allantoin-13C2,15N4 as an internal standard (1ppm), vortexed for 1 min and centrifuged at 16,000 g for 5 min at 4°C. Samples were then filtered through a 0.22 µm nylon filter prior to injection in a Acquity H-Class Ultra-Performance LC system coupled to a Xevo TQD mass spectrometer (Waters, Milford, MA, United States).​

Ammonium and sulfate were quantified as described in N’guyen et al. (2018), by FIA with Quikchem 8500 series2 using the Quikchem method 10-107-06-2-B ammonia in surface water, wastewater, and the Quikchem method 10-116-10-1-C (tubidimetric method), respectively. Data acquisition was carried out with Omnion 3.0 from Lachat Instruments and quantification was performed based on external calibration.

A yeast strain from the prototrophic deletion collection (MATα hoΔ0::KanMX can1Δ:: STE2pr-SpHIS5 his3Δ1 lyp1Δ0) (VanderSluis et al., 2014), a Pseudomonas spp. strain (MJ020) isolated from sap (N’guyen et al., 2020) and the Leuconostoc mesenteroides strain ATCC 23386 were used in growth experiments. Each strain was individually precultured at room temperature overnight without agitation in allantoin media (Allantoin: 1.25 g/L, Yeast Nitrogen Base without amino acid and without nitrogen source (YNB): 1.75 g/L, Sucrose: 20 g/L).

To test the maximum growth capacity of microorganisms in the saps, we used filtered saps harvested in 2018 and 2019. The three strains were individually inoculated in filtered sap adjusted to 2 °Brix. Initial Optical Density (OD) was adjusted for each strain at 0.05 OD in 300 μl in a 96 well microplate and measured in triplicates, including non-innoculated negative controls for each sap. After 4 days incubation at room temperature, final OD was recorded using a Tecan Spark plate reader (Zürich, Switzerland). To assess which nutrient may be limiting, the same growth assays were performed in a subset of nine sap samples, which were spiked with either methionine (15uM), allantoate (1.25 g/L), sulfate (0.5 g/L) or sucrose (20 g/L). Corresponding non-supplemented sap were used as a control to compare each strain maximum growth.

Based on previous observations of maple sap biochemical variation using yeast functional mutants as biological reporter (N’guyen et al., 2018), we classified S_bb groups as Early (S_bb >35.5) or Late (S_bb ≤35.5). In these groups, the probability that the corresponding syrups were of acceptable quality (standard or light off-flavor) was significatively greater for the Early than the Late group (Fisher’s Exact Test, N = 32, Pvalue = 0.0004). We also characterized the syrup colors in the L*a*b* color space (N = 37). We found significant correlations between the L* (Spearman ρ = 0.68, Pvalue <0.0001) and a* (Spearman ρ = -0.55, Pvalue = 0.0004). These result shows that factors predictive of dormancy release in maple trees are also associated with syrup quality and color. Next, we characterized maple sap microbial communities to investigate how they vary along the S_bb index.

Microbial communities of 47 sap samples harvested in 2015, 2016 and 2017 were analyzed. We estimated bacterial and fungal contamination levels through 16S rDNA and ITS amplicon copy numbers determined by ddPCR, respectively. A significant correlation was observed between the S_bb index and overall microbial contamination levels (Figure 1) (Spearman ρ = -0.4243, Pvalue = 0.003). Although, bacterial and fungal contamination were correlated with each other (Spearman ρ = 0.48, Pvalue = 0.0005), only 16S rDNA copy numbers showed a significant correlation with the S_bb index (Spearman ρ = −0.43, Pvalue = 0.0026). Microbial composition and relative abundances were determined by amplicon sequencing and the total copy number measured by ddPCR was used to find estimates of absolute abundance of each community member (Supplementary Data S2). Analysis at the phylum level revealed that Proteobacteria were the most abundant, reaching 7.8 log amplicon copies/mL, but that the strongest association with the S_bb index was with Firmicutes, followed by Actinobacteria (Figure 2). For the 47 saps samples, we identified 550 bacterial and 161 fungal OTUs belonging to 77 and 64 genera, respectively (Figure 3, see Supplementary Figure S1 for an overview by sample). The most frequent and abundant bacteria were Pseudomonas, Rahnella, Janthinobacterium, Leuconostoc, Carnobacterium, Cryobacterium and Sanguibacter. Concurrently, Mrakia, Tausonia, Cadophora, Leucosporidium and Rhodosporidiobolus were the most frequently encountered fungal genera with over 70% prevalence and were also the most abundant. Molds such as Fusarium and Penicillium were also identified, but with a lower prevalence. In general, fungi ranked as high as bacteria in occurrence, but not in abundance.

Microbial community diversity varied with the S_bb index. The alpha diversity, as measured with the Shannon’s diversity was significatively correlated to the S_bb index (Spearman ρ = −0.3820, Pvalue = 0.0081), but was not different between organoleptic conformity groups (T-test Pvalue = 0.28, Figure 4A). To compare the beta diversity between S_bb and organoleptic conformity groups, we performed a PCoA ordination (Figure 4B) and multivariate statistical tests. We found that group dispersion was significatively different between S_bb groups (PERMDISP F-value: 8.7979, N = 47, Pvalue = 0.0048133), the Late group being more heterogenous. Because of this heterogeneity, we performed a non-parametric test to evaluate similarity between groups. We found that the similarity between S_bb groups is significatively lower than the similarity within the groups (ANOSIM R: 0.27239, N = 47, Pvalue <0.001). When comparing syrup organoleptic conformity, no difference in heterogeneity was detected (PERMDISP F-value: 1.0409, N = 32, Pvalue = 0.31577) and the similarity between groups was only marginally lower than similarity within groups (ANOSIM R: 0.10986, N = 32, Pvalue = 0.099). Altogether, these results indicate that microbial community profiles vary along factors influencing the tree dormancy release, but their relationship with maple syrup flavor defects may be more subtle than whole community profiles and binary groups or may require more statistical power to be supported.

Because microorganisms can alter their environment, in this case maple sap composition, the abundance of some species or combinations may be predictive of the maple syrup quality and color produced. To further investigate the relationship between OTU abundance and quantitative descriptors of sap composition, we performed a nonlinear iterative partial least squares analysis of OTUs as predictor (Figure 5). The models were generated for sap composition descriptors that are correlated with the S_bb index; ammonium, organic sulfur compounds, sulfate. For phosphate, which was not correlated to the S_bb index, the cross validation did not identify a number of factors minimizing the Root Mean PRESS. Restricting our focus on OTUs negatively correlated with the S_bb index, the strongest predictor OTUs belonged to the Proteobacteria, then Basidiomycota, Ascomycota, Firmicutes and Actinobacteria phyla. We performed the same analysis for L*a*b* values. The best predictors OTUs of maple syrup color descriptors belonged to Proteobacteria, then Basidiomycota, Ascomycota, and Firmicutes, but variable importance for projection were in general smaller than for sap composition descriptors. A few OTUs were found among the best predictors for both sap composition and syrup color parameters, namely Leuconostoc.sp_1, Pseudomonas.sp_19, Pseudomonas.sp_62 and, Pseudomonas.sp_65.

To identify the top ten predictors of maple syrup organoleptic conformity, we used a bootstrap forest partitioning approach (Figure 6A). Eight of these OTUs were more abundant in the sap of non-standard syrups (Leuconostoc.sp_1, Pseudomonas.syringae_7, Curvibasidium.cygneicollum_2, Cadophora.malorum_2, Cystofilobasidium.capitatum_1, Pseudomonas.sp_48, Carnobacterium.sp_2 and Pseudomonas.sp_62) while two where more abundant in the sap of standard syrups (Pseudomonas.sp_27, Chryseobacterium.polytrichastri_1) (Figure 6B). The fact that OTUs belonging to the Pseudomonas genus are found in both cases highlights that effect may be species- or even strain-dependent.

Among the identified top predictor OTUs of maple sap composition, syrup color and quality (Figures 5, 6A), 16 out of 45 were also significatively negatively correlated to the S_bb index (Spearman correlation, FDR adjusted Pvalue <0.05, boldface in Figures 5, 6A). Thus, only part of the microbial contaminants associated with syrup quality appears to vary along factors influencing the tree dormancy release, leaving room for other factors influencing their abundance patterns, such as microbial interactions.

To better grasp how microorganisms associate and interact in the changing maple sap environment, we constructed microbial association networks of the Early and Late groups of samples (Figure 7A). The Early network was smaller, with 270 nodes and 4211 edges compared to 485 nodes and 6332 edges for the Late network. However, the density was higher in the Early network (0.116) than in the Late network (0.055). Conversely, network heterogeneity was 0.948 in the Late network, compared to 0.651 in the Early network. Most of the OTUs present in the Early network were also present in the Late network, but few edges were shared, indicating circumstantial or context-dependent associations (Figure 7B). Although the ratios of Bacteria:Fungi were similar between networks (3:1), inter- and intra-kingdom associations proportions were significatively different (likelihood ratio Pvalue = 2.8e-131) (Figure 7C). The Early network harbored a higher proportion of inter-kingdom correlations, while the Late network had more Bacteria-Bacteria associations.

The prevalence of negative association was higher in the Late network (149, 2.4%) compared to the Early network (18, 0.4%) implicating mostly bacteria. The few negative hubs in the Early network belonged to the Pseudomonas genus, while in the Late network, members of Caulobacter, Mesorhizobium, Phyllobacterium, Bradyrhizobium, Delftia and Sphingomonas had the most negative associations. In the Early network, the abundant OTUs ranked high in betweenness centrality, but not necessarily in degree (Figure 8). In contrast, in the Late network the most abundant OTUs tended to also be the most connected, with intermediate ranks of betweenness centrality. Interestingly, among the top predictors associated with non-standard syrups, Leuconostoc.sp_1 and Pseudomonas.syringae_7 occupied hub positions in the Late network, but not in the Early network. Conversely, Pseudomonas.sp-27, which is associated with standard samples, was highly connected in the Early, but not in the Late network. These results show that microorganisms, including OTUs predictor of syrup organoleptic conformity, associate with each other’s differently in the two S_bb groups.

To complement the microbial community composition and structure analysis and identify functional features associated with the S_bb index, we performed a predictive functional analysis of maple sap bacterial communities. The results revealed 469 features strongly negatively correlated to the S_bb index (Spearman ρ < −0.50, FDR adjusted Pvalue <0.01) (Supplementary Data S3). Positive correlations with the S_bb index were all comparatively weaker (Spearman ρ < −0.36, FDR >0.01) and are thus not reported. For ease of interpretation and visualization, annotated features were grouped into functional categories (Supplementary Data S3). The top functions belonging to degradation and biosynthesis categories are shown in Figure 9. The results show that features correlated to the S_bb index are mostly unidirectional, that is they involve either biosynthesis or degradation of various metabolites, but rarely both. The strongest correlation is for sucrose degradation, the main carbon source in maple sap. The various other carbohydrates also mostly figure in the degradation categories, except for levan and peptidoglycan. Interestingly, genomic features associated with the biosynthesis of sulfur-containing secondary compounds are negatively correlated to the S_bb index. Since the abundance of organic sulfur compounds has been previously measured in these samples (N’guyen et al., 2018) and shown to be negatively correlated with the S_bb index, the predictive functional association with their biosynthesis revealed in Figure 9 suggest that some bacteria contribute to their presence in maple sap. A few Pseudomonas.sp OTUs shown in Figure 5 are correlated to both the S_bb index and organic sulfur compounds and thus are likely candidates. As for amino acids, they are split between both categories; tryptophan, leucine, isoleucine, arginine and methionine degradation, vs. asparagine, lysin, alanine and citrulline biosynthesis. For other nitrogen sources, the results include significant correlations for features related to urate and allantoin degradation as well as ammonia assimilation but also peptidases.

The increase in microbial contamination levels along the S_bb index could be explained by increasing temperatures which allow for faster growth, or by the variation in sap composition along the S_bb index. To distinguish between these two factors, we performed growth experiments with representative microorganisms (a Pseudomonas, Leuconostoc and a yeast strain) in 42 sterilized sap samples from five production sites collected in 2018 and 2019 to determine their carrying capacity at controlled temperature. The optical density measurement for the yeast (Spearman ρ = −0.55, Pvalue = 0.0002) and Pseudomonas (Spearman ρ = −0.50, Pvalue = 0.0007) strain indeed significatively increased at lower S_bb values, but not for Leuconostoc (Spearman ρ = −0.13, Pvalue = 0.40) (Figure 10A). These results demonstrate that as the tree undergoes dormancy release, the sap collected can support more growth for two of the three microorganisms tested.

In light of the current results and previous work showing that allantoate was an important nutrient for yeast in maple sap and that its quantity increased over the harvest period at the one production site tested (Filteau et al., 2017), we quantified allantoate, allantoin and ammonium in the 2018–2019 samples. We observed significant negative correlations between the S_bb index and allantoate (Spearman ρ = −0.56, Pvalue <0.0001) and ammonium (Spearman ρ = −0.54, Pvalue = 0.0002) (Figure 10B). Their concentration was relatively low and stable in the Early S_bb group and increased rapidly in the Late group. Comparatively, allantoin concentration was low and positively correlated to the S_bb index (Spearman ρ = 0.50, Pvalue = 0.0008), but this correlation was driven mostly by samples from only one production site. Therefore, the change in nitrogen content occurring after S_bb = 35.5 could explain the increased carrying capacity of maple sap.

To further understand which nutrient may influence the growth of Pseudomonas, Leuconostoc and yeast, we performed experiments with a subset of saps from three production sites supplemented with either sucrose, sulfate, methionine or allantoate (Figure 11). The growth of the Pseudomonas strain consistently increased when adding a nitrogen source (allantoate or methionine), except in the latest sample. Meanwhile, the growth of Leuconostoc was almost always increased by methionine supplementation, but not allantoate or sulfate, indicating a requirement that is more complex than sulfur or nitrogen limitation. Interestingly, the addition of sulfate to saps in the Early S_bb group, partly inhibited Pseudomonas growth. As for the yeast strain, results varied between samples, but sulfur appeared as a limiting nutrient in some cases. Altogether, these results demonstrate that Pseudomonas can degrade and use allantoate as a nitrogen source while Leuconostoc can benefit from an organic sulfur source, i.e., methionine, two metabolites that become more abundant when S_bb < 35.5. These observations thus contribute to a mechanistic understanding of the associations between the S_bb index, maple sap microbial contamination and maple syrup quality.