Baseline gene expression in subcutaneous adipose tissue predicts diet-induced weight loss in individuals with obesity

Study design

DiOGenes is a multicenter, randomized, dietary intervention study (trial number NCT00390637), described in detail elsewhere (Larsen et al., 2010). Briefly, 938 healthy adults with overweight or obesity (BMI 27–45 kg/m²) without concomitant medications entered the first phase of the study: a 2-month weight-loss period using a complete meal replacement low calorie diet (LCD) providing 800 kcal/day (Modifast; Nutrition et Santé France, Revel, France). Among the 781 participants who completed the LCD period, 773 achieved >8% weight loss and those individuals were randomly assigned to one of five different weight maintenance diets for the next 6 months, differing in macronutrient compositions (a 2 × 2 factorial combination diet of low/high protein and a low/high glycemic index or a control diet). A total of 548 (71%) participants completed the whole 8-month intervention. Here, we examined data from 281 out of these 548 participants, for whom abdominal subcutaneous adipose tissue RNA sequencing data at baseline was available. Weight development and metabolic health were assessed at baseline, 2, and 8 months.

The ethics committee of each center/country approved the protocol. The committees included (1) Medical ethical commission from Maastricht University, NL (2) Copenhagen ethical research commission, DK (3) Bedfordshire local Research Ethics Committee, Luton and Dunstable Hospital NHS Trust, UK (4) Ethics Committee of the Faculty Hospital, Prague University, CZ (5) Ethical Commission by NMTI, Sofia, BG (6) Ethical Commission University Potsdam, D (7) Ethical Commission Medical University, Navarra, SP (8) Scientific council Heraklion general university hospital, Heraklion, GR and (9) Commission Cantonale d’ éthique de la recherche sur l’ être humain, Canton de Vaud, CH. Furthermore, the protocol was in accordance with the Declaration of Helsinki. All study participants provided written consent.

Subcutaneous adipose tissue biopsy and transcriptomic analyses

Abdominal subcutaneous adipose tissue needle biopsies were collected 6–8-cm laterally from the umbilicus under local anesthesia. Thereafter, biopsy specimens were snap-frozen in liquid nitrogen and stored at −80 °C until further analysis. Total RNA was extracted and analyzed from subcutaneous adipose tissue biopsy specimens as described elsewhere (Armenise et al., 2017; van der Kolk et al., 2019). For each sample, gene expression was then examined using 100-nucleotide-long paired-end RNA sequencing with an Illumina HiSeq 2000 of libraries prepared using the Illumina TruSeq kit according to the manufacturer’s instructions. We used GenomicAlignments to retrieve the number of reads mapping onto 53,343 genes (GRCh37 assembly). Only reads with both ends mapping onto coordinates of a single gene were considered (Armenise et al., 2017; van der Kolk et al., 2019). The data is accessible through Gene Expression Omnibus under the accession no. GSE95640.

Construction, training and testing of weight loss prediction models

A diagram illustrating our method is shown in Figs. 1A and 1B. The RNAseq data from the studied individuals at baseline were split into training and testing subsamples. To prevent overfitting of the prediction models, we used a five-fold cross-validation with 100 resampling (i.e., 100 times our samples were randomly divided to 80% training and 20% testing subsets). This machine learning setup was used for all the prediction models that are reported in this article. Within each training and testing subset, participants were divided into high-weight losers (WL) and low-WL groups based on the cutoff of median weight loss percentage in the subset.

We used support vector machines (SVM) with linear kernel (and cost value of 1) to build classifier prediction models. Prior to the training/testing of the models, we applied the ComBatSeq() function supported by the SVA R library (Johnson, Li & Rabinovic, 2006; Zhang, Parmigiani & Johnson, 2020) on the count matrix to adjust for biases and effects introduced by sex, age and study center. We then applied Log₂ Count Per Million (LCPM) scaling while adjusting the library sizes by normalization factors that were calculated using the Trimmed Mean of M values (TMM) (Robinson, McCarthy & Smyth, 2010). The normalizations were carried out independently for the 80% training and the remaining 20% testing sub-samples. Furthermore, the training and testing of the prediction models were carried out on their corresponding sub-samples. We finally used the ROCR R-package to measure the area under the receiver operating characteristic (ROC) curve, i.e., AUC. The 95% confidence interval (CI) of the AUC of the most optimal models were calculated using the DeLong method supported by the pROC R package (DeLong, DeLong & Clarke-Pearson, 1988).

To ensure that the performances of the models are not affected by unaccounted for biases within the expression data, the performances of all prediction models (based on the studied gene-lists, see below) were compared to those that featured randomly selected genes as their input parameters. In further detail, we created 100 lists of randomly picked genes (each list featuring ten genes) over the whole genome, from the genes with >5 mapped reads across all samples. Next, we built 100 prediction models based on the randomly selected genes. For each prediction model a P-value was measured. The P-value is the fraction of the 100 prediction models that were based on the randomly picked genes that resulted in higher median AUCs compared to the prediction model (i.e., based on a studied gene-list). We considered P < 0.05 as significant.

Selection of genes for the prediction models

We studied differentially expressed genes (DEGs), comparing high-WL to low-WL, by applying DESeq2 on the training samples in each cross-validation run (Love, Huber & Anders, 2014). We adjusted the models for sex, age and study center and we considered the nominal P < 0.01 as significant. We then extracted the biological process Gene Ontology (GO) categories enriched in the significant DEGs, using the enrichGO() function supported by the clusterProfiler R/Bioconductor package (with the default parameter settings of P < 0.05 and q < 0.2 cutoffs) (Yu et al., 2012). This function applies hypergeometric distribution to calculate P-values based on the number of selected genes associated with the GO category, the overall number of genes associated with the GO category, and the overall number of genes in background (Boyle et al., 2004). The Benjamini Hochberg method was used to correct for multiple testing (Benjamini & Hochberg, 1995) and we considered FDR <0.05 significant. Furthermore, we constructed enrichment maps based on the similarities of the ontologies using the enrichplot R package and the redundant GO categories were filtered based on a ≥70% similarity cutoff. We applied the Jaccard similarity coefficient method supported by the pairwise_termsim() function of the enrichplot R package to construct the similarity matrix of the GO categories (Yu & Hu, 2020). The ‘components’ within the enrichment maps, i.e., biological process GO categories that for each pair at least one path exists, were finally extracted. We used the vertices within these components, i.e., classes of interrelated biological process GO categories, to extract the suitable genes for construction of the prediction models.

The biological process GO categories discovered from comparing the high-WL vs. low-WL include fatty acid metabolic pathways (referred to as Lipid), mitotic/meiotic pathways (referred to as Mitosis), and ‘response to virus’ pathways (referred to as Virus). Expression data of the genes that are associated with a pathway (i.e., Lipid, Mitosis or Virus) were extracted and principle component analysis (PCA) was performed. From each list of genes associated with a pathway, ten genes that most contributed to the PC dimensions that collectively explained ≥70% of the variance within the expression data was chosen. Since we used five-fold cross-validation with 100 resampling, the training data varied across the different runs and the discovered genes from a specific pathway across the different runs somewhat varied too. Therefore, from each studied Lipid/Mitosis/Virus gene list, the ten most frequently featured genes during the 100 runs were selected as input parameters to the corresponding prediction models. Finally, to further analyze the performance of the prediction models based on lipid metabolism and ‘response to virus’ genes, we also constructed models based on combined gene lists (named as Lipid + Virus), and prediction models based on genes featured in one but absent from the other (i.e., Lipid—Virus and Virus—Lipid).

Incorporating clinical factors into the prediction models

We studied the effects of anthropometric and clinical factors (Table S1 and S2) by building a prediction model solely based on these factors as well as incorporating them into the input parameters of the Lipid and Virus models. The values of the clinical factors were adjusted for sex and age using the non-parametric method supported by ComBat() function of the sva R/Bioconductor package (Johnson, Li & Rabinovic, 2006; Zhang, Parmigiani & Johnson, 2020). The highly skewed clinical data (i.e., skewness > 1 or skewness < −1) were log_e scaled. For each Lipid and Virus model that incorporated a clinical factor, a P-value was calculated by using paired t-test, by comparing their performances to the original Lipid/Virus prediction models that lack the clinical factor information. We considered P < 0.05 as significant.

Expression heatmaps

We normalized the gene expression within all samples across all three time points using variance stabilizing transformation (VST) supported by DESeq2, whilst adjusting for sex, age and study center in the design model (Anders & Huber, 2010).

Baseline characteristics and overall clinical outcome

Baseline characteristics of the participants are described in Table 1. For descriptive reasons, we created this table based the median weight loss percentage achieved at 8-months (9.9%) from 281 participants. Sixty-five percent of the participants were females, and the participants were on average 42-years-old. At baseline, the BMI of included participants was 34.6 ± 0.2 kg/m². The high-WL and low-WL groups were comparable at baseline in term of weight and insulin sensitivity (Table 1). Upon LCD at the 2-month time point, both groups had achieved significant weight loss (>8% initial body weight); yet the high-WL group had significantly higher weight loss (mean difference 2.3%, P < 0.001). By design, following 8 months, the high-WL group lost significantly more of their initial body weight 15.4 ± 4.6% compared to the low-WL group who lost 5.8 ± 3.4% (P < 0.001).

Prediction models based on biological processes associated with weight-loss

To create a predictive model for diet-induced weight loss success based on global subcutaneous adipose tissue gene expression and biological processes, we extracted the significantly DEGs in the training subsamples used in each cross-validation run by comparing the high-WL to the low-WL subsamples at baseline. Subsequently, we created GO enrichment maps that illustrate the relationships of the discovered significantly enriched biological process GOs based on their pairwise gene overlaps. Throughout the 100 runs, the most frequently discovered biological process GO categories were associated with fatty acid metabolic pathways (i.e., Lipid, discovered in 42 runs), mitosis/meiosis pathways (Mitosis, discovered in 36 runs), and ‘response to virus’ pathways (Virus, discovered in 27 runs) (Table 2). In 34 runs, however, no significant GOs were discovered. The largest nonoverlapping Lipid, Mitosis, and Virus associated biological process classes discovered across the 100 cross validation runs are illustrated in Fig. 2A.

Next, we extracted ten genes from each of the three discovered pathways (Fig. 2B and Table 3). We used the baseline expression of these genes as input parameters to construct weight-loss prediction models. The Lipid model (median AUC = 0.58, P < 0.01) (Fig. 3A) and Virus model (median AUC = 0.59, P < 0.01) (Fig. 3C) predicted the 8-months weight loss of the individuals in our study significantly better than 100 models that feature randomly selected genes as their input parameters or the Mitosis model (median AUC = 0.50, P = 0.21) (Fig. 3B). The maximum AUC achieved by the Lipid models was 0.74 with 95% CI [0.62–0.86], whilst for the Virus models was 0.72 with 95% CI [0.61–0.83] and for the Mitosis models was 0.66 with 95% CI [0.54–0.78].

Notably, the genes in the Virus model were highly in common with the genes in the Lipid models (Fig. 2B) indicating that the discovered ‘response to virus’ pathways have genes with heterogenous functions that are also associated with the lipid metabolism pathways. Therefore, we further examined prediction models that feature combinations of the Lipid and Virus genes. Prediction models based on combined ‘response to virus’ and lipid metabolism gene-lists (Lipid + Virus) did not perform significantly better than Virus or Lipid models (Figs. S1A, S1B). Furthermore, excluding the genes that are annotated with ‘response to virus’ from the Lipid models (Lipid–Virus) did not affect the performance of the Lipid metabolism significantly (Fig. S1C and Table S3). However, prediction models that featured genes associated with ‘response to virus’ but not annotated with lipid metabolism (Virus–Lipid) performed drastically different compared to the Virus models (Fig. S1D). The performance of these models was no longer significantly better than those models that feature randomly selected genes (median AUC = 0.51, max AUC = 0.67, 95% CI [0.51–0.75], P = 0.16) (Fig. S2A and Table S3). This shows that the performance of the Virus models is highly dependent on the genes that they have in common with the Lipid models. In contrast, the performance of the Lipid models are not dependent on the genes that they have in common with the Virus models (Fig. S2B and Table S3).

To unravel the conditions under which the Lipid and Virus prediction models best performed, and to explain their bimodal AUC distributions (Figs. 3A, 3C), we compared the training subsamples in the runs that resulted the AUCs higher than 0.60 vs. those in the runs that resulted lower AUCs. Note that the overall weight loss during the 8-months correlated highly with the weight loss during the 6-months weight maintenance phase (r = 0.90, P < 2.2e−16, Pearson correlation test) (Fig. 3D), but our results show that the models perform worst when the training dataset include more samples from individuals who gained weight during the 6-months weight maintenance phase despite being in the high-WL group. In detail, 30 high-WL individuals gained weight during the 6-months weight maintenance of which seven samples regained more than 4% (i.e., 5^th percentile) weight from the LCD-induced weight loss (Fig. 3D). We realized that runs with AUCs higher than 0.6 included significantly less samples from those seven individuals (P < 0.05, Jonckheere–Terpstra test) in their training sets compared to other runs. This was true for both, Lipid (P = 0.048, Jonckheere–Terpstra test) and Virus (P = 0.034, Jonckheere–Terpstra test) prediction models. These results indicate that training the models on data from high-WL individuals who regained weight considerably can lead to less efficient prediction models.

Incorporating clinical factors to improve the prediction models

Whilst most clinical parameters did not display strong separation between the High-WL and Low-WL groups (Table 1), we still hypothesized that the combination of easily accessible clinical parameters could predict weight loss. We constructed a prediction model based on 30 anthropometric and clinical factors, but this model did not predict the weight loss status of the samples significantly better than the models that feature randomly selected genes (median AUC = 0.51, max AUC = 0.75, 95% CI [0.60–0.90], P = 0.12) (Fig. S3).

We then hypothesized that adding these 30 anthropometric and clinical factors, one by one, to the input parameters of Lipid and Virus prediction models may yield to improved predictions. The following factors improved the performance of the Lipid models: systolic blood pressure (median AUC = 0.59, max AUC = 0.74, 95% CI [0.63–0.86], P = 0.037) and fibrinogen (median AUC = 0.59, max AUC = 0.77, 95% CI [0.67–0.88], P = 0.049) (Table S1). The performance of the Virus models was improved by adding systolic blood pressure (median AUC = 0.61, max AUC = 0.76, 95% CI [0.65–0.86], P = 0.018), diastolic blood pressure (median AUC = 0.59, max AUC = 0.74, 95% CI [0.65–0.86], P = 0.037) or waist-to-hip ratio (median AUC = 0.59, max AUC = 0.74, 95% CI [0.67–0.88], P = 0.037) (Table S2).

Notably, in most of the runs, no significant changes of performance were observed upon adding anthropometric and clinical factors (ΔAUC < 0.1) (Figs. 4A–4E). However, in one to six out of 100 runs where the Lipid and Virus models had performed poorly (AUC < 0.5), adding these parameters significantly improved the prediction model performances. In contrast, in a few runs, these parameters worsened the prediction performances where the Lipid and Virus models solely based on genes had performed well (AUC > 0.5). Overall, these results indicate that the effect of incorporating anthropometric and clinical factors into the weight-loss prediction model is negligible and should be carried out with caution.