To examine how TB disease perturbs immune responses, we stimulated whole blood from active TB patients and LTBI controls with relevant stimuli that included TB Ag, BCG and IL-1β, plus a non-stimulated control (Null). Principal component analysis (PCA) of 622 immune genes measured showed clustering of stimuli-specific responses in both LTBI and TB, with 54% of the total variance captured by the first three principal components ( Figure 1A and Supplementary Figure S1A ). As previously reported (16), several genes were differentially expressed in active TB disease compared to LTBI in the absence of stimulation, namely 32% (200 genes, q <0.001) of the immune genes examined.

Focusing on differential genes that were revealed only by immune stimulation, we compared LTBI and TB groups using a t-test with a stringent cut-off value of q<0.001. This resulted in 114 genes (18% of total) with differences between LTBI and TB groups only after specific stimulation (TB Ag; 71, BCG; 28, IL-1β; 15) ( Figure 1B and Supplementary Table S2 ), with TB Ag stimulation revealing significantly more differences than the other stimulations ( Supplementary Table S3 ). Examples of stimuli-specific induced differences for each condition and the Null are illustrated in Figure 1C ; Supplementary Figure S1B and Supplementary Table S4 .

Whole blood transcriptomic approaches, based on signatures that primarily comprise interferon stimulated genes (ISG), are being developed for TB diagnosis, with promising strategies based on a 3-gene signature (Sweeney3 signature) in which a TB score is calculated using the formula (GBP5 + DUSP3)/2 – KLF2 (17), or the 11-gene blood TB risk signature, RISK11 (18). As TB Ag stimulation induced a high number of differences between LTBI and TB patients, we tested whether standardised immune stimulation could improve the performance of these signatures by comparing ROC curve analysis from Null and TB Ag stimulated samples ( Figures 1D, E ; Supplementary Figure S1D , S1E ). For the Sweeney3 signature, TB Ag stimulation accentuated the pre-existing differences at baseline, which translated into a superior area under the receiver operating characteristic curve (AUC) for Null_AUC of 0.82 [95% CI (0.66-0.94)] and for TB Ag_AUC of 0.95 [95% CI (0.88-1)] improving the ability to classify LTBI and TB cases (p=0.02) ( Figure 1E ). Similar results were obtained with the RISK11 signature, with Null_AUC of 0.82 [95% CI (0.69-0.94)] and TB Ag_AUC of 0.97 [95% CI (0.93-1)] (p=0.02) ( Supplementary Figure S1E ). These differences were not present after successful antibiotic treatment of the TB patients in the Sweeney3 signature ( Figure 1F ; Null_AUC = 0.63 [95% CI (0.45-0.82)] and TB Ag_AUC = 0.51 [95% CI (0.32-0.71)], (p=0.4) or the RISK11 signature ( Supplementary Figure S1F ; Null_AUC = 0.57 [95% CI (0.38-0.77)] and TB Ag_AUC = 0.63 [95% CI (0.45-0.81)] (p=0.7), confirming the specificity of these signatures.

In summary, immune stimulation revealed multiple differences between latent M.tb infection and active TB disease that likely reflect relative differences in ISG expression in blood leukocytes, which may help to better understand perturbed immunity during TB pathogenesis and potentially improve new diagnostic strategies.

The TB Ag stimulation illustrated in Figures 1D, E captures differences in antigen-specific responses of CD4+ T cells, and subsequent downstream responses, between LTBI and TB patients. To explore the broader impact of TB disease in mounting efficient immune responses to complex stimuli, we investigated the response to BCG. To identify immune pathway differences, we performed correlation analysis between BCG-induced cytokines at the protein level (Fold-change>1.3 in the BCG condition respect to the Null) in the LTBI and TB groups ( Figure 2A ). This revealed negative correlations in the control LTBI group, that were absent or markedly weaker in TB patients, suggesting altered regulation of cytokine production in TB disease. This included several key cytokines important for immune responses against TB such as IL-1α, IL-1β, IL-10, GM-CSF (CSF2) and TNF, which all showed a positive correlation with the IL-1 receptor antagonist (IL-1ra) in TB patients ( Figure 2A right panel), in contrast to a negative correlation in LTBI ( Figure 2A left panel).

Because IL-1 signalling is essential for M.tb control (20) and three key members of the IL-1 family (IL-1α, IL-1β and IL-1ra) showed significant differences, we focused on understanding the distinct secretion of IL-1 family cytokines between LTBI and TB. We started by investigating levels of IL-1ra and observed that TB patients presented higher concentrations of IL-1ra both at baseline and after BCG stimulation ( Figure 2B ). Here, we explore underlying mechanisms of differential IL-1ra secretion at baseline ( Figures 2 , 3 ) and after BCG stimulation ( Figure 4 ). To better understand the drivers of IL-1ra secretion at baseline in TB, we performed correlation analysis of the protein data comprising 40 cytokines/chemokines in the Null tube and showed that the strongest association (r=0.52, p=0.0001) was with IL-8 ( Figure 2C ), a cytokine which also had significantly higher levels in TB patients ( Figure 2D ). Expression levels of both IL-8 receptors (CXCR1 and CXCR2) were also higher in TB patients ( Figure 2E ), however stimulation of blood from healthy controls with IL-8 did not induce IL-1ra ( Figure 2F ). We next explored gene expression signatures commonly associated with these two cytokines and found genes related to the peroxisome proliferator-activated receptor (PPAR)γ pathway ( Figure 2G ). We validated our findings using a published dataset of transcriptional analysis of human monocyte-derived macrophages (hMDMs) infected with different bacterial strains (19). This showed a weak, yet significant positive correlation (r=0.34, q=0.03) between IL1RN (the transcript for IL-1ra) and PPARG transcripts in hMDMs infected with mycobacteria, but not with non-mycobacterial species (r= - 0.28, q=0.1) ( Figure 2H ). These results suggest an immuno-metabolic link between the PPARγ pathway and the receptor antagonist IL-1ra.

PPARγ is a ligand-activated transcription factor that plays an essential role in metabolism and energy homeostasis, particularly in adipogenesis (21, 22). To explore potential metabolic pathways involved in IL-1ra secretion, we performed unsupervised mass spectrometry analysis on supernatants from unstimulated LTBI and TB patients’ blood and examined which metabolites correlated with IL-1ra protein at baseline in the TB group. Interestingly, 5 of the top 6 ranked metabolites showing significant correlations (r>0.6, q ≤ 0.05) with IL-1ra in TB patients belonged to pregnane steroid derivatives ( Figure 3A and Supplementary Table S5 ). These correlations were observed in TB patients, and not in LTBI ( Figures 3A, B ). Interestingly, levels of these steroids did not differ between the LTBI and TB groups ( Supplementary Figure S2A ), in contrast to PPARG expression, which was higher in TB ( Supplementary Figure S2B ). Pregnane steroids have been shown to be potential ligands of PPARγ (23), and PPARγ is mainly expressed by CD14+ monocytes (Human Cell Atlas). To investigate whether activation of PPARγ could drive IL-1ra secretion, we incubated CD14+ monocytes isolated from healthy donors, with combinations of the PPARγ agonist rosiglitazone and the PPARγ antagonist GW9662, in the presence of BCG. We validated our approach by measuring surface expression of CD36, a fatty acid translocase that is up-regulated by PPARγ (24–26). CD14+ monocytes stimulated with rosiglitazone upregulated CD36, and this could be prevented with GW9662 ( Supplementary Figure S2C ). We observed that PPARγ activation through rosiglitazone increased IL-1ra secretion by CD14+ monocytes, which could be inhibited by GW9662 ( Figure 3C ). To further confirm the role of PPARγ in IL-1ra secretion, we silenced PPARγ using small interfering RNA (siRNA). We confirmed PPARγ knockdown in CD14+ monocytes by measuring levels of PPARγ protein by flow cytometry ( Supplementary Figure S2D ). Both at baseline and after BCG stimulation, silencing of PPARγ resulted in decreased IL-1ra secretion ( Figure 3D ), with this knockdown impacting both conditions equally ( Figure 3E ). This suggests that, in healthy volunteers, the PPARγ pathway is not involved specifically in BCG-induced IL-1ra, although these data confirm PPARγ as a potential regulator of baseline IL-1ra production in TB disease.

TB patients secreted higher protein levels of IL-1ra both at baseline and after BCG stimulation compared to LTBI ( Figure 2B ). The differences in IL-1ra secretion between TB and LTBI after BCG stimulation could not be explained by PPARγ ( Figure 3E ). Therefore, we investigated the drivers of IL-1ra secretion in the context of BCG stimulation. Previous studies have identified TNF, IL-1 and type I and II IFNs as drivers of IL-1ra secretion (27–32). To test this in our system we stimulated whole blood from healthy individuals with TNF, IL-1β, IFNγ, IFNβ, IFNα, and IL-8 and confirmed TNF and IFNβ (q=0.001) as the strongest inducers of IL-1ra ( Figure 4A ). Furthermore, TNF and type I IFNs acted synergistically in inducing IL-1ra, as the combination of these cytokines (TNF + IFNα/β) resulted in a 4-fold higher induction as compared to stimulation with each cytokine alone ( Figure 4B ). We then investigated the effect of blocking IFN and TNF signalling through the use of blocking antibodies. Blocking of TNFR (q=0.001), but not IFNAR (q=0.7), resulted in significantly less IL-1ra production ( Figure 4C ). To obtain further insight into the potential contribution of TNF and type I IFN in IL-1ra production upon BCG stimulation, we examined a previously published kinetic study of BCG stimulated whole blood (33). This consisted of five healthy donors at 13 time points within 30h of BCG stimulation from which we measured IFNα and IFNβ to analyse with existing TNF and IL-1ra data. A time-dependent linear regression on IL-1ra using IFNα, IFNβ and TNF protein secretion as predictors revealed a strong association with TNF (q=0.002) but not type I IFNs (IFNα q=0.6, IFNβ q=0.6) ( Figure 4D ). IFNβ was not induced upon BCG stimulation, and IFNα was induced by BCG from the 10h time point ( Figure 4D ).

To apply these findings to TB patients, we measured TNF protein after BCG stimulation and observed that levels of TNF were surprisingly higher in LTBI compared to TB ( Figure 4E ). We next examined expression levels of the TNF receptor subunit (TNFRSF1B) and observed significantly higher levels (q=0.0001) in TB patients ( Figure 4F ), which may reflect increased signalling. To test this hypothesis, we calculated a gene score based on unique cytokine-induced gene signatures (see Materials and Methods and Supplementary Table S6 ) as previously described (15). The TNF gene score, which is a surrogate for quantifying signalling by this cytokine, was higher in the TB group compared to LTBI after BCG stimulation ( Figure 4G ), and also positively correlated (Rs=0.32, p=0.03) with IL-1ra protein concentrations ( Figure 4H ). These results suggest a role for TNF in BCG-induced IL-1ra, a response that is augmented in active TB disease.

BCG stimulation revealed dysregulation of multiple arms of the IL-1 pathway in active TB disease ( Figure 2A ). As well as observing differentially induced IL-1ra between LTBI and TB ( Figure 2B ), we also saw significant differences in IL-1α and IL-1β secretion ( Figure 5A ). In contrast to IL-1ra, both IL-1α and IL-1β protein were secreted at significantly higher levels in LTBI upon BCG stimulation. However, these differences were not mirrored at the transcript level ( Figure 5B ). We therefore hypothesized that post-transcriptional/post-translational modifications might be impaired during active TB disease. IL-1β requires post-translational cleavage to signal and engage functional outcomes, whereas IL-1α can manifest basal levels of activity in its unprocessed form (34). As this protein processing can occur through either inflammasome-dependent or independent manners (34, 35) we explored whether both pathways were engaged in our whole blood stimulation system. We stimulated blood from healthy donors with BCG in combination with either the selective NLRP3 inflammasome inhibitor MCC950, the Syk inhibitor R406 (e.g. it inhibits the NLRP3 inflammasome and C-type-lectin receptors [CLRs]), a serine-protease inhibitor (3,4-Dichloroisocoumarin, DCI) capable of blocking a range of enzymes involved in inflammasome-independent cleavage of IL-1 immature proteins, including granzymes and cathepsins. Treatment with R406 and DCI, but not MCC950, resulted in a partial but significant reduction of IL-1β secretion after BCG stimulation ( Figure 5C ), suggesting a role for Syk-proteins outside the NLRP3 inflammasome, and serine-proteases in the secretion of IL-1β in stimulated whole blood cultures.

We then examined differences in these pathways between LTBI and TB patients. The transcripts of components from different inflammasomes (i.e. NLRP3, NLRC4, AIM2) NLRP3, SYK, CARD9 and CASP1 were significantly more expressed after BCG stimulation in TB patients compared to LTBI ( Figure 5D ). Cathepsin G, a neutrophil serine-protease that cleaves pro-IL-1β, was also expressed at higher levels in the TB group ( Figure 5E ), whereas no differences were observed in cathepsin C transcripts. In contrast, other genes that participate in inflammasome-independent IL-1 cleavage, such as granzymes A, B and K, were expressed at higher levels in the LTBI group ( Figure 5F ) suggesting that this pathway was more relevant for BCG induced IL-1β. Concentrations of granzyme B in the BCG supernatants were higher in the LTBI group compared to TB, however, no significant differences were observed with granzyme A ( Figure 5G ). It was notable that most granzyme A levels were below the limit of detection of the assay. Levels of granzyme B strongly correlated with IL-1α (R=0.53, p=0.0001) and there was a moderate but significant (R=0.33, p=0.02) correlation between IL-1β protein and Granzyme A ( Figure 5H ). These data support our hypothesis of aberrant IL-1 induction in TB disease, with lower levels of granzyme protein associating with reduced cleavage of IL-1 precursors, and thus less active forms of IL-1 cytokines.

While BCG stimulation revealed perturbed secretion of IL-1 cytokines in active TB disease ( Figures 2B , 5A ), it did not allow examination of potential differences in IL-1 signalling. To investigate this, we stimulated whole blood from TB patients and LTBI controls with IL-1β and defined a signature of 107 induced genes (Null vs IL-1β q<0.001, Supplementary Table S7 ). Of this IL-1β driven gene signature, 72 genes were differentially expressed (q<0.001) between TB and LTBI pre-treatment ( Figure 6A , left panel), with 95% of the differential genes showing higher levels of expression in TB than in LTBI. Among these differentially expressed genes were key signal transducers and transcription factors such as NKFB1/2 and IRAK2/3 ( Figure 6B ), suggesting altered downstream signalling. Analysis of 18 TB patients after successful antibiotic treatment ( Figure 6A , right panel) showed that the perturbed IL-1β signalling response was completely reset and showed no differences compared to controls ( Figure 6B and Supplementary Table S8 ). Our data shows that IL-1β signalling is strongly perturbed in active TB disease, and this response is restored after successful treatment.

Adult patients with sputum Xpert MTB/RIF-positive TB disease who tested HIV-negative were identified and recruited at clinics in Worcester, South Africa (62). Healthy QuantiFERON-TB Gold (QFT) positive (M.tb-infected) adults were recruited from communities living in or around Worcester and those who matched TB cases by age and gender were enrolled as LTBI controls ( Table 1 ). Blood was collected prior to treatment initiation in active TB cases and again 12-18 months later, after successful completion of treatment defined by clinical cure (V2, n=18). For the QFT+ controls, blood was also collected at a second time-point, 12-18 months after the initial visit (V2, n=19). The TB clinical study protocols and informed consent forms were approved by the Human Research Ethics Committee of the University of Cape Town (ref: 234/2015). Healthy donor blood was obtained from the CoSImmGEn cohort of the Investigation Clinique et Accès aux Ressources Biologiques (ICAReB) platform, (Centre de Recherche Translationnelle, Institut Pasteur, Paris, France) or from the Etablissement francais du sang (EFS, France). Blood was collected in sodium-heparin tubes. Written informed consent was obtained from all study participants.

TruCulture tubes were prepared in batch as previously described (14, 15) with the following stimuli: QFT antigens (ESAT-6, CFP-10, TB7.7), Bacillus Calmette–Guérin (BCG; Sanofi Pasteur, 10⁵ units/ml), and IL-1β (Peprotec, 25 ng/ml). They were resuspended in 2 ml of buffered media and maintained at -20°C until use. Blood was collected in sodium-heparin tubes (50 IU/ml final concentration). Within 30 min of collection, 1 ml of whole blood was distributed into pre-warmed TruCulture tubes, inserted into a dry block incubator, and maintained at 37°C (+/- 1C), room air for 22 hr (+/- 15 min). After incubation, a valve was inserted to separate cells from the supernatant and to stop the stimulation reaction. Upon removal of the liquid supernatant, cell pellets were resuspended in 2 ml Trizol LS (Sigma), vortexed for 2 min, and rested for 10 min at room temperature before -80°C storage. Cytokine kinetic secretion data sets were previously described (33). All stimulations were performed for 22h unless otherwise specified.

For monocyte experiments, blood was collected in buffy coats from healthy donors sampled at the EFS, France. First, PBMCs were isolated by density centrifugation from peripheral blood leukocyte separation medium (Ficoll-Paque™ plus; GE Healthcare). Human monocytes were purified by negative selection with the classical monocyte isolation kit (Miltenyi). Monocytes were cultured at 1x10⁶ cells/mL in RPMI 1640 (Invitrogen, Gaithersburg, MD) (R10) containing 10% heat-inactivated fetal bovine serum, penicillin (100 U/ml), streptomycin (100 μg/ml; 1% pen-strep), and 1 mM glutamine (HyClone, Logan, UT).

Supernatants from TruCulture tubes were thawed on ice and tested by Luminex xMAP technology for a total of 32 proteins including cytokines, chemokines and growth factors as previously described (14). Samples were measured on the Myriad RBM Inc platform (Austen, Texas, US) according to CLIA guidelines (set forth by the USA Clinical and Laboratory Standards Institute). The least detectable dose (LDD) for each assay was derived by averaging the values obtained from 200 runs with the matrix diluent and adding 3 standard deviations to the mean. The lower limit of quantification (LLOQ) was determined based on the standard curve for each assay and is the lowest concentration of an analyte in a sample that can be reliably detected and at which the total error meets CLIA requirements for laboratory accuracy. IL-1β, IL-1ra, IL-1α, IL-6, IP-10 and Granzyme A/B were also measured using Human ELISA Kit (ThermoFischer Scientific) or Human Custom Procarta-plexes (ThermoFischer Scientific), following manufacturer’s instructions.

To detect low concentrations of cytokines homebrew Simoa digital ELISA were used for IFNα (66) and IFNβ (67), as previously described.

Cell pellets in Trizol LS were thawed on ice for 60 min prior to processing. Tubes were vortexed twice for 5 min at 2000 rpm and centrifuged (3000 x g for 5 min at 4°C) to pellet the cellular debris generated during lysis. RNA was isolated from whole blood samples using the NucleoSpin 96 RNA tissue kit protocol (Macherey-Nagel) with some modifications as previously described (15).

RNA integrity was assessed (Agilent RNA kits for the 2100 Bioanalyzer System). The NanoString nCounter system was used for the digital counting of transcripts. RNA was quantified using the Qubit RNA HS Assay Kit (Thermo Fischer Scientific) and 100ng of total RNA hybridized with the Human Immunology v2 (plus 30 additional genes ( Supplementary Table S1 ) relevant for TB that were included) Gene Expression CodeSet according to the manufacturer’s instructions. Samples were processed in 7 batches (4 for Visit 1, 3 for Visit 2), within which the samples were randomized, and the same lot of reagents was used for all samples. All samples were normalized together following background subtraction of the negative control probes, using positive control probes and housekeeping genes (SDHA, HPRT1, POLR2A, RPL19, G6PD, TBP) selected by the GNorm method as previously described (68). This was done using the nSolver™ analysis software (NanoString technologies). Quality control for our data involved checking the following metrics: fields of view counted (flag if < 0.75), binding density (flag if not in the 0.05-2.75 range), linearity of positive controls (flag if R2 < 0.9) and limit of detection for positive controls (flag if 0.5 fM positive control < 2 SDs above the mean of the negative controls).

Supernatants from TruCulture Null tubes were thawed and tested by Ultrahigh Performance Liquid Chromatography-Tandem Mass Spectroscopy (UPLC-MS/MS) for a total of 696 metabolites. Samples were measured on the Metabolon Inc platform (Morrisville, North Carolina, US). For metabolite quantification, peaks were quantified using area-under-the-curve. A data normalization step was performed to correct variation resulting from instrument inter-day tuning differences, and values were rescaled to set the medial equal to 1. Missing values were imputed with the minimum.

The NLRP3-specific inhibitor MCC950, the inhibitor of Syk phosphorylation R406 (In vivogen) and the serine-protease inhibitor 3,4-Dichloroisocoumarin (Sigma Aldrich) were used at 1 μM, 5 μM and 10 μM, respectively. The PPARγ agonist Rosiglitazone (Sigma-Aldrich) and the antagonist GW9662 (Sigma-Aldrich) were used at 10 μM. All drugs were added 1h before BCG stimulation.

Whole blood from healthy donors was pre-incubated with antibodies or a relevant isotype control at 10 µg/mL for 1h (unless otherwise stated), using the following: anti-human CD120a (Clone MABTNFR1-B1, BD Bioscience) antibody, anti-human Interferon alpha/beta receptor 1 antibody (Anifrolumab), anti-human Interferon alpha/beta receptor 2 antibody (PBL, clone MMHAR-2). Whole blood was then diluted 1/3 in TruCulture media and stimulated with BCG or with selected human recombinant cytokines for 22h at 37°C, as follows: IFNα (IntronA; 1000 IU/mL), IFNβ (Betaferon; 1000 IU/mL), IFNγ (Peprotec; 25 ng/mL), IL-1β (Peprotech; 25 ng/mL), TNF (Miltenyi; 10 ng/mL), IL-8 (Biolegend; 25 ng/mL).

Monocytes were seeded at 2x10⁵ cells/200μl in 96-well plates and incubated at 37°C. Cyclophilin B (control) and PPARγ siRNA (SMARTPool, Dharmarcon) were diluted in DOTAP (1,2-dioleoyl-3-trimethylammonium-propane; Roche Applied Sciences). The mix was gently mixed and incubated at room temperature for 15 min. After incubation, the mix was added to cells in culture at a final concentration of 160 nM. Cells were then incubated at 37°C for 24 hours before adding BCG stimulation for 16 hours.

Cells were washed in PBS and resuspended in PBS containing 2% fetal calf serum and 2 mM EDTA and stained with the extracellular mix using anti-human CD14-BV421 antibody (clone M5E2, BD Bioscience) and anti-human CD36-BUV605 antibody (clone CD38, BD Bioscience) at 1:200. For PPARγ intracellular staining, a Fixation/Permeabilization Solution Kit (BD Cytofix/Cytoperm) was used according to the manufacturer’s protocol. Briefly, the cells were fixed for 10 min at 4°C with 100 µL of the Fixation/Permeabilization solution and then washed and stained in 100 µL of the BD Perm/Wash Buffer containing anti-human polyclonal PPARγ antibody (PA3-821A, Invitrogen) at 1:25 for 1 hour at 4°C. Goat anti-rabbit IgG-Alexa Fluor 700 was used as a secondary antibody at 1:1000 (ThermoFisher Scientific). Data acquisition was performed on a FACS LSR flow cytometer using FACSDiva software (BD Biosciences, San Jose, CA). FlowJo software (Treestar, Ashland, OR) was used to analyze data.

ANOVA (Kruskal-Wallis) testing was performed for multi-group comparisons and paired t-tests for intra-patient analysis with Qlucore Omics Explorer, v.3.2 (Qlucore) or GraphPAD Prism. To correct for multiple testing we report false discovery rate (FDR)-adjusted ANOVA p values; q values. Dot plot graphs were compiled with GraphPad Prism v.6.0, heatmaps and Principal Component Analysis (PCA) plots with Qlucore. ROC curves were calculated with R (v.3.4.4) and results drawn with graphical package ggplot2 (v.3.1.0). UpSet plots were drawn with UpSetR (v1.3.3), and correlations plots corrplot (v0.84).

Time series analysis on the protein secretion measurements, from 5 different donors, was conducted using a linear mixed-model approach. The time dependency was modelled by incorporating the parameter time as a continuous linear predictor, alongside the other protein predictors, and the donors were modelled as a random effect. The model was implemented using the R package ‘lme’’ (v1.1-20).

For calculation of cytokine-specific genes score, the genes uniquely induced by each cytokine were identified from our previous study in healthy donors (15). The first Principal component of the expression matrix was computed using the PCA() function of the FactoMineR package (version 1.42) with the options “scale” set to TRUE, and the option “ncp” set to 1. Position of each sample on the first component was then extracted and used for further analysis (R v.3.6.1).

Pearson correlations were performed in R (v3.6.0) between IL-1Ra protein concentrations quantified by Luminex and the metabolites quantified through mass spectrometry for either TB or LTBI cases. To down-size the number of metabolites analyzed, only metabolites that belong to metabolic groups composed of at least 7 metabolites and that are not labelled as ‘Drug’ or ‘Chemical’ were considered. The metabolite and protein correlation heatmap was plotted with the ggplot2 R package (v.3.3.1).