A defined microbial community reproduces attributes of fine flavour chocolate fermentation

Farms selection

Cocoa plantations used in this work were distributed throughout Colombia and were separated into three agroecological zones on the basis of climatic conditions, topography and soil composition. The selected agroecological zones were: (1) Santander, a mountainous region with the largest production rates of cocoa in the country, (2) Huila, an inter-Andean dry valley region and (3) Antioquia, a Pacific region. For each zone, one farm was selected for sampling: Santander (3,888 ft above sea level), Huila (3,640 ft above sea level) and Antioquia (3,993 ft above sea level) (Extended Data Fig. 1a). The Huila and Antioquia farms were located ~430 km apart, and ~550 km and 198 km away from the analysed farm in Santander, respectively (Extended Data Fig. 1a). For the farms selection, we also considered best agricultural practices, a well-established infrastructure and characteristics of cocoa fermentation. All regions had a tropical and humid climate, and the monthly average temperature during sample collection was from 22.8 °C to 28.7 °C.

Cocoa bean fermentation in the farms

All cocoa bean fermentations were performed on the farms using the farmers’ traditional practices and were analysed during the mid (May) and main harvests (October–November) on the three farms except for Antioquia, where only the main harvest period was characterized. Briefly, mature ripe cacao pods were harvested and opened manually in the fields. Cocoa beans and surrounding pulp were scooped out by hand and placed into pre-washed wooden fermentation boxes. The beans (200–400 kg) were then covered with banana leaves and/or jute bags to control the environmental temperature in the boxes. Natural cocoa bean fermentation proceeded at ambient temperatures ranging from an average minimum temperature of 22 °C (night-time) to a maximum of 33 °C (daytime). In all cases, the beans were turned at 48 h and 96 h after fermentation began. Fermented beans were removed from the fermentation box 144 h or 168 h after the start of the fermentation on the basis of the temperature, pH and cut test results.

Temperature, pH and bean colour analyses

To ensure the reproducibility of our observations, we only evaluated fermentation events that followed the traditional practices of local farmers who use wooden boxes to ferment cocoa beans. By maintaining these traditional protocols, we minimized the risk of introducing experimental bias (other types of practices such as sacks, bamboo baskets, plastic baskets, styrofoam and others) that could have altered the natural fermentation trajectory and influenced the final flavour profiles. This approach allowed us to preserve the authenticity of the fermentation process, ensuring that our results reflected the typical outcomes of the region’s traditional methods. Consequently, this consistency helped in achieving reproducibility of the fermentation process and flavour profiles under the experimental conditions. The temperature of the fermentation mass within the boxes was recorded daily using a Brannan digital thermometer (Extended Data Fig. 1b). Measurements were taken at three different zones in the box (bottom left corner, middle and top right corner) at two depths: (1) 7 cm below the top surface of the beans and (2) midway through the fermenting mass. For pH measurement, three beans were collected from the boxes, 7 cm below the top surface in each of the three selected fermentation zones. Testa (seed coats) covered with pulp were separated from the cotyledons of the beans and both tissues (testa/pulp and cotyledons) were macerated in 10 ml distilled water using a mortar and pestle. Then, the pH of the suspensions was determined using a Hanna Checker H198103 pH tester. The colour changes of the beans were determined using images taken daily with a Samsung SM-G9600 camera in automatic mode. RGB values were extracted from the bean images using ImageJ v.1.54d. A minimum of eight points on each bean image were selected for each analysed time point. Greyscale and luminance values were derived using the formulas (R + G + B)/3 and 0.299 R + 0.587 G + 0.114B, respectively. Results from the temperature and pH measurements, along with the colour changes of the beans, were collectively used to determine the end of the fermentation.

Sample collection for microbial community analyses

Samples for the microbial community analyses were collected from the fermenting beans daily using a Zymo Collection Swab (R1104) (Extended Data Fig. 1b). For the collection, we removed the top 7 cm of beans from the fermenting mass, at the centre of the fermentation box to create a small cavity, and samples were collected in duplicate by swabbing the surfaces of the beans at the bottom of the cavity. The swab buds were then placed in Zymo DNA/RNA Shield Lysis and Collection tubes (R1104), and the tube contents were mixed by vigorous shaking for 10 s. A total of n = 66 fermentation samples were collected across the three farms. Using the same protocol, swab samples were also collected from various environmental sources on the farms, including the surface of cacao leaves and pods, the inner surface of the fermentation box, the hands of farm workers involved in scooping beans and transferring them to the fermentation box, and their pod cracking tools. Soil samples were collected by discarding the top 7 cm of soil and transferring 150–250 mg of soil into Zymo DNA/RNA Shield Lysis and Collection tubes using a clean spatula. To analyse the microbiota on fruit flies around the fermentation, fruit flies were caught and incubated in 0.5 ml Zymo DNA/RNA Shield Lysis solution with agitation for 5 min. The flies were then removed, and the solution was transferred into a Zymo DNA/RNA Shield Lysis and Collection tube with an additional 0.5 ml of lysis solution. In total, n = 70 farm environment samples were collected for metagenomic analysis across the three farms. Clean swabs were used as control samples without DNA.

DNA extraction, library preparation and whole-metagenome shotgun sequencing

Samples for microbial community analysis were homogenized using the SuperFastPrep-2 cell lysis homogeniser (MP Biomedicals) at maximum speed for 10 s and were subsequently centrifuged for 1 min and 30 s at 9,800 × g. DNA was extracted using the ZymoBIOMICS DNA Miniprep kit (Zymo, D4300) following manufacturer instructions, and the resulting DNA concentration was determined using a Qubit fluorometer (Thermo Fisher). DNA sequencing libraries were prepared using the Rapid PCR Barcoding kit (SQK-RPB004) from Oxford Nanopore Technologies (ONT). In brief, 1 µl fragmentation mix was added to 3 µl DNA (2–10 ng µl⁻¹), and the reaction was mixed by gentle finger-flicking. The tube was placed in a miniPCR mini16 thermal cycler (Amplyus) and the DNA was fragmented using the following conditions: 30 °C for 1 min, then 80 °C for 1 min. The fragmented DNA was cooled and amplified in a PCR reaction containing 20 µl nuclease-free water, 25 µl LongAmp Taq 2× master mix (New England Biolabs (NEB), M0287L), 4 µl fragmented DNA and 1 µl barcode adaptor. The reaction was gently mixed and amplified using the following conditions: 95 °C for 3 min, 20 cycles of denaturation at 95 °C for 15 s, annealing at 56 °C for 15 s and extension at 65 °C for 6 min, and a final extension of 65 °C for 6 min. The resulting DNA library was purified using 0.6× Agencourt AMPure XP beads (Beckman Coulter, A63881) and eluted in 10 µl 10 mM Tris-HCl pH 8.0 and 50 mM NaCl. The library concentration was determined using a Qubit fluorometer (Thermo Fisher). Equimolar quantities of individual barcoded sample libraries were pooled and the volume adjusted to 10 µl using 10 mM Tris-HCl pH 8.0 and 50 mM NaCl. Subsequently, 1 µl of Rapid Adapter solution was added to the pooled library and the tube was incubated at room temperature for 5 min. Then, 34 µl sequencing buffer, 25.5 µl loading beads and 4.5 µl nuclease-free water were added to the tube, and the contents were mixed gently. The prepared pooled library was added to a verified and primed FLO-MIN106 R9.4.1 flow cell (ONT, FLO-MIN106D) in a MinION DNA sequencer (ONT) following manufacturer instructions. DNA sequencing was conducted with default parameters using MinIT (ONT) with MinKNOW v.2.1.12 (ONT). Fast5 files were base called with Guppy v.4.0.15 using the ‘template_r9.4.1_450bps_hac.jsn’ high-accuracy model (ONT).

Temperature, pH and bean colour analysis

Temperature and pH distributions were compared using the Kolmogorov–Smirnov test with the ks.test() function from the stats v.4.3.0 package in R and plotted with ggplot2 v.3.4.2. To explore the correlations between temperature and pH in the testa/pulp and cotyledons, scatterplots with correlation coefficients and P values were generated using the stat_cor() function from the ggpubr v.0.6.0 package. To assess dissimilarities in bean colour at the different time points during fermentation, principal component analysis (PCA) was performed with the prcomp() function in R, and the first two principal components were plotted using ggplot2 v.3.4.2. In addition, scatterplots with correlation coefficients and P values were employed to explore correlations between bean colour values, temperature and pH.

Processing and filtering of whole-metagenome shotgun sequence data

We obtained a total of 29,714,777 base-called reads (97.7 Gbp) from the whole-metagenome shotgun sequencing across the three farms. The initial dataset underwent demultiplexing, and primer and barcode sequences were trimmed using qcat v.1.1.0 (ONT). Reads with ambiguous barcode assignments were excluded from further analysis. The reads were filtered with NanoFilt (v.2.8.0)⁵⁹ to discard low quality sequences (Q-score < 9) and sequences <100 bp. Reads were mapped to the Theobroma cacao Criollo v.2.0 reference genome⁶⁰ as well as the Homo sapiens reference genome GRCh38.p14 (RefSeq GCF_000001405.40) using minimap2 (v.2.17)⁶¹ to identify and remove potential contaminating DNA in samples. Reads that mapped were removed using SAMtools (v.1.9)⁶² and Seqtk v.1.3 scripts. Following these processing steps, we retained 24,300,820 reads (80.8 Gbp) with an average read length of 3,326.4 bp (Extended Data Fig. 1h) and a mean read quality of Q13.2.

Profiling microbiota community composition

We used the Kraken v.2.1.2 pipeline⁶³ for classifying the whole-metagenome shotgun sequencing reads. The reads were classified using the Kraken 2 archaea, bacteria, viral, plasmid, human, UniVec_Core, protozoa and fungi reference databases (k2_pluspf_20220607). To estimate relative abundances, the Bracken v.2.7 pipeline⁶⁴ was applied to the classification results. Subsequently, Pavian v.1.0 facilitated the extraction of abundance and taxonomic tables. Functions in phyloseq v.1.44.0 with microbiome v.1.22.0 and microbiomeutilities v.1.0.17 were used to filter the dataset and remove samples with low read depth, remove unidentified taxa and singletons, transform abundance values using rarefaction, subset and merge sample and taxonomic groups, and perform other dataframe manipulations. To assess alpha diversity across the samples, we calculated the Shannon diversity index using phyloseq v.1.44.0. We used analysis of variance (ANOVA) to test for significant differences in Shannon diversity indices between groups, and means were separated using Tukey’s honestly significant difference (HSD) test from the agricolae v.1.3.5 R package. For beta diversity, Bray–Curtis dissimilarity matrices were calculated using the phyloseq v.1.44.0 ‘bray’ method, and the variances explained by fermentation time, farm location and harvest season were estimated by performing permutational multivariate analysis of variance (PERMANOVA) using the adonis2() function in the vegan v.2.6.4 R package. Unconstrained and constrained ordination of beta diversity was plotted using principal coordinate analysis (PCoA) and canonical analysis of principal coordinates (CAP), respectively, on the basis of Bray–Curtis dissimilarity matrices calculated with vegan v.2.6.4. We visualized differences in fermentation time, farm location and harvest season with the CAP analysis, using the following models:

The relative abundance of taxa was plotted as a stacked bar representation using phyloseq v.1.44.0. The tax_glom() function in phyloseq v.1.44.0 was used to agglomerate taxa, and the aggregate_rare() function in microbiome v.1.22.0 was used to aggregate rare groups. Mantel correlations between bacterial and fungal Bray–Curtis dissimilarity matrices were performed using the mantel() function of vegan v.2.6.4, with the Pearson method and 10,000 permutations. We used DESeq2 (v.1.40.0)⁶⁵ to calculate the enrichment profiles at Santander by fitting a generalized linear model (GLM) with the following design:

We extracted the following comparisons from the fitted model: 24 h vs 0 h, 48 h vs 0 h, 72 h vs 0 h and 96 h vs 0 h. Taxa at the family, genus and species levels were considered significant if they had a false discovery rate (FDR)-adjusted P value (q value) < 0.05. The results of the GLM analysis were rendered in heat maps coloured on the basis of the log₂ fold change output by the GLM. Significant differences between comparisons with a q value < 0.05 with log₂ fold change > ±2 were highlighted with black squares.

Single-nucleotide polymorphism (SNP) genotyping of cacao varieties

To elucidate the genetic backgrounds of the cacao varieties cultivated across the three cocoa farms, we first conducted a survey focusing on diverse fruit morphologies to gauge the diversity present. The parameters assessed encompassed fruit characteristics such as form, basal constriction, apex, rugosity, ridging, length, diameter, wall thickness and the colour of the mature unripe fruit. Following this morphological survey, young and healthy leaf samples were collected from each distinct morphotype identified (n_total = 24; Santander n = 12, Huila n = 5 and Antioquia n = 7). The leaf samples were washed and dried, and DNA isolation was carried out using a ZR Plant/Seed DNA MiniPrep kit (Zymo, D6020) with modifications detailed in ref. ⁶⁶. The samples were genotyped at 96 SNP sites (Supplementary Table 1) on a Fluidigm Juno System using a Juno 96.96 Genotyping IFC (Standard BioTools) in accordance with manufacturer instructions. The SNP sites analysed were selected from the Theobroma cacao global reference SNP panel^67,68. Briefly, genotyping assays were prepared using the Juno SNP Type Genotyping Reagent kit (Standard BioTools, 100-8364) and specific SNP type genotyping assays manufactured by Standard BioTools. The Juno 96.96 thermal cycling protocol included a multiplex specific target amplification (STA) step before the SNP genotyping to enrich the template molecules. STA thermal cycling conditions comprised 95 °C for 2 min, followed by 14 cycles at 95 °C for 15 s, and 60 °C for 4 min. For SNP genotyping, reactions were initiated at 95 °C for 10 min, followed by 4 cycles at 95 °C for 15 s, 64–61 °C (1 °C decrease with each cycle) for 45 s, and 72 °C for 15 s. This was followed by 39 cycles at 95 °C for 15 s, 60 °C for 45 s, and 72 °C for 15 s. Fluorescence intensity was quantified using the Fluidigm EP1 software (Standard BioTools), and genotypic calls were automatically made using Fluidigm SNP Genotyping Analysis software v.4.1.3 (Standard BioTools). SNP data generated are available in Supplementary Table 2.

Drying, roasting and sensory evaluation of liquor samples

The fermented beans were spread on a wooden surface in a 3–4-cm layer for sun drying. The drying mass was mixed every 1.5 h for the first 3 days of drying and every 3 h from the fourth day onward. Drying was carried out carefully to ensure that off-flavours did not develop. The beans were covered overnight and during rainy conditions. All batches underwent drying until reaching a final moisture content of 7%, taking ~7–8 days, after which they were stored in jute bags. Quality assessments were conducted on 100-g samples from each bean lot using cut tests following the procedure outlined in ref. ⁴¹. The evaluation included identifying characteristics such as underfermented (purple/violet), overfermented (grey/slaty), insect-damaged, chopped/broken, germinated, mouldy, double, or flat bean defects. Beans categorized as fully fermented with no defects and those that were partially purple were considered desirable/high quality. For the roasting process, the beans were placed on shallow perforated trays and roasted at 140 °C for 25 min a convection oven (Sheldon Manufacturing). Subsequently, the beans were cooled to ambient temperature, broken and winnowed to produce cocoa nibs. The nibs were transformed into cocoa liquor on a granite-wheeled melangeur (CocoaTown), reducing the particle size to 16–18 μm. Liquor samples (60 °C) were evaluated through coded, randomized tastings by 3–6 trained sensory panelists in duplicate or triplicate. The panelists consisted of members of the Food Technology, Quality and Sensory Evaluation team at the Cocoa Research Centre, Trinidad. The panel members were trained in accordance with the Cocoa of Excellence Programme guidelines⁶⁹ under the supervision of the panel leader, and are experienced in cocoa sensory evaluation. Flavour descriptors assessed by the panel members were based on the cocoa liquor sensory evaluation template of E.S. Seguine and D.A. Sukha⁷⁰ and expressed as numerical values between 0 and 10. Reference liquors from Madagascar (fine or flavour cocoa) and Ivory Coast and Ghana (bulk cocoa) were included in the sensory analysis. To neutralize palates between tastings, soda crackers and mouth rinsing with still water were employed.

Analysis of cacao genotypes

SNP genotypic data were generated at 96 SNP sites for each morphotype in the study (n_total = 24; Santander n = 12, Huila n = 5 and Antioquia n = 7). This dataset was then combined with the SNP profiles from 228 cacao reference accessions sourced from the Cocoa Research Centre, Trinidad, SNP database. The reference SNP profiles were primarily generated from cacao accessions at the International Cocoa Genebank Trinidad and were selected across the 10 cacao genetic clusters identified in ref. ⁷¹ (n_total = 228; Amelonado n = 28, Contamana n = 15, Criollo n = 15, Curaray n = 22, Guiana n = 24, Iquitos n = 22, Marañon n = 27, Nacional n = 17, Nanay n = 23 and Purús n = 5). In addition, 30 Amelonado-Criollo hybrid accessions were included. The combined dataset was filtered by removing SNPs with >10% missing data and monomorphic SNPs. The result was a final dataset of 84 high-quality SNP markers, with a missing data range between 0 and 3.57% and a mean of 0.52% across all accessions. For phylogenetic analysis, SNP profiles were converted into DNA strings, aligned using DECIPHER v.2.24.0 and transformed into a distance matrix with seqinr v.4.2.16. A neighbour-joining tree was constructed with ape v.5.6.2, and the resulting tree was visualized using ggtree v.3.8.0 with ggtreeExtra v.1.10.0. The genetic distances were further analysed through a PCoA for visualization. The PCoA involved converting genetic distances to 2 dimensions using classical multidimensional scaling with the stats v.4.3.0 package, and plotting with ggplot2 v.3.4.2. Ancestry was inferred using STRUCTURE (v.2.3.4)⁷², employing structure-threader (v.1.3.10)⁷³ for parallelized runs across multiple CPU cores. To facilitate the analysis, reference accessions were replicated to ensure a minimum representation of 60 individuals for each of the 10 cacao genetic groups. Simulations were calculated using the admixture model with alpha inferred and independent allele frequency with 200,000 burn-ins and 500,000 Monte Carlo Markov Chain repetitions without any previous genetic or geographic origin information. The number of clusters (K) was set from 8 to 12 with 30 iterations for each K value. CLUMPAK (Cluster Markov Packager Across K) (v.1.1)⁷⁴ was utilized to assess the congruence among independent STRUCTURE runs for each K value, and the optimum K value was determined according to ref. ⁷⁵.

Community-wide microbial source tracking analysis

To explore how the surrounding microbial environmental sources in the cocoa plantations may be contributing to cocoa fermentation communities, we used FEAST v.0.1.0 to perform community-wide microbial source tracking analysis. The results were plotted with ggplot2 v.3.4.2.

Statistical analysis of bean quality and cocoa liquor sensory profiles

Bean quality assessments were analysed using Fisher’s exact test with the fisher.test() function of the stats v.4.3.0 package. For the analysis of sensory attributes in the cocoa liquors, we estimated the variances explained by farm location, harvesting period and sensory panelists by performing PERMANOVA using distance matrices with the adonis2() function in the vegan v.2.6.4 R package. This analysis allowed us to examine how location, harvest and panelist, and their interactions, contributed to the variation in the sensory data (Fig. 3a). A constrained ordination of the sensory attributes was plotted using CAP with vegan v.2.6.4 using the following model:

To illustrate the sensory characteristics of individual cocoa liquors, the mean scores for each sensory attribute were calculated across panelists for each liquor sample. These scores were then transformed to a scale between 0 and 6 using the rescale() function of the scales v.1.2.1 R package. The transformed scores were visualized on a heat map generated with ggplot2 v.3.4.2. Hierarchical clustering of sensory attributes was applied using the ward.D2 or single method within the hclust() function in R. The clustering was based on Euclidean distances calculated using the dist() function on the transformed scores.

Extracting abiotic kinetic features and random forest analysis

To identify abiotic features associated with the sensory attributes of the cocoa liquors, we utilized the Practical Program for Forces Modeling (PPFM 2020)⁷⁶ tool to model the kinetics of the temperature changes during bean fermentation across the three locations. This involved randomly selecting a minimum of 15 temperature versus time data points and fitting the temperature curve using a 5-parameter general model equation. From the model, we derived several key features including: (1) maximum growth rate (peak rate at which the system grows during a specified period); (2) time to maximum kinetic energy (duration for the system to reach its maximum kinetic energy level); (3) temperature at maximum kinetic energy (the specific temperature value at the point of maximum kinetic energy); (4) exponential phase duration (period of rapid increase in numbers or activity); (5) linear phase duration (phase where the kinetics rate becomes relatively constant); (6) exponential decay phase duration (timeframe when the system starts to decline after reaching its maximum kinetics); (7) temperature change during exponential phase (change in temperature during rapid exponential kinetics); (8) temperature change during linear phase (alterations in temperature during stable kinetics); (9) temperature change during exponential decay (variations in temperature during the decline following exponential kinetics); (10) rate of temperature change during the exponential phase (speed at which temperature changes during rapid exponential kinetics); (11) rate of temperature change during the linear phase (speed at which temperature changes during stable kinetics); (12) time to inflection point (duration for the system to reach the inflection point, indicating a shift in kinetics pattern); and (13) inflection point (point on the kinetics curve where the curvature changes, signifying a transition in kinetics rate or pattern). The entire process was repeated at least three times for each farm’s fermentation. In addition, we extended our analysis to model the kinetics of the inverse cotyledon pH, extracting similar curve features. The feature values were normalized using the rescale() function in the scales v.1.2.1 R package. The mean normalized feature values were then visually represented on a heat map using ggplot2 v.3.4.2. Subsequently, Pearson correlation coefficients and corresponding P values between these features were computed using the rcorr() function in the Hmisc v.5.0.1 package. The results of the correlation analysis were graphically presented using ggplot2 v.3.4.2, where the colour of the plots reflected the correlation coefficient values. Significant correlations (P < 0.05) were emphasized with black squares on the plots. Furthermore, the coefficient of variation for the feature values was calculated and depicted using ggplot2 v.3.4.2. The 3 plots were integrated on the basis of the hierarchical clustering of the Pearson correlation coefficients of the features. The clustering employed the ward.D2 method within the hclust() function in R, utilizing Euclidean distances calculated using the dist() function. For each cluster identified, we selected the feature with the highest coefficient of variation as a representative for the cluster. Following this, for each sensory attribute, we employed the randomForest v.4.7.1.1 R package to construct a random forest model. This was done to pinpoint the most significant features associated with each sensory attribute. Subsequently, we visualized the percentage increase in mean squared error (%IncMSE) for each feature by generating a heat map using ggplot2 v.3.4.2. Hierarchical clustering of the feature importance was applied using the ward.D2 method within the hclust() function in R. The clustering was based on Euclidean distances calculated using the dist() function.

Extracting taxonomic kinetic features and random forest analysis

To identify taxonomic markers associated with the sensory attributes of the cocoa liquors, we began by pinpointing the pivotal bacteria and fungi responsible for the notable variances in beta diversity, specifically focusing on those exerting the most influence on beta diversity disparities observed among fermentation time points and across the three farm locations. To identify the top bacteria and fungi driving the differences in beta diversity across fermentation time and farm location, we calculated PERMANOVA coefficients of the taxa using the adonis() function in the vegan v.2.6.4 R package at the genus level, and assessed their prevalence with microbiome v.1.22.0. The results were visualized using ggplot2 v.3.4.2. Bacteria with coefficients ≥10 and fungi with coefficients ≥5 were selected for further analysis. We verified that the selected bacteria and fungi could recapitulate the differences in the beta diversity of the entire community by performing a PERMANOVA with the selected taxa, as well as a CAP, on the basis of Bray–Curtis dissimilarity matrices calculated with vegan v.2.6.4 using the following models:

Subsequently, the relative abundances of the selected bacteria and fungi were extracted from the metagenomic dataset for each fermentation across the three farms and different fermentation time points. Following this, we utilized the gcplyr v.1.5.2 R package to extract growth features of the selected taxa. The extracted features encompassed: (1) first local maxima (the initial peak density achieved during growth before a subsequent decline); (2) initial density (starting density of taxa, corresponding to the first local minima in taxa density); (3) area under the curve (overall taxa growth); (4) maximum density (peak taxa density, offering insights into the taxa carrying capacity within a specific environment, or alternatively, measures of taxa growth yield or efficiency); (5) time to maximum density (duration taken to reach the maximum taxa density in the environment); (6) midpoint (the moment when the density initially reaches half of the maximum density); and (7) inflection point (the instance when the derivative of the growth curve attains its maximum value). The feature values were normalized using the rescale() function in the scales v.1.2.1 R package. The mean normalized feature values were then visually represented on a heat map using ggplot2 v.3.4.2. Hierarchical clustering of the taxa was performed using the hclust() function with the ward.D2 method, on the basis of Euclidean distances calculated with the dist() function. Subsequently, Pearson correlation coefficients and corresponding P values between these features were computed using the rcorr() function in the Hmisc v.5.0.1 package. The results of the correlation analysis were graphically presented using ggplot2 v.3.4.2, where the colour of the plots reflected the correlation coefficient values. Significant correlations (P < 0.05) were emphasized with black squares on the plots. Furthermore, the coefficients of variation for the feature values were calculated and depicted using ggplot2 v.3.4.2. The 3 plots were integrated on the basis of the hierarchical clustering of the Pearson correlation coefficients of the features. The clustering employed the ward.D2 method within the hclust() function in R, utilizing Euclidean distances calculated using the dist() function. For each cluster identified, we selected the feature with the highest coefficient of variation as a representative for the cluster. Following this, for each sensory attribute, we employed the randomForest v.4.7.1.1 R package to construct a random forest model. This was done to pinpoint the most significant features associated with each sensory attribute. Subsequently, we visualized the percentage increase in mean squared error (%IncMSE) for each feature by generating a heat map using ggplot2 v.3.4.2. Hierarchical clustering of the feature importance was applied using the ward.D2 method based on Euclidean distances.

Validation of abiotic and taxonomic markers linked to sensory attributes in cocoa liquors

To assess the robustness of the association between abiotic and taxonomic markers and sensory attributes in cocoa liquors, we modelled the kinetics of temperature and pH changes during bean fermentation across 19 independent fermentations, conducted on cocoa farms in diverse agroecological regions of Trinidad between 2018 and 2022. The raw data, including bean temperature, pH and corresponding bean flavour profiles from fermentations, was obtained from the internal database of the Cocoa Research Centre, The University of the West Indies. Farms and fermentation events were selected to capture the full spectrum of cocoa flavour profiles found in Trinidad, a country renowned for producing high-quality fine or flavour cocoa beans. From the kinetic curves, we extracted the following features: temperature inflection point, time to temperature inflection point, duration of the temperature exponential phase, duration of the temperature exponential decay phase, rate of temperature change during the exponential phase, pH exponential decay phase duration, and pH change rate during the exponential phase. In addition, the relative abundances of the selected bacterial and fungal taxa were extracted from metagenomic data of 11 fermentations. Subsequently, growth curves were plotted as described earlier, and the following growth features were extracted: area under the curve, inflection point, initial density and midpoint. As before, we constructed random forest models to identify the most important abiotic and taxonomic features associated with each sensory attribute. Feature importance was visualized using heat maps, displaying the percentage increase in mean squared error (%IncMSE).

Construction of metagenome assembled genomes

Various strategies were utilized to construct the MAGs. Initially, a single-sample assembly and binning approach was adopted, where reads from individual samples were assembled into contigs using metaFlye⁷⁷ in the Flye v.2.9 package with default mode. Reads from each sample were subsequently mapped to the respective assembly using minimap2 (v.2.17)⁶¹, and the corresponding abundance files were generated using SAMtools (v.1.12)⁶². The abundance files were used for metagenomic binning of the contigs using two different binning tools: MaxBin (v.2.2.4)⁷⁸ with default parameters and MetaBAT (v.2.15)⁷⁹ with specific parameters (percentIdentity=85, minContigLength=1000, minContigDepth=1). To help capture low-abundance microbes, a co-assembly and binning approach was implemented. This involved pooling reads from fermentation samples within each region (Santander, Huila, Antioquia), assembling contigs and generating metagenomic bins. In addition, a concatenation of reads from fermentation samples across all three farms was performed, followed by contig assembly and binning. These diverse strategies resulted in the construction of 1,591 MAGs. Subsequently, MAGs underwent dereplication using dRep (v.3.4.0)⁸⁰, with genome filtering options set at 10,000 bp minimum length, 10% minimum completeness, 10% maximum contamination and 95% average nucleotide identity (ANI) threshold for species-level dereplication (see ref. ⁸¹ for species-level definition). The dereplicated MAGs were evaluated using CheckM (v.1.1.6)⁸² to determine their genome completeness and contamination levels. MAGs were assigned to be either low, medium or high quality based on the completeness and contamination levels recommended by ref. ⁸³ (low-quality: completeness 0–50%, contamination <10%; medium-quality: completeness 50–90%, contamination <10%; high-quality: completeness >90%, contamination <5%). After excluding MAGs with contamination levels >10% and those with completeness <50% for bacteria or 30% for fungi, 55 MAGs were retained for further analysis. The completeness and contamination statistics for the final MAG set were visualized using ggplot2 v.3.4.2 (Supplementary Fig. 4a,b), and additional quality and genome statistics can be found in Supplementary Table 4. To assess how well the MAGs represented the fermentation and farm environment communities, sequence reads were mapped to the MAGs using minimap2 (v.2.17)⁶¹, and mapped reads were identified with SAMtools (v.1.9)⁶² and Seqtk v.1.3 (Supplementary Fig. 4c).

MAG classification, annotation and abundance

The taxonomic classification of the MAGs was performed using the CAT v.8.22 taxonomic classification pipeline⁸⁴. This process entailed identifying open reading frames (ORFs) within each contig, followed by mapping the predicted ORFs against the NCBI NR protein database. The taxonomic assignment of the MAGs was determined on the basis of the consensus classification of individual ORFs. Taxonomic classification of the MAGs can be found in Supplementary Table 4. To visualize the relationships among the MAGs, a dendrogram was constructed using the neighbour-joining approach, utilizing marker gene sequences identified across the genomes of the 55 MAGs from the CheckM (v.1.1.6)⁸² tool. The marker gene sequences for each MAG were initially merged into a string and aligned using Clustal W in the msa v.1.32.0 package. The alignment was trimmed with microseq v.2.1.6, transformed into a distance matrix with seqinr v.4.2.16, and a neighbour-joining tree was constructed with ape v.5.6.2. The resulting tree was visualized using ggtree v.3.8.0 with ggtreeExtra v.1.10.0 (Fig. 4a). The ribosomal RNA (rRNA) genes in the MAGs were identified using Barrnap v.0.9 (https://github.com/tseemann/barrnap). Open reading frames from each MAG were predicted using FragGeneScanRs (v.1.1.0)⁸⁵ with default settings. Functional annotation of predicted proteins was performed using eggNOG-mapper (v.2.1.9)⁸⁶ with the eggNOG v.5.0.2 database⁸⁷ with Diamond (v.2.0.11)⁸⁸ and MMseqs2 release 12-113e3 (ref. ⁸⁹). To assess the relative abundance of the MAGs during fermentations, reads from each fermentation sample across the three regions were mapped to the MAGs using minimap2 (v.2.17)⁶¹, and CoverM (v.0.6.1)⁹⁰ was used to extract the relative abundance counts with the ‘genome’ mode and transcripts per million coverage method. Subsequently, DESeq2 (v.1.40.0)⁶⁵ was used to determine the enrichment profiles of the MAGs in the three farms by fitting a GLM with the design:

We extracted the following comparisons from the fitted model: 48 h vs 0 h, 72 h vs 0 h and 96 h vs 0 h. Significance was determined with an FDR-adjusted P value (q value) < 0.05. The results of the GLM analysis were rendered in a heat map coloured on the basis of the log₂ fold change. Significant differences between comparisons (q value < 0.05) with log₂ fold change > ±2 were highlighted with black squares (Supplementary Fig. 4d).

Enrichment of microbial biological functions during fermentation

To identify metabolic processes that were enriched within the microbial communities during the fermentation, we selected the contigs assembled from individual samples (380,365 contigs) and subsequently subjected them to deduplication using the dedupe.sh tool in BBTools v.38.76 to eliminate redundancies. Next, we determined the relative abundance of the contigs during fermentations by mapping the reads from the samples to the contigs using minimap2 (v.2.17)⁶¹, and extracting the relative abundance counts using CoverM (v.0.6.1)⁹⁰ in the ‘contig’ mode and reads_per_base coverage method. Taxonomic classification of the contigs was performed using the CAT v.8.22 taxonomic classification pipeline⁸⁴. Subsequently, the contigs were filtered to retain only bacterial and fungal sequences. DESeq2 (v.1.40.0)⁶⁵ was utilized to determine the contig enrichment profiles in the three farms by fitting a GLM with the following design:

We extracted the following comparisons from the fitted model: 48 h vs 0 h, 72 h vs 0 h and 96 h vs 0 h. Contigs meeting the criteria of an FDR-adjusted P value (q value) < 0.05 and a log₂(fold change) > ±2 were selected for further analysis. Open reading frames encoded within the contigs were predicted using FragGeneScanRs (v.1.1.0)⁸⁵ with default settings. This was followed by functional annotation of the predicted proteins using the eggNOG-mapper (v.2.1.9)⁸⁶ pipeline with the eggNOG v.5.0.2 database⁸⁷ with Diamond (v.2.0.11)⁸⁸ and MMseqs2 release 12-113e3 (ref. ⁸⁹). Enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) biochemical reactions, along with associated enzymes and proteins, were identified on the basis of an adjusted P-value threshold of <0.01 and visualized using a heat map generated with ggplot2 v.3.4.2. The genes annotated with Gene Ontology (GO) classifications were subsequently extracted, and a GO enrichment analysis focusing on biological processes was conducted. This involved employing adaptive GO clustering in conjunction with Mann–Whitney U testing, utilizing the GO_MWU tool as previously described⁹¹. In this analytical approach, genes were ranked on the basis of signed log₂ fold change values. Significantly enriched and depleted GO categories were determined by an adjusted P value < 0.05 (Supplementary Table 5). The most prominent enriched and depleted GO categories shared across comparisons were visualized in ggplot2 v.3.4.2 and coloured on the basis of the square root transformed delta rank values (enrichment score) of the GO categories (Supplementary Fig. 5). Similarly, enriched biochemical reactions annotated in the KEGG database, and proteins, were identified using a generalized linear model with an adjusted P-value threshold of <0.05 (Supplementary Tables 6 and 7). Prominently enriched and depleted categories shared across comparisons were visualized as heat maps (Supplementary Fig. 6).

Metabolic network analysis and identification of a defined microbial community

For each MAG, we combined the predicted coding DNA sequence and corresponding translated amino acid sequences identified with FragGeneScanRs (v.1.1.0)⁸⁵ with the functional annotations predicted by eggNOG-mapper (v.2.1.9)⁸⁶, into genbank-formatted files using emapper2gbk (v.0.3.0)⁹² in ‘genes mode’. Subsequently, these files were utilized to generate the metabolic network of the fermentation community using the Metage2Metabo (v.1.5.3)^92,93 graph-based metabolic analysis pipeline. Briefly, genome-scale metabolic networks (GSMNs) were reconstructed for the 44 MAGs detected in the fermenting bean using the Metage2Metabo (v.1.5.3)^92,93 pipeline with Pathway Tools (v.26.0)⁹⁴. The networks were then analysed to determine individual metabolic capabilities and, subsequently, the collective metabolic capabilities of the community. Metabolites known to be present in cocoa pulp (artificial pulp components: 0.14% (w/v) high-viscosity carboxymethyl cellulose, 0.77% (w/v) low-viscosity carboxymethyl cellulose, 1.09% (w/v) pectin, 2.5% (w/v) sucrose, 4% (w/v) glucose, 5% (w/v) fructose, 1% (w/v) citric acid, 0.5% (w/v) yeast extract, 0.5% (w/v) peptone, 0.1% (w/v) calcium lactate pentahydrate, 0.1% (v/v) Tween 80, 0.05% (w/v) magnesium sulfate heptahydrate and 0.02% (w/v) manganese sulfate monohydrate) (see refs. ^27,50) were used as metabolic precursors to seed the network. The cooperation potential between GSMNs was assessed by calculating the added value of metabolic cooperation within the community. The added value of cooperation was used as the metabolic target to compute the key species and a defined community. Defined microbial communities were then identified by simplifying the complexity of the full community into a defined community with equivalent metabolic capabilities. The metabolites reachable by each MAG, identified on the basis of the cocoa pulp metabolic precursors seeded in the network, were compiled into a data matrix. This matrix was visualized using ComplexHeatmap v.2.12.1 with circlize v.0.4.15 (Supplementary Fig. 7). To visualize changes in the abundance of microbes with different metabolic potential over time, the data were transformed into a distance matrix using the dist() function in R⁹⁵ with the ‘euclidean’ method and converted to two dimensions using classical multidimensional scaling of the dissimilarities with the stats v.4.3.0 package. Results were plotted with ggplot2 v.3.4.2 and coloured on the basis of the relative abundance of each MAG during the fermentation at each time point for each region (Extended Data Fig. 3a). Finally, our metabolic network analysis identified 10 MAGs possessing metabolic capabilities equivalent to the full community. The metabolites reachable by the 10 MAGs, based on the cocoa pulp as the precursor, were compiled and transformed into a distance matrix, then converted to two dimensions and plotted with ggplot2 v.3.4.2 as described above. We used ggvenn v.0.1.10 to visually represent the count of metabolites generated by the microbial communities, as well as to demonstrate the overlap of metabolites shared among them.

Isolation of bacterial and fungal strains from fermenting cocoa beans for in vitro studies

We isolated bacterial and fungal strains from fermenting cocoa beans from the Cocoa Research Centre’s fermentation facility at the University of the West Indies in Trinidad. For this, mature ripe cacao pods were harvested from the International Cocoa Genebank, Trinidad (ICGT), and opened manually in the field by trained staff. The beans and surrounding pulp were scooped out by hand, placed into clean plastic bags and transported to the Cocoa Research Centre fermentation facility. The beans were placed into a pre-washed wooden fermentation box, covered with jute bags to minimize heat loss, and allowed to undergo natural fermentation at ambient temperatures ranging from a minimum of 22 °C (night-time) to a maximum of 34 °C (daytime). The beans were turned at 48 h and 96 h after fermentation initiation. Using sterile surgical gloves, two beans were collected daily from the fermentation mass at a depth of 7 cm at the centre of the fermentation box. These beans were cut into small pieces with a sterile scalpel blade. A few pieces of each bean were placed into tubes containing 80% glycerol, and the tubes were stored at −80 °C until further used. To culture the isolates, the bean–glycerol mixture was homogenized using a sterile scalpel. Samples (1 ml) of the homogenate were serially diluted (1/100, 1/1,000, 1/10,000, 1/100,000) with 10 mM MgCl₂, from which 100-µl aliquots of the dilutions were plated on different selective agar media. The selective media included: (1) acetic acid medium (AAM): 1% (w/v) d-glucose (Sigma-Aldrich, G7021-5KG), 0.5% (v/v) ethanol, 0.3% (v/v) acetic acid, 1.5% (w/v) bacteriological peptone (Millipore) (Sigma-Aldrich, 91249-500 G), 0.8% (w/v) yeast extract (Millipore) (Sigma-Aldrich, 70161-500 G), 2% (w/v) agar (Sigma-Aldrich, A6686-500G), pH 4.5 with nystatin (30 mg l⁻¹) (Sigma-Aldrich, N3503-25MU) and penicillin (50 mg l⁻¹) (Sigma-Aldrich, 13752-5G-F); (2) De Man Rogosa Sharpe agar (MRS) (Millipore)⁹⁶ (Sigma-Aldrich, 69964-500 G) with 0.1% (v/v) Tween 80 (Sigma-Aldrich, P4780-100ML) and nystatin (30 mg l⁻¹); (3) yeast peptone glucose agar (YPG): 1% (w/v) yeast extract, 2% (w/v) peptone, 2% (w/v) glucose, 2% (w/v) agar, pH 5.6 with chloramphenicol (100 mg l⁻¹) (Sigma-Aldrich, C1919-25G); and (4) nutrient agar (NA) (Millipore) (Sigma-Aldrich, 70148-500G). These media were optimized for the culture of the main taxonomic groups present in fermenting cocoa beans. The plates were incubated at various temperatures (25 °C, 30 °C, 37 °C and 42 °C) for 1–3 weeks. Colonies with distinct morphologies based on colony appearance, colour, optimal growth temperature and growth rate were selected and purified through successive subculturing until no visible signs of contamination were observed. The purified isolates were grown in liquid culture, mixed at a ratio of 1:1 with 80% (v/v) glycerol and stored at −80 °C for future use.

Identification of isolates using Sanger sequencing

To identify the isolates, we conducted amplification and sequencing of the 16S rRNA gene for bacteria and the internal transcribed spacer (ITS) region for fungi that were cultured. A single colony of each isolate was inoculated into LB media (500 µl) and incubated overnight at 28 °C with agitation at 200 r.p.m. Subsequently, 10 µl of the culture was heated at 95 °C for 5 min and then centrifuged at 10,000 × g for 1 min. The supernatant (1 µl) was used in the following PCR reaction mix: 5.4 µl Milli-Q water, 2 µl 5× Phusion HF buffer, 0.6 µl 25 mM MgCl₂, 0.2 µl 10 mM dNTP, 0.2 µl DMSO, 0.25 µl of each primer (20 pmol µl⁻¹), 0.05 µl Phusion High-Fidelity DNA polymerase (NEB, M0530L) and 1 µl DNA. The reactions were initially heated to 94 °C for 3 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at the optimized primer temperature (16S rRNA: 58 °C, ITS: 64 °C) for 30 s, extension at 72 °C for 1 min and 30 s, and a final extension at 72 °C for 10 min. The V1–V9 region of the bacterial 16S rRNA gene was amplified with the 8F forward primer (5′-AGAGTTTGATCCTGGCTCAG-3′) and fD1 reverse primer (5′-ACGGCTACCTTGTTACGACTT-3′), while the ITS region was amplified using the ITS1 forward primer (5′-TCCGTAGGTGAACCTGCGG-3′) and ITS4 reverse primer (5′-TCCTCCGCTTATTGATATGC-3′). To confirm the specific amplification of target DNA regions, half of the PCR volume was visualized on an agarose gel via electrophoresis. The resulting amplicons were prepared for sequencing by combining 14 µl of Milli-Q water, 1 µl of the PCR reaction and 2 µl of sequencing primer (10 pmol µl⁻¹). Sanger sequencing of the PCR products was conducted using the 8F and ITS1 primers for bacteria and fungi, respectively, and amplicons were sequenced at Eurofins Genomics.

DNA extraction, library preparation and genome sequencing of selected isolates

To characterize the metabolic potential of the defined community employed in the in vitro fermentation experiment, we sequenced the full genomes of the selected isolates. This allowed us to characterize the individual metabolic capabilities of each isolate and, in turn, understand their collective metabolic capabilities. First, we extracted the total DNA from the isolates. To accomplish this, each isolate was cultivated in its respective medium. This involved inoculating 3 ml of the medium with a single colony and incubating under optimal conditions until saturation. The cells were collected through centrifugation at 4,000 × g for 10 min, and the supernatant was discarded. The collected cells were resuspended in 3 ml 10 mM MgCl₂, subjected to centrifugation as before and finally resuspended in 1 ml 10 mM MgCl₂. The suspension was transferred to a 2-ml tube with glass beads (150–212 μm and 425–600 μm in size) and centrifuged at 13,000 × g for 2 min. The supernatant was discarded and 1 ml of DNA extraction buffer (50 mM Tris-HCl pH 8.0, 5 mM EDTA pH 8.0, 350 mM sorbitol, 1% N-lauryl sarcosine, 71 mM NaCl, 0.1% CTAB) with 1 µl of Monarch RNase A (NEB, T3018L) was added. The cells were lysed in a Qiagen TissueLyser II Bead Mill (QIAGEN), which involved shaking at 30 Hz for 10 min. Subsequently, the sample was incubated at 60 °C for 20 min; then an equal volume of chloroform (Scientific Laboratory Supplies, CHE1574) was added. The sample was mixed by inverting several times and centrifuged at 13,000 × g for 5 min. The aqueous layer (top layer) was transferred to a new tube and an equal volume of ice-cold isopropanol was added. The tube was inverted several times and incubated at −20 °C overnight. Afterwards, the tube was centrifuged at 13,000 × g for 5 min, the supernatant was discarded, and the tube was inverted on tissue paper to allow the DNA pellet to air dry for 10 min. Following this, the DNA was resuspended in 50 µl of Milli-Q water. The suspension was centrifuged at 13,000 × g for 5 min, and the supernatant containing the DNA was transferred to a 1.5-ml tube and quantified using a Qubit fluorometer (Thermo Fisher). For the preparation of DNA libraries, the DNA was digested with NEBNext dsDNA fragmentase (NEB, M0348L) in the following reaction mix: 200 ng of DNA in 16 µl Milli-Q water, 2 µl 10× fragmentase reaction buffer v.2 and 2 µl NEBNext dsDNA fragmentase. The reactions were incubated at 37 °C for 20 min, and the process was halted by adding 5 µl 0.5 M EDTA pH 8. The volume was adjusted to 50 µl with Milli-Q water, and DNA fragments between 300 and 500 bp were selectively isolated using double-sided DNA selection with Agencourt AMPure XP beads (Beckman Coulter, A63881). Subsequently, the fragments were end repaired using a mixture comprising 30 µl DNA, 2.5 µl 3 U µl⁻¹ T4 DNA polymerase, 0.5 µl 5 U µl⁻¹ Klenow DNA polymerase, 2.5 µl 10 U µl⁻¹ T4 PNK, 5 µl 10× T4 DNA ligase buffer with 10 mM ATP, 0.8 µl 25 mM dNTP mix and 8.7 µl Milli-Q water. After incubation at 20 °C for 30 min, the fragments were purified again using Agencourt AMPure XP beads. Following this, the DNA fragments were adenylated in a mix containing 34 µl of the end-repaired DNA, 3 µl 5 U µl⁻¹ Klenow exo-, 5 µl 10× Enzymatics Blue buffer, 1 µl 10 mM dATP and 9 µl Milli-Q water. The mixture was incubated at 37 °C for 30 min, followed by 70 °C for 5 min, and then purified using Agencourt AMPure XP beads. Individual samples were indexed through ligation using a mix comprising 10.25 µl DNA, 1 µl 600 U µl⁻¹ T4 DNA ligase, 12.5 µl of 2× Rapid Ligation buffer and 1.25 µl 2.5 µM indexing adapter from the KAPA Dual-Indexed Adapter kit (Kapa Biosystems, KK8722). Samples were incubated at 25 °C for 15 min; then 5 µl 0.5 M EDTA pH 8 was added. The fragments were purified twice with Agencourt AMPure XP beads and then enriched in the following reaction: 20 µl DNA, 25 µl of 2× KAPA HiFi HS Mix (Kapa Biosystems, KK2602), 2.5 µl 5 μM I5 primer (5′-AATGATACGGCGACCACCGAGATCTACAC-3′) and 2.5 µl 5 μM I7 primer (5′-CAAGCAGAAGACGGCATACGAGAT-3′). The reactions were initially heated to 98 °C for 45 s, followed by 14 cycles of 98 °C for 15 s, 60 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 1 min. The resulting DNA libraries were purified using Agencourt AMPure XP beads, quantified on a Qubit fluorometer (Thermo Fisher), and equimolar quantities of individual barcoded DNA libraries were pooled and sequenced (PE150 bp) on an MGI Tech MGISEQ-2000 sequencing platform at Beijing Genomics Institute.

Preparation of a defined community inoculum

A glycerol stock sample of each isolate was plated on the respective selective agar media and incubated at 30 °C for 72–96 h. Subsequently, a 15-ml tube containing 3 ml of selective medium was inoculated with a single colony from the agar plate. The tube underwent incubation at 30 °C with agitation at 80 r.p.m. for fungi, or 200 r.p.m. for bacteria, in a shaking incubator for 72–120 h. Following this incubation period, cells were collected by centrifugation at 4,000 × g at 4 °C for 8 min and subjected to three washes with 10 mM MgCl₂ to remove residual media and cellular debris. The cells were then resuspended in 10 mM MgCl₂, and the optical density at 600 nm (OD₆₀₀) was measured to estimate cell concentrations. A pooled inoculum containing all isolates of the defined community was prepared, with the final concentration of each isolate in the pool set at 10⁹ colony-forming units per millilitre (c.f.u.s ml⁻¹) assuming that 1 OD₆₀₀ unit is equal to 10⁹ c.f.u.s ml⁻¹. In addition, individual strains were systematically removed (single-strain dropout) from the 9-member microbial community to evaluate how the absence of each strain impacts overall community structure and function.

In vitro fermentation set-up and sampling

Mature, ripe, healthy and undamaged cacao pods were harvested from the ICGT. The pods were thoroughly washed with water to dislodge any debris and surface sterilized using 20% hypochlorite solution containing 0.05% Triton X-100 for 10 min in a sterile hood. Afterwards, the pods were rinsed with sterile water and swabbed with 70% alcohol. The pods were carefully opened and the beans with pulp were extracted, maintaining sterile conditions. All beans from all the pods were pooled, mixed to homogeneity and divided into six sterile microboxes (Sac O2, TP1600 + TPD1600), each containing ~1 kg of beans. Three of the microboxes received inoculation with 100 µl of the defined synthetic community inoculum (SYNCOM), while the remaining three microboxes were inoculated with only 100 µl 10 mM MgCl₂, serving as the No SYNCOM control. The beans were incubated for 96 h in a temperature-controlled incubator: 0–48 h at 30 °C, 48–72 h at 35 °C, and 72–96 h at 45 °C. Daily pH measurements were taken by collecting three beans from each fermentation. The testa/pulp were separated from the cotyledons, macerated in 10 ml distilled water using a mortar and pestle, and the pH of the suspensions was determined. For monitoring the microbial community during fermentation, swab samples of the beans were collected at 0, 48 and 96 h using a Zymo Collection Swab (R1104). Samples were collected in duplicate for each fermentation box and placed in Zymo DNA/RNA Shield Lysis and Collection tubes (Zymo, R1104). The tube contents were vigorously shaken for 10 s and the tubes were stored at −20 °C until used. Following the fermentation, the beans were spread on foil trays and placed in an oven at 35 °C for 5 days to reduce the moisture content to <7%. The beans were stirred on drying days 1, 2 and 3 to prevent bean clumping.

For the single-strain dropout experiments, 125 g of dried unfermented cocoa beans was sterilized with a 70% ethanol solution containing 1% Tween 20 and then rehydrated with sterile water. Excess water was discarded and the beans were transferred to a sterile microbox (Sac O2, TP1600 + TPD1600) containing 200 ml of sterile artificial pulp. The artificial pulp contained the following components: 0.14% (w/v) high-viscosity carboxymethyl cellulose (Sigma-Aldrich, C5678-1KG), 0.77% (w/v) low-viscosity carboxymethyl cellulose (Sigma-Aldrich, C5013-1KG), 1.09% (w/v) pectin (Sigma-Aldrich, P9135-500G), 2.5% (w/v) sucrose (Sigma-Aldrich, S0389-1KG), 4% (w/v) glucose (Sigma-Aldrich, G7021-5KG), 5% (w/v) fructose (Sigma-Aldrich, F0127-1KG), 1% (w/v) citric acid (Sigma-Aldrich, C0759-1KG), 0.5% (w/v) yeast extract (Millipore) (Sigma-Aldrich, 70161-500G), 0.5% (w/v) peptone (Millipore) (Sigma-Aldrich, 91249-500G), 0.1% (w/v) calcium lactate pentahydrate (Sigma-Aldrich, C8356-250G), 0.1% (v/v) Tween 80 (Sigma-Aldrich, P4780-100ML), 0.05% (w/v) magnesium sulfate heptahydrate (Sigma-Aldrich, M2773-1KG) and 0.02% (w/v) manganese sulfate monohydrate (Sigma-Aldrich, M7899-500G), adjusted to pH 3.6. The experimental design included a full synthetic microbial consortium (T1, 9-member SYNCOM) and modified versions where individual strains were removed (T2–T10). Control groups consisted of non-inoculated beans (T11 and T12) and beans inoculated with a randomly selected 9-member microbial consortium (T13). For each treatment, four independent fermentations were performed. The beans were incubated for 120 h in a temperature-controlled incubator under the following conditions: 0–48 h at 30 °C, 48–72 h at 35 °C, and 72–120 h at 45 °C. However, for T12, the beans were maintained at a constant temperature of 30 °C for the entire 120-h incubation period. pH measurements of the testa/pulp and cotyledons were recorded daily from a single bean of each fermentation. Swab samples for microbial community analysis were collected at 0, 24 and 48 h (n = 156), while 5 beans from each treatment replicate were sampled at 0, 48 and 120 h for metabolomic analysis (n = 156). After fermentation, beans from each treatment replicate (n = 52) were oven dried at 35 °C, as previously described, to produce cocoa liquors. In addition, total cell counts in the 9-member SYNCOM were measured using a BlauBrand Thoma counting chamber (Brand). To assess potential growth limitations, whether due to nutrient deficiencies in the pulp or environmental stresses from temperature and/or pH, each isolate was cultured in artificial pulp at pH values of 3.6, 4.6, 5.6 and 6.6. Cultures were incubated under three fermentation temperature conditions (30 °C, 35 °C and 45 °C) for 70 h with 200 r.p.m. agitation. Growth was monitored across all conditions by measuring OD₆₀₀.

DNA extraction from in vitro fermentation samples

The DNA samples were placed in a Qiagen TissueLyser II Bead Mill (QIAGEN) and homogenized at 30 Hz for 10 min. Following this, DNA was extracted using the ZymoBIOMICS DNA Miniprep kit (Zymo, D4300) following manufacturer instructions, and the resulting DNA concentration was determined using a Qubit fluorometer (Thermo Fisher).

Bacteria 16S rRNA library preparation and sequencing

We amplified the V3–V4 highly variable region (~480 bp) of the bacterial 16S rRNA gene using the 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) universal primer sequences. Unique frameshifting tags were added to the 5′ end of both primers following the method outlined in ref. ⁹⁷ to enhance library diversity and enable efficient multiplexing of samples for sequencing. Each sample was amplified in triplicate, and for each 96-well PCR plate of reactions, three unique sets of frameshifting tag combinations were employed with both the forward and reverse primers. This approach facilitated the effective multiplexing of samples for sequencing across multiple plates. The reaction mix for each sample included 1 µl DNA, 5 µl 2× KAPA HiFi HS Mix (Kapa Biosystems, KK2602), 0.25 µl of 338F forward primer frameshift mix (10 pmol µl⁻¹), 0.25 µl of 806R reverse primer frameshift mix (10 pmol µl⁻¹) and 3.5 µl Milli-Q water. The amplification protocol involved an initial heating step at 94 °C for 3 min, followed by 24 cycles of 94 °C for 30 s, 50 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 10 min. The PCR products from the triplicate reactions were combined and purified using Agencourt AMPure XP beads (Beckman Coulter, A63881). Subsequently, the PCR products were indexed using 96 unique reverse indexing primers. The indexing mix for each sample included 4.5 µl PCR product DNA, 5 µl 2× KAPA HiFi HS Mix (Kapa Biosystems, KK2602), 0.25 µl forward enrichment primer (10 pmol µl⁻¹) and 0.25 µl reverse enrichment-indexing primer (10 pmol µl⁻¹). The forward enrichment primer used was (5′-AATGATACGGCGACCACCGAGATCTACACGCCTCCCTCGCGCCATCAGAGATGTG-3′), and the reverse enrichment-indexing primer was the TruSeq Read 2-annealing reverse Illumina adapter compatible with the Illumina MiSeq platform. The indexing procedure involved an initial heating step at 94 °C for 3 min, followed by 9 cycles of 94 °C for 30 s, 60 °C for 30 s, and 72 °C for 30 s, with a final extension at 72 °C for 10 min. The DNA libraries were purified using Agencourt AMPure XP beads (Beckman Coulter, A63881) and quantified with a Qubit fluorometer (Thermo Fisher). Subsequently, the libraries were pooled in equal amounts and diluted to 10 pM for sequencing. The sequencing process (PE300) was conducted on an Illumina MiSeq instrument using the Reagent Kit V3 600-cycle (Illumina) at the DeepSeq Sequencing Facility at the University of Nottingham.

Fungi ITS library preparation and sequencing

For fungal profiling, the ITS2 region was amplified and sequenced using the ITS3-F (5′-GCATCGATGAAGAACGCAGC-3′) and ITS4-R (5′-TCCTCCGCTTATTGATATGC-3′) universal primer sequences described in ref. ⁹⁸. To enhance library diversity and facilitate sample multiplexing for sequencing, unique frameshifting tags were incorporated at the 5′ end of both primers using the methodology outlined in ref. ⁹⁷. Every sample was subjected to triplicate amplification. In addition, for each 96-well PCR plate of reactions, three distinct sets of unique frameshifting tag combinations of the forward and reverse primers were used to enable the multiplexing of samples for sequencing from multiple plates. The reaction mix for each sample comprised 1 µl DNA, 5 µl 2× KAPA HiFi HS Mix (Kapa Biosystems, KK2602), 0.25 µl ITS3-F forward primer frameshift mix (10 pmol µl⁻¹), 0.25 µl ITS4-R reverse primer frameshift mix (10 pmol µl⁻¹) and 3.5 µl Milli-Q water. The amplification protocol began with an initial heating step at 94 °C for 3 min, followed by 24 cycles of 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, concluding with a final extension at 72 °C for 10 min. Amplicons from the triplicate reactions were consolidated and purified using Agencourt AMPure XP beads (Beckman Coulter, A63881). Subsequently, the samples were indexed using the TruSeq Read 2-annealing reverse Illumina adapter, pooled as previously described (see ‘Bacteria 16S rRNA library preparation and sequencing’) and sequenced on an Illumina MiSeq platform using the 600-cycle V3 Reagent kit (Illumina) at the University of Nottingham’s DeepSeq Sequencing Facility.

Processing of cocoa beans and liquors for sensory analysis and metabolomics

For the validation of the defined community experiment, dried beans from the three 9-member SYNCOM inoculated fermentations were combined into a single pool, while the three No SYNCOM control batches were pooled separately. The pooled beans were then processed into cocoa liquors and subjected to sensory evaluation as described previously (see ‘Drying, roasting and sensory evaluation of liquor samples’), and well as to gas chromatography–mass spectrometry (GC–MS) analysis. For the single-strain dropout experiment, cocoa beans collected at 0, 48 and 120 h of fermentation were freeze dried for 72 h and ground into a fine powder for liquid chromatography–mass spectrometry (LC–MS) and GC–MS analysis. In addition, the dried fermented beans were processed into cocoa liquors and analysed using sensory evaluation and GC–MS.

Characterization of volatile compounds in cocoa beans and liquors

The aroma and other volatile compounds in the cocoa bean and liquor samples were analysed at the International Flavour Research Centre (IFRC) at the University of Nottingham. Milled powder (1 g) of the cocoa samples was mixed with 10 µl 3-heptanone (0.01 μg µl⁻¹) internal standard in hermetically sealed 20-ml vials and incubated for 5 min at 50 °C in a thermostatic agitator. A 50/30 μm DVB/CAR/PDMS SPME Fibre (Supelco) was used to extract volatile compounds from the headspace of each sample. The SPME fibre extracted for 15 min at 50 °C and desorbed for 0.5 min at 240 °C. The volatiles were analysed by GC–MS using splitless injection into a TRACE 1300 series gas chromatograph coupled with a single quadrupole mass spectrometer (Thermo Fisher). A ZB-WAX-plus column of 30 m length, 0.25 mm internal diameter and 0.250 μm film thickness (Phenomenex) was used with the following time–temperature programme: 40 °C for 2 min, followed by a temperature increase from 40–240 °C at a rate of 6 °C min⁻¹, and then held at 240 °C for 5 min. A minimum of 3 replicates per liquor sample were analysed with randomized sample injections for the validation of the defined community experiment, including 9-member SYNCOM-inoculated samples, No SYNCOM samples, Santander, Huila, Antioquia and reference liquors. For the single-strain dropout experiment, individual biological replicates were utilized. The SPME fibre was conditioned for 3 min at 240 °C between samples. The quality of the headspace GC–MS runs was assessed by running the internal standard after 5–20 consecutive sample runs and estimating the variations in retention time and peak areas. Volatile compounds were identified by comparing each mass spectrum with either the spectra from standard compounds or with spectra in reference libraries (NIST/EPA/NIH Mass Spectral Library). The relative abundance of volatiles was calculated from GC peak areas by comparison with the peak area of the internal standard.

Characterization of non-volatile compounds in cocoa beans

A 100-mg portion of the powdered sample was weighed and placed into a 1.5-ml microcentrifuge tube. The sample was defatted by adding 800 µl ice-cold hexane, vortexing and incubating in a sonic water bath. The mixture was centrifuged and the supernatant was discarded. This defatting process was repeated twice more, and the defatted pellet was dried using nitrogen gas with a sample concentrator (Techne). To extract metabolites, 460 µl 80% methanol was added to the dried pellet, and the mixture was vortexed, sonicated and centrifuged. The supernatant containing the metabolites was saved and the extraction was repeated twice more. The combined supernatant was centrifuged again and the final supernatant was transferred to a new tube. The extract was dried using a Savant SpeedVac SPD140DDA vacuum concentrator (Thermo Scientific) and then stored at −20 °C. Before analysis, the dried pellet was reconstituted with 50% aqueous ethanol, followed by sonication. The sample was centrifuged and the supernatant was transferred to a new 1.5-ml tube. For each time point (0, 48 and 120 h), aliquots from biological replicates of each treatment were pooled separately, transferred to LC–MS vials, capped and stored for further analysis. The generation of the untargeted metabolic profiles was performed using an Agilent 1260 Infinity II Ultra High-Performance Liquid Chromatography system coupled to an Agilent 6546 tandem quadrupole time-of-flight mass spectrometer (Agilent Technologies). Chromatographic separation was performed with an Acuity UPLC HSS T3 column (2.1 × 100 mm, 1.8 μm; Waters) fitted with a KrudCatcher pre-filter (Phenomenex). The flow rate of the mobile phase (A: 5% acetonitrile, versus B: 95% acetonitrile, both with 0.1% formic acid v/v) was at 0.3 ml min⁻¹, with the analytical gradient starting at 5% solvent B, increasing to 15%, 25%, 35%, 45% and 65% at 2, 4, 8, 10 and 12 min, respectively, followed by column washing and re-equilibration (total run time 22 min). Quality control (QC) samples were made by pooling samples from all treatments. After injection of ×10 QC samples to condition the system, each sample was randomized to ×5 injections across the batch. QC samples were injected after every 10 consecutive runs to assess system performance across the batch. Data were collected in MS1 mode scanning 50–1,700 m/z. Reference masses were continuously injected for mass correction.

Analysis of 16S rRNA and ITS regions from Sanger sequencing for isolate identification

To identify the bacterial and fungal isolates cultured, we sequenced the 16S rRNA gene from the bacteria and the ITS region from the fungal isolates. Initially, low-quality bases were trimmed from the sequences, and the results were searched on the National Centre for Biotechnology Information (NCBI) nucleotide database using the Basic Local Alignment Search Tool (BLAST v.2.12.0)⁹⁹ to determine the taxonomy of the species. We performed multiple sequence alignments with the 16S rRNA and ITS sequences for the bacteria and fungi, respectively, using DECIPHER v.2.24.0. The alignments were trimmed with microseq v.2.1.6, transformed into distance matrices with seqinr v.4.2.16, and neighbour-joining trees were constructed with ape v.5.6.2. The resulting trees, based on the bacterial 16S rRNA sequences (Extended Data Fig. 3b) and fungi ITS sequences (Extended Data Fig. 3c), were visualized using ggtree v.3.8.0 with ggtreeExtra v.1.10.0. To assess the overall relative abundance of our collection representing the cocoa fermentation microbiome, we used the tax_glom() function within phyloseq v.1.44.0 to aggregate taxa from the microbiome dataset of the three Colombian fermentations to the family level. Subsequently, we computed the mean relative abundance of each family at each fermentation time point and plotted the results using ggplot2 v.3.4.2.

Genome assembly, annotation and construction of the isolates metabolic network

We used Cutadapt (v.4.6)¹⁰⁰ to eliminate primer and barcode sequences, as well as low-quality sequences, from the paired-end reads of the sequenced genomes of the isolates. Subsequently, the high-quality filtered reads were de novo assembled into a draft genome for each isolate using SPAdes (v.3.15.5)¹⁰¹ with default parameters. The assembled genomes were evaluated for contiguity and completeness using BUSCO (v.5.6.1)¹⁰². Open reading frames in the genomes were predicted with FragGeneScanRs (v.1.1.0)⁸⁵ with default settings. Functional annotation of the predicted proteins was carried out using the eggNOG-mapper (v.2.1.9)⁸⁶ pipeline, utilizing the eggNOG v.5.0.2 database⁸⁷ with Diamond (v.2.0.11)⁸⁸ and MMseqs2 release 12-113e3 (ref. ⁸⁹). The predicted coding DNA and translated amino acid sequences, along with the predicted functional annotations, were combined into genbank-formatted files using emapper2gbk (v.0.3.0)⁹² in ‘genes mode’. Subsequently, the Metage2Metabo (v.1.5.3)^92,93 pipeline was used to generate the metabolic network of the isolates used in the defined community. Metabolites reachable by each of the isolates in the network, based on the cocoa pulp metabolites as the precursor, were compiled into a data matrix and visualized as described previously (see ‘Metabolic network analysis and identification of a defined microbial community’).

16S rRNA and ITS amplicon sequence data processing

Raw reads were demultiplexed and trimmed with Cutadapt (v.4.6)¹⁰⁰. Subsequently, the processed sequences were denoised and collapsed into amplicon sequence variants (ASVs) using the DADA2 v.1.24.0 pipeline. In brief, paired reads were filtered by removing sequences with uncalled bases, eliminating reads with >2 expected errors, and truncating reads when the average quality score dropped to <2. Error rates for forward and reverse reads were determined separately through the learnErrors() function. These error rates were then utilized to infer ASVs individually for both the forward and reverse reads, and the forward and reverse sequences were subsequently merged. The merged ASVs were used to construct an ASV sequencing table, and chimaeras were removed. Bacteria ASVs were classified using the SILVA 138 database¹⁰³, while the fungi ASVs were classified using the UNITE v.9 database¹⁰⁴. Functions in phyloseq v.1.44.0 with microbiome v.1.22.0 and microbiomeutilities v.1.0.17 were used to filter the dataset and remove samples with low read depth, remove unidentified taxa and singletons, transform abundance values using rarefaction, subset and merge sample and taxonomic groups, and perform other dataframe manipulations. To assess alpha diversity across the samples, we calculated the Shannon diversity index using phyloseq v.1.44.0. We used ANOVA to test for significant differences in Shannon diversity indices between groups, and means were separated using Tukey’s HSD test in the agricolae v.1.3.5 R package. For beta diversity, Bray–Curtis dissimilarity matrices were calculated using the phyloseq v.1.44.0 ‘bray’ method, and the variances explained by treatment and fermentation time were estimated by performing PERMANOVA using the adonis2() function in the vegan v.2.6.4 R package. Constrained ordination of beta-diversity was plotted using CAP on the basis of Bray–Curtis dissimilarity matrices calculated with vegan v.2.6.4. We visualized differences in treatment and time with the CAP analysis, using the following models:

PCoA based on Bray–Curtis dissimilarities was used to visualize shifts in microbial community composition across fermentation treatments (T1–T13). Bar plots showed the Euclidean distance between each treatment centroid and T1 (the full synthetic community), reflecting the degree of dissimilarity from the baseline. The relative abundance of taxa was plotted as a stacked bar representation using phyloseq v.1.44.0. The tax_glom() function in phyloseq v.1.44.0 was used to agglomerate taxa, and the aggregate_rare() function in microbiome v.1.22.0 was used to aggregate rare groups.

Analysis of volatile compounds in cocoa beans and liquors

The data were first pre-processed, followed by analysis in R. PCA was performed using Euclidean distances with the prcomp() function, while PCoA was conducted using a dissimilarity matrix computed with the Manhattan distance in the dist() function. Classical multidimensional scaling was then applied using the cmdscale() function. The results were visualized using ggplot2 v.3.4.2. A bar plot of Manhattan distances between each treatment centroid and T1 (the full synthetic community) was generated to quantify dissimilarity from the baseline. To identify enriched volatile compounds among samples, we employed DESeq2 (v.1.40.0)⁶⁵, fitting a GLM with the following design:

From the fitted model, we extracted key comparisons. A volatile compound was deemed significant if it exhibited an FDR-adjusted P value (q value) < 0.05. The GLM analysis results were visualized in a heat map, with the colours representing the log₂ fold change generated by the GLM. Black squares were used to highlight significant differences (q value < 0.05) with log₂ fold change > ±2 between the aforementioned comparisons.

Analysis of non-volatile compounds in cocoa beans

For data analysis, the total ion chromatograms of repeat QC injections were visually assessed to check comparability of runs throughout the dataset. Initially, global MS1 features (peak height >20,000) were first extracted using Mass Profiler software (MP; v.10 Agilent Technologies) and exported to a common .CEF file for each polarity. Thereafter, files for each replicate group were time aligned to the central QC sample in Profinder software (v.10, Agilent); then features were extracted (peak height >5,000) in ‘batch targeted’ mode using the global MS1 features .CEF file as a reference library. Following this pre-processing, PCoA was conducted using Euclidean distances in vegan v.2.6.4. The dissimilarity matrix was computed with the vegdist() function and then used to perform PCoA. The results were visualized using ggplot2 v.3.4.2. A bar plot depicting the Euclidean distances between each treatment centroid and T1 (the full synthetic community) was created to measure dissimilarity from the baseline. We employed DESeq2 (v.1.40.0)⁶⁵ to discern enriched compounds among samples. This was achieved by fitting a GLM with the design:

From the fitted model, we identified key comparisons between treatments and fermentation times, highlighting significant differences where q values were <0.05 and log₂ fold changes exceeded ±2. The results of the GLM analysis were visualized in a heat map, with colours representing the log₂ fold change values generated by the model.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.