Robust experimental strategy. We implemented a well-plate based workflow suited to the analysis of a large number of C. albicans isolates from different clades (Fig. 1, supplemental Table 1). A trial run was performed to assess the feasibility of this high throughput approach, in terms of cell growth and possible cross contamination. From the results of the trial run (Fig. 1b), we observed low cell densities for some isolates at the end of incubation, therefore we decided to increase the inoculum optical density (OD measured at 600 nm) from 0.1 to 0.2 in the culture, to obtain higher metabolic signal. Additionally, as the metabolomics studies are prone to batch effect, we selected one isolate in each of the five most represented C. albicans clades (Clades 1, 2, 4, 9, and 11) and in Clade 13. These six representative isolates were selected as they had similar OD at the end of the run and were included in duplicates in all subsequent experiments. Consequently, these isolates could be further used for possible statistical batch effect correction.
With this workflow, we have been able to grow 95 out of 96 isolates. One isolate (CEC5136 from cluster C) did not grow. This translates to more than 600 samples, as the representative isolates were cultured in duplicate in all batches (n = 14) and the other isolates were cultured in duplicate in three successive culture batches (n = 6). The isolates CEC3558 and CEC5120 have grown only in two out of three batch cultures (n = 4). Even with different inoculum OD in the final culture, the growth profiles obtained for both trial and experimental cultures were similar (Fig. 1b and 1c), with the Spearman correlation coefficient between OD averages of isolates at ρ = 0.66 (supplemental Fig. 1).
From the cell culture, we separated the supernatant and analyzed the extracellular metabolites through 1H NMR. During the analysis, the culture batch grouping was maintained but the sample analysis sequence was randomized. The spectra were then binned to extract metabolic fingerprints. The quality of data from the generated metabolic fingerprints was assessed through unsupervised analysis (principal component analysis, PCA). The C. albicans samples were firstly plotted along with control media and blanks samples and then, only samples results but labeled as batch runs number (batches 1–7) (Fig. 2). The resulting scatter plots show that samples were neatly categorized by sample type (Fig. 2a) and secondly that no clustering between batches could be discerned (Fig. 2b). These indicate that there was no cross contamination during the culturing and sample preparation and variations do not arise from the technical run. Considering that batch effect was negligible, batch effect correction using the selected representative isolates was not applied.
Using the same metabolic fingerprint dataset, 40 metabolites were annotated and 39 were quantified using Chenomx NMR Suite (version 8.6). Lactate was not quantified due to low selectivity detection, caused by interference from the threonine spectrum. We have annotated and quantified around 65% of the peaks covering the spectral region from 0.2 to 10 ppm (supplemental Fig. 2). With the annotated and quantified dataset, it was also possible to determine the variations arising from either technical or biological parameters. Therefore, using the metadata provided for each isolate (isolate name, clade affiliation, country of origin and infection pathways) and culture specific parameters (final OD and batch run), a variance partitioning analysis (VPA), an analysis designed to handle data generated from high-throughput experiments (31), was performed. Through the VPA, the contribution of each factor (parameter) in driving samples variabilities can be determined. Based on the VPA, there were five main factors that driven the separation among samples (Fig. 3a) and each of these factors have differential influence on each of the quantified metabolite (Fig. 3b). High residual effect indicates that most of the variations were due to information not captured in the metadata. This was also a consequence of comparing closely related samples, as was the case for C. albicans isolates. The VPA also showed that batch effect was negligible, consistent with the results obtained for the binned NMR spectra dataset using PCA analysis. Taking out residual and batch effects, the three remaining sources of variation were all biological factors: Strain (isolate), OD (as a proxy for biomass) and Clade. For other parameters that did not appear prominently in the VPA, such as infection pathways and country of origin, we have run an unsupervised analysis (PCA) on these, with missing values classified as one category. The resulting PCA was consistent with the VPA as no apparent clustering based on these parameters were observed (supplemental Fig. 3).
The effect of OD was also evident in the metabolic fingerprint dataset. A PCA of the metabolic fingerprints dataset with OD as label showed an OD gradient and formation of two distinct Clade 13 clusters (supplemental Fig. 4). To compensate the OD effect, we employed various strategies to normalize the datasets. Particularly for the quantified metabolite concentrations, we have done ∆ concentration normalization (normalized concentration = [isolate metabolite concentration – control metabolite concentration]/OD) and min max normalization (normalized value = metabolite concentration/[isolate OD/lowest isolate OD]). However, comparison of RSD (relative standard deviation) values between the raw concentrations (non-normalized) and the two normalized concentrations showed that these strategies gave poorer results than without normalization (supplemental Table 2). One of the reasons that the normalization did not work was that some metabolites were not utilized by some isolates and induce a bias. When the concentrations of these metabolites fluctuate around the control, normalization increases the standard deviations tremendously. Further, when these isolates have low OD, normalization produced falsely high values. Thus, we used the raw concentrations for the analysis. The OD gradient, including the bifurcation observed for Clade 13, was only evident in the second principal component of the PCA (Supplemental Fig. 4). Further, it should also be noted that bifurcation of Clade 13 isolates based on OD was due to one isolate (CEC4943). For all Clade 13 isolates, clade-based delineation from the other isolates was apparent in the first principal component, hinting that other inherent properties was exerting more effect than the OD.
Clade delineation. With metabolic fingerprinting and quantified metabolite datasets, we ran a supervised analysis (partial least square discriminant analysis, PLS DA) with clade as the response (grouping) variable (Fig. 4). While the cumulative explained variance (R2Y) as well as the model predictive accuracy (Q2Y) for metabolic fingerprints (0.34 and 0.23) and quantified metabolite (0.25 and 0.17) datasets are low, these models were still able to distinguish Clade 13 from all other clades. The similarity between the metabolic fingerprinting and quantified metabolite datasets reflects that the quantified metabolites capture the high variation within samples, although the reduced R2Y values in the metabolite quantified dataset (Fig. 4b) suggest that the non-quantified features in the metabolic fingerprint datasets (Fig. 4a) contribute to sample variations as well. Further, such values were expected as variability attributed to clade affiliation was lower than those ascribed to isolate and OD.
To investigate further the possible drivers of clade delineation, the VPA was used to determine which metabolites were most responsive to clade affiliation (Fig. 3b). Based on the VPA, we examined the mean clade concentrations of choline and trehalose and found that these were lower than that of the control (non-inoculated YPD media), except for isolates from Clade 13 (Fig. 5a and 5b). These reduced metabolite concentrations were both significant (Kruskal-Wallis Test: choline, p < 0.001 and trehalose, p < 0.001) and have large effect sizes (choline, ε2ordinal = 0.36 and trehalose, ε2ordinal = 0.35), indicating that clade affiliation was a determinant of metabolite concentration. For Clade 13, the concentration values for trehalose and choline, although slightly lower, were not significantly different from that of the control (pHolm adjusted value = 1.00 for both metabolites). Thus, concentration values of both metabolites in the five Clade 13 isolates were compared with the control. Choline concentrations in Clade 13 were unusual as the values of each isolate varied significantly with that of the control (Fisher’s ANOVA, p < 0.001) and have registered large effect size (ε2ordinal = 0.64) (Fig. 5c). Two of the isolates showed significantly higher concentrations (CEC4104 and CEC5020, pHolm adjusted value < 0.001 for both) while one isolate had significantly lower concentrations (CEC4943, pHolm adjusted value < 0.001). On the other hand, all Clade 13 isolates have trehalose concentrations that were not statistically different with that of the control (Kruskall-Wallis Test, p = 0.13) (Fig. 5d).
Since an instance of clade delineation was observed, clade-based analysis was further pursued through supervised (PLS DA) pairwise comparison of metabolic fingerprints. With this comparison, we could demonstrate that the three largest clades are delineated from each other (Fig. 6). These pairwise comparison show high cumulative explained variance (R2Y values higher than 0.85) and high model predictive accuracy (Q2Y values higher than 0.78). Thus, while clade delineation was dominated by the distinct Clade 13 metabolic phenotype, pairwise clade delineation showed that globally, isolates within a clade share collective metabolic signatures distinguishable from that of the isolates from another clade.
Isolate specific variation. The most prominent biological factor driving variation based on the VPA could be attributed to factors inherent to the isolate. When the effect of this inherent isolate variability was assessed on individual metabolites, it was observed that fumarate and glycerol are the most responsive (Fig. 3b). Both metabolites showed increase in concentrations when compared with the control. These increases were most visible for fumarate, which ranges from 30 to 127 µM, while this metabolite was detected in minute quantities in the control media. Also, a striking feature of the dataset was the pronounced increase in the glycerol production of the clade 2 isolate CEC3558 (12296 µM) which is seven times the levels in control and almost three times more than the mean concentration observed for isolates within Clade 2 (Fig. 7b). The CEC3558 is a commensal isolate originating from Belgium (supplemental Table 1). This isolate survived only in two out of three successive cultures and has shown poor growth (OD600 = 0.4 ± 0.03 (mean ± SD), n = 4 replicates) compared to the other isolates (Fig. 1). Further, the differences in the concentration values of fumarate and glycerol were apparent between isolates (Fig. 7). Notable, however, is the low variations observed within isolate replicates. The RSD is below 20% for fumarate and 30% for glycerol, indicating that culture conditions do not vary much between batch runs. We performed within clade comparison of isolate concentrations for these two metabolites, to ascertain the high isolate specific variations observed in the VPA (supplemental Table 3). For both metabolites, isolate concentrations were all significant in all comparisons. The computed effect sizes were also large for choline (ε2ordinal ranging from 0.27 to 0.78) and very large for trehalose (ε2ordinal ranging from 0.87 to 0.95), indicating that differences in the concentration values are caused by elements inherent to the isolates.
Considering that metabolite production can be a function of biomass, we evaluated if culture OD affects metabolite concentration as OD effect is prominent in the VPA. We performed correlation analysis on both fumarate and glycerol (Supplemental Fig. 5). Although highly significant, the positive correlation between OD and fumarate (Spearman rank correlation, p < 0.001 and ρ = 0.18) and glycerol (p < 0.001 and ρ = 0.22) concentrations, respectively, are very weak. This indicates that while OD affects metabolite concentrations, its contribution is overall negligible.
To determine the relative utilization of the annotated metabolites we constructed a heat map (Fig. 8). From the heat map, it was clear that there was a marked reduction of metabolites used as precursor molecules for the central carbon metabolism. Consequently, there is an increase of metabolites produced as by-products of the glycolysis and as intermediates of the TCA (Tricarboxylic acid)/Krebs cycle. Except for the close proximity of four out of five Clade 13 isolates, these utilization patterns do not correspond to clade nor OD. This strongly suggests that metabolite utilization is isolate-specific.