Dulbecco’s modified Eagle medium (DMEM)/F-12 medium (1:1) enriched with L-glutamine and 15 mM 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES), fetal bovine serum (FBS) and penicillin streptomycin solution were purchased from Gibco by Life Technologies (Thermo Fisher Scientific Inc., Waltham, Massachusetts). TRIzol reagent, PureLink™ RNA Mini Kit and SuperScript™ III one-step RT-PCR system with Platinum™ Taq DNA polymerase were purchased from Invitrogen (Thermo Fisher Scientific Inc., Waltham, Massachusetts). Rosiglitazone was purchased from Cayman Chemical (Ann Arbor, Michigan). Allicin was purchased from Solarbio Life Sciences^® (Beijing, China).

All other chemicals used in the experiment and not listed above were purchased from Sigma-Aldrich (Darmstadt, Germany).

Human SGBS cells were grown Dulbecco’s Modified Eagle Medium Nutrient Mixture F-12 (DMEM/F-12) supplemented with 10% fetal bovine serum (FBS), 3.3 mM biotin, 1.7 mM panthotenate and 1% penicillin/streptomycin solution, at 37°C, 5% CO₂ and 95% relative humidity. Cells were platted in Petri dishes (100mm) in duplicate. Once the cells reached approximately 90% confluence, differentiation was induced by feeding the cells with serum-free growth medium supplemented with 10 µg/ml transferrin, 0.2 nM triiodothyronine (T₃), 250 nM hydroxycortisone, 20 nM human insulin, 25 nM dexamethasone, 250 µM 3-isobutyl-1-methylxanthine (IBMX) and 2 µM rosiglitazone (day 0 of differentiation). After 4 days, the differentiation medium was replaced with maintenance medium composed by serum-free growth medium supplemented with 10 µg/ml transferrin, 0.2 nM T₃, 250 nM hydroxycortisone and 20 nM human insulin. Fresh maintenance medium was added every 2 days.

Treatment with allicin began on day 0 of differentiation (D0), and continued until analysis on day six (D06) of differentiation ( Figure 1 ). Allicin concentrations of 5, 12.5, 25, and 50 µM were tested. Prior to treatment, allicin was diluted in dimethyl sulfoxide (DMSO) and a stock solution was prepared. Stock solutions were prepared so that the volume of DMSO in the treatment medium did not exceed 0.5%. Control cells (CTRL) were incubated with the same volume of 0.5% DMSO in the differentiation medium for 6 days. As a positive control (cAMP), SGBS cells were treated with 500 µM dibutyryl cAMP sodium salt (a cyclic nucleotide derivative that mimics endogenous cAMP) for 24 hours before the sixth day of differentiation (24) ( Figure 1 ).

Cell viability was determined using the 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay. Cells were plated in a 96-well plate and treated with 5.0, 12.5, 25.0, and 50 µM allicin (25). Before incubation with 5 mg/mL MTT in HBSS, cells were rinsed with 1X phosphate buffer saline (PBS) 1X. Incubation with the MTT solution was performed at 37°C for 4 hours. The resulting formazan was dissolved in dimethyl sulfoxide (DMSO) and incubated overnight (O/N) at 37°C. Optical density was used as an indicator of cell viability and was measured at 590 nm.

Cells for BODIPY™ staining and subsequent confocal imaging were cultured on ibiTreat 8-well μ-slides (Ibidi GmbH, Planegg/Martinsried, Germany). Cells were fixed in a 2% formalin solution diluted in PBS 1X at room temperature (RT) for 15 minutes. Subsequently, after washing three times in PBS 1X, the cells were incubated in a solution of BODIPY™ 493/503 in PBS 1X to fluorescently label the lipid droplets. Incubation was performed for 45 minutes in the dark at RT. The slides were then washed three times in PBS 1X.

Fluorescence images were acquired using a Leica SP8 confocal microscope (Leica Microsystems Srl, Milan, Italy) and LAS X 3.1.5.16308 software. Slides were viewed with the HCX PL APO lambda blue 63x/1.40 OIL objective. DAPI fluorescence was detected with a 405 diode laser (410/480 nm), while BODIPY fluorescence was detected with a white light laser (503/588 nm). Images were acquired using a photomultiplier tube (PMT) that allowed point scanning of the region of interest (ROI) with the selected laser and produced 1024 × 1024 px images.

The MRI_Lipid Droplets Tool (http://dev.mri.cnrs.fr/projects/imagej-macros/wiki/LipidDroplets_Tool), a macro in the ImageJ 1.50b software (http://rsb.info.nih.gov/ij/), was used to measure LD area (26). Individual cells were defined by regions of interest (ROIs) and images were analyzed as previously described (27). For each cell, the LD area (in μm²), maximum Feret diameter (MFD, in μm), and integrated optical density (IOD, dimensionless) were measured. The MDF is used as a measure of the diameter of irregularly shaped objects, whereas the IOD is related to both triglyceride accumulation and the size of LDs (28).

SGBS cells cultured on ibiTreat 8-well μ-slides (Ibidi GmbH, Germany) were incubated at 37°C with 100 nM MitoTracker^® Orange CMTMRos (Thermo Fisher Scientific, USA) for 30 minutes. The stained cells were washed with PBS 1X and fixed with 2% formalin at RT for 15 minutes. After fixation, cells were rinsed three times with PBS 1X, mounted in DAPI-containing mounting medium (Cayman Chemical Company, USA), and imaged using a Leica SP8 confocal microscope (Leica Microsystems, Germany) and LAS X 3.1.5.16308 software. Slides were viewed with the HCX PL APO lambda blue 63x/1.40 OIL objective. DAPI fluorescence was detected with a 405 diode laser (410/480 nm), while MitoTracker^® fluorescence was detected with a white light laser (550/605 nm). Images were acquired using a photomultiplier tube (PMT) that allowed point-by-point scanning of the region of interest (ROI) with the selected laser and produced images with a resolution of 1024 x 1024 px.

To quantify mitochondrial morphology on standard confocal fluorescence microscopy images of CTRL, ALLI, and cAMP-treated cells, the Mitochondrial Analyzer based on adaptive thresholding and the ImageJ/Fiji open-source image analysis platform were used (29). Scale was set for magnification and the global check box in the Set Scale dialog box was selected. After 2D threshold optimization, the images were thresholded with a block size of 1,350µm and a C-value of 5. Subsequently, the images were also processed using the MiNa (30) and Micro2P (31) tools.

The Mitochondria Analyzer tool was used to measure counts (number of mitochondria in the image), total area (sum of the area of all mitochondria in the image), mean area (total area/mitochondria number), total perimeter (sum of perimeter of all mitochondria in the image), mean perimeter (total perimeter/mitochondria number), mean aspect ratio (shape descriptor measuring elongation), and mean form factor (shape descriptor measuring round to filamentous shape). In addition, parameters describing the connectivity of the mitochondrial network were calculated, including the number of branches, the total length of branches, the mean length of branches, the branch junctions, the end points of branches, and the mean diameter of branches. The branches consist of point-shaped objects without branching junctions and minimal length, long single tubular objects without branching junctions, but the highest branch length and complex objects with multiple branches and junctions. The number of branches, total branch length, branch junctions and branch end points were also expressed as normalization to either the number of mitochondria or total area (29).

The MiNa tool, a macro of the ImageJ1.53o software (http://rsb.info.nih.gov/ij/), was also used to quantify mitochondrial morphology (30). Threshold images were processed using the Tophat option (32), as was the MiNa interface. The macro detected ‘individual’ mitochondrial structures in a skeletonized image, such as punctate, rod-shaped, and large/round structures without branching, and ‘networks’ identified as mitochondrial structures with a single node and three branches. All parameters were used in the discriminant analysis. Among the nine parameters calculated by MiNa, the number of individuals (punctate, rod-shaped, and large/round mitochondria), the number of networks (objects with at least one branch), and the mean rod/branch length, which refers to the average length of all mitochondrial rods/branches, were considered for statistical comparisons. Other parameters such as the mean number of branches per network, i.e., the mean number of mitochondrial branches per network, the mean length of branches, and the mitochondrial footprint, which refers to the total area of mitochondria, were included in the discriminant analysis.

Mitochondria were also analyzed and classified using MicroP software, a useful tool validated in CHO-K1 cells (31). The software classifies six morphological types of mitochondria, such as small globules, round-shaped mitochondria, that may have arisen by fission; large globules with a larger area; simple tubules, i.e. straight, elongated mitochondria without branches; twisted tubules, elongated tubular mitochondria with a non-linear development; donuts, like elongated tubules mitochondria but with fused ends; branched tubules, complex interconnected mitochondria with a network-like structure. On each image, the total number of mitochondria and their area were calculated as the ratio of mitochondria and area for each subtype in the different SGBS cells treated. These data with the number of objects and total area were used for discriminant analysis.

The experiment was set up with 2 biological replicates for the 3 experimental conditions.

After removing the culture medium from the Petri dishes, 1ml/10cm² of TRIzol reagent was added to each plate and repeatedly pipetted to induce a severe breakdown of the cell structures. These samples were immediately processed further using the PureLink™ RNA Mini Kit according to the manufacturer’s instructions.

The concentration of extracted total RNA was quantified using a spectrophotometer (NanoDrop 1000 Spectrophotometer, ThermoScientific, Wilmington, Delaware), and the purity of the RNA samples ranged from 1.8 to 1.9. RNA integrity was assessed by observing the 18S and 28S ribosomal bands after electrophoresis on 1% agarose gel, in the presence of GelRed. β-actin expression was used as an internal control, and confirmed the complete integrity of the RNA.

The purified total RNA was subjected to deep sequencing analysis. First, the isolated RNA was quantified using Agilent Bioanalyzer 2100 with the RNA integrity number (RIN) greater than 8.0 before sequencing using Illumina Genome Analyzer (GA). Generally, 2-4 ug of the total RNA was used for library construction. Total RNA was reverse transcribed into double-stranded cDNA, digested with NlaIII and ligated to an Illumina specific adapter containing a recognition site of MmeI. After MmeI digestion, a second Illumina adapter, containing a 2-bp degenerate 3’ overhang was ligated. The obtained sequences were aligned on GRCh38 human genome (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.39) using STAR software (33).

Raw data were uploaded to the R package (v0.92) Differential Expression and Pathway analysis (iDEP951) that is a web-based tool available at http://bioinformatics.sdstate.edu/idep/(34, 35). In the pre-processing step, genes expressed at very low levels across samples were filtered out, and genes expressed at a minimum of 0.5 counts per million (CPM) in a library were further analyzed. To reduce variability and normalizecount data, EdgeR log2(CPM+c) was chosen with pseudocount c = 4 transformation,. Next, the DESeq2 package in the R language was used to identify differentially expressed genes (DEG) between ALLI_ cAMP, ALLI_CTRL and cAMP_CTRL using a of false discovery rate (FDR) threshold ≤ 0.05 and fold-change > |1.0|. Heatmaps, principal component analysis (PCA), k-means cluster and enrichment analyses were also performed in iDEP951.

Gene set enrichment analysis to determine the shared biological functions of differentially regulated genes based on significant GO terms (36), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (37) and TF. target.TRED analyses were performed (38).

Venn diagrams were created by web tool available at http://bioinformatics.psb.ugent.be/webtools/Venn/.

On the basis of the online tool Search Tool for the Retrieval of Interacting Genes (STRING; https://string-db.org/), PPI networks of the up regulated and the down regulated DEGs in each comparison were generated with a confidence level of ≥ 0.4, and the PPI network was visualized using Cytoscape software (version 3.9.1, https://cytoscape.org/). Then, the PPI networks of DEGs in each comparison were analyzed using the Cytoscape CytoHubba plugin to select the top 10 hub nodes according to the Degree algorithm (39). The Molecular Complex Detection (MCODE) plug-in (40) in the Cytoscape suite was used to examine the significant modules in the PPI network of overlapping DEGs that are the target ofshared TF networks between comparisons. Degree cutoff = 2, K-core = 2, and node score cutoff = 0.2 were set as options. Enrichment analysis of DEGs in modules with a score ≥ 5 was then performed.

Estimation of the proportion of brown adipocytes in each sample was analyzed based on read counts using the PROFAT tool, which automatically performs hierarchical cluster analysis to predict the browning capacity of mouse and human RNA-seq datasets (41).

To identify a larger number of potential targets, PharmMapper (http://www.lilab-ecust.cn/pharmmapper/; 42, 43), the similarity ensemble approach (SEA, https://sea.bkslab.org/), the STITCH database (http://stitch.embl.de/; 44), Swiss Target Prediction (http://www.swisstargetprediction.ch/; 45), and GeneCard (https://www.genecards.org/) were used. The 2D structure and canonical SMILES of allicin (CID_65036), diallyl sulphide (CID_11617), diallyl disulfide (CID_16590), and diallyl trisulfide (CID_16315) were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). The sdf files were uploaded to the PharmMapper server, and the search was started using the maximum generated conformations of 300 by selecting the option ‘Human Protein Targets Only (v2010, 2241)’ and the default value of 300 for the number of reserved matching targets. for the other parameters, the ‘default mode’ was selected. Canonical SMILES were uploaded to the other tools. The predicted targets were entered into the UniProt database (https://www.uniprot.org/) with the species set to Homo sapiens to determine their gene IDs. A Venn diagram was used to find common targets among the allicin compounds. Genes related to ‘adipocyte’, ‘browning’, ‘non shivering thermogenesis’, ‘cold-induced thermogenesis’, ‘brown adipose thermogenesis,’ and ‘adaptive thermogenesis’ were downloaded from GeneCard, and the intersection of the targets was determined using the Venn tool. The resulting common genes associated with browning and adipocytes were then crossed with common putative targets of allicin, and a Venn diagram was generated. Subsequently, the overlapping targets were uploaded to GeneMANIA (https://genemania.org/) (46) to perform functional gene analysis and generate a PPI network.

CytoNCA, another Cytoscape plugin, was applied to the network to perform topological analysis evaluating the centrality measures of the network (47). Then, the Cytoscape intersectional merge function was used to isolate the PPI subnetworks. Key node functions were determined by analyzing GO terms, KEGG and Reactome pathways.

By entering the screened key nodes into the online tool VarElect (48), the correlation between nodes and ‘cold induced thermogenesis’ was investigated.

All measurement results are given as means ± SD and were analyzed with XLSTAT (49). Measurements of LD area surface/cell, MFD/cell, and IOD/cell obtained from 15 biological replicates were compared along with Mitochondrial Analyzer, MiNa, and Micro2P results using the Kruskal-Wallis statistical test, followed by pairwise comparisons using the Mann-Whitney approach with Bonferroni correction (p < 0.0167).

All Mitochondrial Analyzer and MiNa parameters, as well as ratios of parameters obtained with the Micro2P tool, were calculated together to perform a canonical discriminant analysis (DA) that integrates morphological mitochondrial parameters into a single multivariate model with the aim of maximizing differences between treatments and calculating the best discriminant components between treatments (49).

To investigate possible adverse effects of allicin, a viability assay was used to investigate possible adverse effects of ALLI extract on SGBS cells treated with doses of 5, 12.5, 25, and 50 µM. In particular, the analysis showed that the viability of cells treated with ALLI extract decreased significantly (p < 0.001) in a dose-dependent manner ( Figure S1 ), with the 50 µM dose always significantly different from the other doses. Nevertheless, viability remained above 85% up to 25 µM and 12.5 µM ALLI and did not differ from 5 and 25 µM. Based on these findings, 12.5 µM ALLI was selected for further experiments.

Next, We investigated whether allicin has an effect on early adipogenic differentiation. Thus, we performed LD analysis in SGBS cells after ALLI treatment during the induction of adipogenesis. Figure 1 shows statistically significant differences in area of LDs/cell, MFD/cell, and IOD/cell between treatments and illustrates Bodipy staining in SGBS cells after 6 days of treatment with allicin (ALLI), CTRL, and dibutyryl cAMP (cAMP). The area of LDs/cell was significantly lower in cells treated with cAMP and ALLI compared with cells from CTRL (p < 0.0001). No significant differences were observed between cAMP- and ALLI-treated cells ( Figure 2A ). A significant (p < 0.0001) decrease in MDF/cell was observed in cells treated with cAMP and ALLI compared with CTRL (p < 0.0001) ( Figure 2B ), while IOD/cell was significantly lower in cells treated with cAMP compared with cells treated with ALLI (p = 0.016) and CTRL (p = 0.0001) ( Figure 2C ). A significant increase in the number of LDs/cell (p = 0.015) was observed between ALLI treated cells and CTRL cells ( Figure 2D ).

Accurate analysis of mitochondria is critical for determining mitochondrial dynamics, so three different tools were used to characterize mitochondrial morphology. Figure S2 shows the identification and classification of objects using adaptive thresholding (radius = 1.350 µm, C = 5) on images of ALLI- and cAMP-treated cells and CTRL. Figures S2A, B show the effects of the different treatments measured with the three tools, on the total number and area of mitochondria in SGBS cells. No statistical differences in mitochondrial distribution patterns were detected between treatments, although treatments with allicin and cAMP caused a shift toward a greater number of mitochondria ( Figure S2A ) and a greater total area of mitochondria compared with cells from CTRL ( Figure S2B ). These results suggest that ALLI treatment increases the number of mitochondria by inducing mitochondrial fission or biogenesis and thus also increases the total area of mitochondria.

Therefore, more than 26000 mitochondria were analysed using Micro2P software to classify six morphological subtypes and calculate the average proportion of mitochondrial subtypes within each treatment ( Figure 3A, B ). ALLI and cAMP treatments resulted in 13.6% (p < 0.01) and 11.5% (p < 0.05) more small globules, respectively, compared with CTRL. ALLI treatment reduced the percentage of mitochondria with large globules to 32.5% and that of the simple tube subtype to 19.5% compared with CTRL (p < 0.01) ( Figure 3A ). Accordingly, the ratio of surface area was significantly higher (p < 0.05) in ALLI- and cAMP-treated cells compared with CTRL at 30.3 and 31.7%, respectively, whereas the area of large globules in ALLI-treated cells decreased significantly (p < 0.05) ( Figure 3B ).

Specifically, the average ratio of mitochondrial amount in the different treatments was 53.4% in ALLI-treated cells, 28.6% in cAMP-treated cells, and 18.02% in CTRL. The relative percent area of total mitochondrial content between treatments was 51.9% in ALLI-treated cells, 27.8% in cAMP-treated cells, and 20.3% in cells from CTRL. Considering all treatments together, small globules (65.7%) and simple tubules (23.3%) were the most representative subtypes, followed by branched tubules (5.42%), large globules (2.3%), donuts (1.7%), and twisting tubules (1.6%). Accordingly, the percentage area of mitochondrial subtypes was 34.9% for small globules, 33.9.8% for simple tubules, 16.67% in branching tubules, 7.4% for donuts, 3.6% for large globules, and 3.5% for twisting tubules.

Mitochondrial Analyzer detected significantly higher numbers of mitochondria in ALLI-treated cells ( Figure 3C ), whereas the mean perimeter and area ( Figure 3D ) were significantly lower compared with cells from CTRL. No statistical differences were observed in the cells treated with cAMP. The shape of mitochondria, characterized by aspect ratio and form factor, decreased significantly under ALLI and cAMP compared with cells from CTRL ( Figures 3E, F ). Mitochondria Analyser tool was also used to quantify the morphological complexity of mitochondria. The mean length of branches and the total number of branches per mitochondrion were lowest in cells treated with ALLI and showed a significant difference in cells treated with cAMP (data not shown).

Network parameters calculated with MiNa confirmed a significant (p < 0.05) increase in the number of individuals (puncta, rods and large) in ALLI-treated cells compared with CTRL cells ( Figure 4A ). The number of networks showed no significant differences between treatments ( Figure 4B ), but the mean value of rod/branch length significantly decreased in cells treated with ALLI (p < 0.05) and cAMP (p < 0.05) ( Figure 4C ).

To determine the extent of treatment-specific variation in an unbiased manner, all Mitochondrial Analyzer measurements, six MiNa measurements from, and the percentage of the number and area of each mitochondrial subtype, number of objects, and total area calculated by Micro2P were integrated into a single multivariate model (DA) that maximized differences between treatments ( Figure 5 ).

The model achieved clear separation of samples. With few exceptions, the treatments were grouped and divided into 3 clusters. Mitochondrial values of ALLI-treated cells were projected in the quadrant with a positive value for the F1 and F2 components, CTRL was projected in the quadrant with a negative value for F1 and a positive value for the F2 component. The projection of cAMP-treated cells was observed in the quadrant with negative value for F1 and positive value for F2. The significant morphological characteristics that distinguished the different experimental groups are shown in Table 1 . This analysis showed that most of the variance between treatments in mitochondrial dynamics was based on elongation index, area of small globes, solidity of simple tubes, number of punctate and rod-shaped mitochondria (individuals), and number of simple tubes ( Table 1 ). Discriminant functions were used to classify the treatments into the correct groups. To check the functionality and robustness of the classification model, a cross-validation test was performed, in which the success rate of correctly classified samples was 88.6%.

After removal of low abundance reads, the final mapping rate of the filtered transcript reads was 71.2%. Hierarchical clustering was performed on the initial analysis of the RNAseq results, that showed that the transcriptomic data was well-clustered according to the treatment. In addition, PCA showed the overall variability in the expression profile of the samples and treatments. Overall, there was a significant difference between CTRL and the treatments with cAMP and ALLI along the first principal component, which explained 61% of the variance, with smaller differences 28% along the second component ( Figure S3 ).

Analysis of DEGs was performed between cells treated with allicin and cAMP compared to control cells and between allicin and cAMP treatments, based on a log₂fold change of |1| and an FDR adjusted p value of ≤ 0.05.

iDEP951 expression analysis significantly identified (p < 0.05) 820 up regulated genes between cells treated with ALLI vs cAMP, 1417 up regulated genes between cells treated with ALLI vs CTRL and 1647 between cells treated with cAMP vs CTRL cells. Significantly down regulated (p < 0.05) genes were 640 between cells treated with ALLI compared with cAMP, 1085 genes between cells treated with ALLI vs CTRL and 1521 between cells treated with cAMP compared with CTRL cells ( Supplementary Table 1 ).

To further investigate the function of the DEGs, GO term enrichment analysis was performed. The up- and down regulated DEGs were significantly enriched in biological processes (BP), molecular functions (MF) and cellular components (CC).

In all comparisons, 30 BP, CC and MF GO terms were found to be significantly enriched ( Supplementary Tables 2 – 4 ). Notably, BP terms such as ‘Oxoacid metabolic process’, ‘Small molecule metabolic process’, ‘Fatty acid metabolic process’, ‘Cellular respiration’, ‘Cellular lipid metabolic process’, ‘Energy derivation by oxidation of organic compounds’, ‘Monocarboxylic acid metabolic process’ were down regulated in the ALLI_cAMP comparison, but up regulated in the ALLI_CTRL and cAMP_ CTRL comparisons ( Supplementary Table 2 ). Interestingly ‘Mitochondrion’, ‘Mitochondrial matrix’, ‘Mitochondrial inner membrane’, ‘Mitochondrial membrane’, ‘Mitochondrial envelope’, ‘Mitochondrial protein-containing complex’ CC terms were down regulated in ALLI_cAMP comparison, but up regulated in ALLI_CTRL and in cAMP_ CTRL comparisons ( Supplementary Table 3 ).

In contrast, significantly enriched MF terms such as ‘Calcium ion binding’, ‘Cell adhesion molecule binding’, ‘Collagen binding’, ‘Extracellular matrix structural constituent’, ‘Glycosaminoglycan binding’, ‘Growth factor binding’, ‘Heparin binding’, ‘Integrin binding’ and ‘Signalling receptor binding’ were up regulated in the ALLI_cAMP comparison, but down regulated in the ALLI_CTRL and in the cAMP_CTRL comparisons ( Supplementary Table 4 ). Interestingly, significantly enriched MF terms such as ‘Active transmembrane transporter activity’, ‘Electron transfer activity’, ‘Oxidoreductase activity’ and ‘Transmembrane transporter activity’ were up regulated in the ALLI_CTRL and in the cAMP_CTRL comparisons ( Supplementary Table 4 ).

Figure 6 shows the results of the KEGG analysis in the form of a graphical representation of the scatter plots. Each figure shows the KEGG enrichment of 15 identified pathways for each treatment comparison with the corresponding GeneRatio, adjusted p-value, and number of enriched genes in the corresponding pathways. The GeneRatio is defined as the number of enriched candidate genes compared with the total number of annotated genes included in the corresponding pathway in the KEGG analysis. Therefore, a higher GeneRatio indicates a higher enrichment of candidate genes in the corresponding pathway. KEGG analysis showed that DEGs were significantly down regulated within the ‘PPAR pathway,’ Fatty acid metabolism, elongation and degradation, and in ‘Citrate cycle (TCA cycle)’ in ALLI_cAMP comparison. Of note, DEGs were over-expressed in pathways such as ‘cAMP signalling pathway,’ ‘Calcium signalling pathway,’ and ‘CGMP-PKG signalling pathway’ in ALLI_cAMP comparison ( Figure 6 ). ALLI_CTRL and cAMP_CTRL had common DEGs significantly enriched in ‘Oxidative phosphorylation’ and ‘Thermogenesis’, whereas DEGs within the ‘PPAR signalling pathway’ were down regulated in ALLI vs cAMP and up regulated only in cells under cAMP treatment ( Figure 6 ).

These results suggest that the browning effect of ALLI is only evident when compared with CTRL cells, so the ALLI-cAMP contrast was not discussed further.

The PPIs of all up regulated and down regulated DEGs with combined scores greater than 0.4 were constructed from the three comparisons, and each entire PPI network was analysed using Cytohubba. The ten most highly regulated hub genes with a high degree of connectivity between nodes are listed in Table 2 . The highly regulated hub genes in the ALLI_CTRL and cAMP_CTRL comparisons shared 6 genes such as PPARG, FASN, SREBF1, SCD, PPARGC1A, and ACLY. Both comparisons, referring to CTRL, were similarly enriched in ‘Fatty acid synthase complex,’ ‘acetyl-CoA carboxylase complex,’ ‘AMPK signalling pathway,’ ‘PPARA activates gene expression,’ ‘Regulation of small molecule metabolic process,’ and ‘Thermogenesis’.

Four down regulated hub genes, such as FN1, THBS1, COL1A2, and CCN2, are common between ALLI- and cAMP-treated cells compared with CTRL cells ( Table 3 ). Enrichment of both comparisons included ‘AGE-RAGE signalling pathway in diabetic complications’, ‘PI3K-Akt signalling pathway’, ‘Focal adhesion’, ‘ECM-receptor interaction’, and ‘TGF-beta signalling pathway’.

TRED analysis (http://rulai.cshl.edu/TRED) allows to know interaction data between transcription factors (TFs) and the promoters of their target genes, including binding motifs (38). In the ALLI_cAMP comparison, the target genes of only 1 TF (c-MYC) were down regulated and 15 were up regulated; in the ALLI_CTRL comparison, 4 TFs were down regulated and 3 (AR, PPARA, PPAG) were up regulated; in the cAMP_CTRL, 15 were down regulated and the same 3 of the ALLI_CTRL comparison were up regulated ( Supplementary Table 5 ).

The up regulated genes were enriched in the target of AR in all comparisons (14 genes), whereas genes enriched in the target of PPARG (37 genes) and PPARA (18 genes) were enriched only in the ALLI_CTRL and cAMP_ CTRL comparisons ( Supplementary Table 5 ). The target genes of JUN, SP1, and TP53 were significantly down regulated in the ALLI_CTRL and cAMP_CTRL comparisons but up regulated in the ALLI_cAMP comparison ( Supplementary Table 5 ). Down regulated DEGs were enriched in target genes of EGR1, ETS1, SMAD1, SMAD3, SMAD4, and TFAP2A in the cAMP_CTRL comparison, but up regulated in ALLI_cAMP ( Supplementary Table 5 ).

Enrichment analysis of the 14 genes targeting to AR and common to all comparisons showed significant up-regulation of ‘Response to hormone, ‘Regulation of lipid metabolic process’, ‘Regulation of small molecule metabolic process’, ‘Cellular response to hormone stimulus’, ‘Zinc finger nuclear hormone receptor-type’ and ‘PPARA activates gene expression’. The MCODE plugin cluster analysis did not filter any cluster with a score ≥ 5 for these genes.

The 37 DEGs targeting AR that were common only between ALLI and cAMP treatments compared with CTRL cells and the 51 common DEGs targeting PPARG resulted in only one MCODE cluster. No cluster with a score ≥ 5 was found for common DEGs targeting PPARA.

Enrichment analysis revealed the potential function of genes in each module. Shared DEGs targeted by AR and PPARG and over-expressed in ALLI_CTRL and cAMP_CTRL were enriched, among others, in ‘PPARG signalling pathway’ ( Figures 7A, B red), ‘Positive regulation of cold-induced thermogenesis’ ( Figures 7A, B , brown), ‘Fatty acid metabolic process’, ‘AMPK signalling pathway’ ( Figures 7A, B , green), ‘AMPK signalling pathway’ ( Figure 7A , blue) and ‘PPARA activates gene expression’ ( Figure 7B , blue).

DEGs targeted by TP53, JUN and SP1 were down regulated in both ALLI_CTRL and cAMP_CTRL comparisons, but up regulated in ALLI_cAMP comparison. Shared DEGs targeted by SP1 and down-expressed in ALLI_CTRL and cAMP_CTRL were enriched in ‘Interleukin-4 and Interleukin-13 signalling’ ( Figure 7C , red), and in ‘Extracellular matrix organization’ ( Figure 7C , red). Interestingly, DEGs targeted by SP1 were found over-expressed in ALLI_cAMP comparison. Down regulated DEGs targeted by TP53 and JUN common to ALLI_CTRL and cAMP_CTRL were enriched in ‘Interleukin-4 and Interleukin-13 signalling’ and ‘IL-18 signalling pathway’ (data not shown).

Allicin stimulation favors the differentiation into beige adipocyte The PROFAT tool (41) generated the heatmap of marker expression starting from normalized reads counts of SGBS cells. The estimation of BAT phenotype in ALLI- and cAMP-treated cells increased significantly (p < 0.0001) in comparison to CTRL cells. In contrast, WAT phenotype decreased significantly (p < 0.0001) in ALLI- and cAMP-treated cells compared with CTRL ( Figure 8 ). These results evidenced that SGBS cells exhibit a gene expression pattern similar to that of brown cells during 6 days of differentiation under allicin treatment.

Because allicin is rapidly converted in vitro to its related fat-soluble organosulfur compounds such as DATS, DADS and DAS, the potential targets of these compounds and of allicin were screened by computational target fishing from the PharmMapper, STITCH, Swiss Target Prediction and GeneCard databases. By overlapping the highest ranked common targets of allicin and related organosulfur compounds with the 315 overlapped ‘adipocyte-browning’ genes, 26 common targets between allicin compounds and adipocyte-browning were used to create a GeneMania network ( Supplementary Table 5 ). The results of the analysis showed that these 26 targets correlated with 20 others and a total of 407 different links were predicted to construct a network linking these 46 genes ( Figure 9A ). The constructed network had 33.99% physical interactions and 23.56% predicted functional relationships between genes. In addiction 20.58% shared the same protein domain and 13.85% shared similar co-expression characteristics, other results were pathways (5.24%) and colocalization (2.77%) as shown in Figure 9A . The molecular functions of the top ranked targets, filtered by their FDR score, were reported as GO categories. The preliminary network illustrated that the genes, depicted by different colours in Figure 9A , were involved in ‘ligand-activated transcription factor activity’, ‘intracellular receptor signalling pathway’, ‘temperature homeostasis’, ‘regulation of cold induced thermogenesis’, ‘reactive oxygen species metabolic process’, ‘positive regulation of lipid metabolic process’, and ‘cold-induced thermogenesis’.

A topological analysis of the preliminary network was performed using the CytoNCA plugin in Cytoscape to find the core proteins that form the preliminary network. The mean the degree (17.70), betweenness (48.96), closeness (0.49), eigenvector (0.106), LAC (12.40), network (10.95) values of the preliminary network was calculated, and the nodes of the preliminary PPI network that were above this mean were sorted out to build the corresponding subnetworks. Using the intersectional merge function in Cytoscape a core PPI subnetwork was extracted ( Figure 9B ) containing 6 key nodes (ESR1, NR1H3, NR1H4, PPARA, RXRA, RXRG) and 15 edges. Among these genes, Nuclear Receptor Subfamily 1 Group H Member 4 (NR1H4) and PPARA were significantly up regulated in the ALLI_CTRL comparison of SGBS cells ( Supplementary Table 1 ), whereas estrogen receptor 1 (ESR1), Nuclear Receptor Subfamily 1 Group H Member 3 (NR1H3) and Retinoid X receptor alpha (RXRA) are common targets of allicin and related fat-soluble organosulfur compounds.

GO terms from biological process, cellular component, and molecular functions were examinedand the most enriched GO terms from biological process were ‘intracellular receptor signalling pathway’, ‘cellular response to lipid’, ‘hormone-mediated signalling pathway’, ‘response to steroid hormone’, ‘response to lipid’, whereas the most enriched GO terms from cellular component and molecular functions were ‘transcription regulator complex’, ‘nuclear receptor activity’, and ligand-activated transcription factor activity, respectively (data not shown).

KEGG pathways were analysed with a redundancy cut-off of 0.7, 17 pathways were statistically significant (FDR < 0.05) ‘PPAR signalling pathway’, ‘Adipocytokine signalling pathway’, ‘Thyroid hormone signalling pathway’, ‘Non-alcoholic fatty liver disease, ‘Insulin resistance’ and ‘Lipid and atherosclerosis’(data not shown).

The pathways enriched by Reactome analysis were ‘Nuclear Receptor transcription pathway’, ‘PPARA activates gene expression’, ‘SUMOylation of intracellular receptors’, Regulation of lipid metabolism by PPARalpha’ and ‘Mitochondrial biogenesis’, which were shown as bubble plot combined with a Sankey diagram ( Figure 10 ).

In addition, the score of the 6 key nodes identified by the topological analysis was scored using VarElect. The score indicates the strength of the association between the target and the ‘cold induced thermogenesis’ phenotype ( Table 4 ). ESR1, PPARA, NR1H3, and NR1H4 scored > 6.