Genetic Polymorphisms, Mediterranean Diet and Microbiota-Associated Urolithin Metabotypes can Predict Obesity in Childhood-Adolescence.

Environmental and genetic factors are associated with pandemic obesity since childhood. However, the association of overweight-obesity with these factors, acting as a consortium, has been scarcely studied in children. We aimed here to assess the probabilities of being overweighed-obese in a randomly recruited cohort of Spanish children and adolescents (n = 415, 5-17 years-old) by estimating the odds ratios for different predictor variables, and their relative importance in the prediction. The predictor variables were ethnicity, age, sex, adherence to the Mediterranean diet (KIDMED), physical activity, urolithin metabotypes (UM-A, UM-B and UM-0) as biomarkers of the gut microbiota, and 53 single-nucleotide polymorphisms (SNPs) from 43 genes mainly related to obesity and cardiometabolic diseases. A proportional-odds logistic ordinal regression, validated through bootstrap, was used to model the data. While every variable was not independently associated with overweight-obesity, however, the ordinal logistic model revealed that overweight-obesity prevalence was related to being a young boy with either UM-B or UM-0, low KIDMED score and high contribution of a consortium of 24 SNPs, being rs1801253-ADRB1, rs4343-ACE, rs8061518-FTO, rs1130864-CRP, rs659366-UCP2, rs6131-SELP, rs12535708-LEP, rs1501299-ADIPOQ, rs708272-CETP and rs2241766-ADIPOQ the top-ten contributing SNPs. Additional research should confirm and complete this model by including dietary interventions and the individuals' gut microbiota composition.

Introduction

Obesity is an aetiological condition associated with some types of cancer and cardiometabolic diseases such as type-2 diabetes, metabolic syndrome, non-alcoholic steatohepatitis, and hypertension[1]. Nowadays, it is widely accepted that the combination of high-energy diets, genetic make-up, sedentary lifestyles and gut dysbiosis (impaired composition and functionality of the gut microbiota) are involved in the obesity pandemic[2-4]. The prevalence of these significant health threats has risen to shocking proportions worldwide, including countries like Spain with ancestral adherence to the Mediterranean diet[5]. The increase in the rate of obesity has been mainly attributed to the ‘Westernization’ of the diet, the decrease of physical activity from childhood, currently aggravated by the abuse of playing videogames, and the increasingly early access to digital devices such as smartphones[6]. However, other potential variables could be participating in this pandemic.

Many studies have associated obesity in children and adults with single nucleotide polymorphisms (SNPs). For example, in genome-wide association studies (GWAS), the rs9939609 SNP in the fat-mass-and-obesity-associated (FTO) gene has been reported to account for a modest, but a statistically significant, increase of 0.4 kg/m2 body mass index (BMI) units for each risk allele (A)[7]. Although it is not fully understood yet, this association is partially mediated via controlling feeding behaviour[8]. However, not all the studies report this association, including those conducted in the child population[9,10]. Indeed, more than one million SNPs have been detected in the human genome[11], and thus, instead of only one or few specific SNPs, it is more conceivable to expect the complex action of a consortium of SNPs potentially interacting with many other variables and associated with different conditions, including obesity.

The gut microbiota is also involved in the pathophysiology of obesity, although the associated mechanisms are not fully known yet[3]. Nevertheless, a number of pathways have been identified such as the translocation of lipopolysaccharides (LPS) from the gut to the bloodstream[12], the regulation of gut hormones, energy harvest, inflammatory responses, lipogenesis and immune interactions[13] as well as the regulation of white adipose tissue inflammation via microRNAs[14]. Recently, a potential nexus between the dissimilar metabolism of some dietary constituents by the microbiota and obesity has been proposed. This has been suggested for the metabolism of the polyphenols isoflavones[15] and ellagitannins[16-18] that yield specific metabolites[19], that is, equol and(or) O-desmethylangolensin (ODMA) in the case of isoflavones[20], and different urolithin combinations in the case of ellagitannins[17]. These particular types of metabolisms, in the context of polyphenols, give rise to specific metabolizing phenotypes: the so-called ‘metabotypes’[16,19]. A metabotype is characterized by the production of specific metabolites derived from the gut microbiota (which are characteristic of the precursor polyphenol metabolism), and also by the associated microbial ecology in terms of composition and activity[16,19]. Therefore, the capacity to excrete high or low amounts of microbial-derived metabolites is not within the concept of metabotype[19]. The definition of metabotypes in the metabolism of polyphenols, only fits currently for the equol and(or) ODMA ‘producers’ vs ‘non-producers’ metabotypes in the case of isoflavones[20], and also for the urolithin metabotypes associated with the metabolism of ellagitannins, i.e., metabotype A (UM-A; individuals that produce only urolithin A), B (UM-B; production of isourolithin A, urolithin B and also urolithin A) and 0 (UM-0; urolithin non-producers)[17,19] (Supplementary Fig. 1). The occurrence of specific gut microbiota metabotypes is behind the inter-individual variability upon polyphenol consumption[21,22] and could be indirect markers of gut dysbiosis reflecting the individuals’ gut microbiota composition, richness, diversity, and functionality[19,23-25]. Although the gut microbiota associated with UM-B and UM-0 individuals show a dysbiotic-prone pattern[25], however, the unequivocal association between these metabotypes with obesity has not been confirmed so far due to its multifactorial aetiology[25,26].

In the present study, we aimed to assess the probabilities of being overweighed or obese in a cohort of children and adolescents from the Southeast of Spain by estimating the odds ratios (ORs) for different predictor variables and their relative importance in the prediction of the response. In this proof-of-concept, we considered as predictor variables the urolithin metabotypes as biomarkers of the gut microbiota, ethnicity, age, sex, the adherence to the Mediterranean diet, physical activity, and a consortium of 53 SNPs from 43 genes mainly related to obesity and cardiometabolic diseases.

Methods

Study population

This research (‘The PolyMicroBio study’) was included in the Spanish National Project AGL2015-64124-R and complied with the ethical guidelines outlined in the Declaration of Helsinki and ethical principles for medical research involving human subjects (Seoul, Korea, 2008). The study was conceived to stratify the participants according to their urolithin metabotypes after three days of walnuts or pomegranate juice consumption[26] and was not intended to modify any variable in the children. The trial was registered at clinicaltrials.gov (NCT03318042), and the Spanish National Research Council’s Bioethics Committee (Madrid, Spain) approved the protocol. Inclusion criteria were ages from 5 to 17 years old and good health status. Exclusion criteria were diagnosed pathology, previous gastrointestinal surgery, chronic medication and antibiotic intake one month before participating. A total of 415 children and adolescents were randomly recruited. Children within the 5 to 12 years old group (n = 202) were recruited from the public primary school ‘CEIP Jara Carrillo’ (Alcantarilla, Murcia, Spain) and adolescents aged from 13 to 17 (n = 213) from the public high school ‘IES Alcántara’ (Alcantarilla, Murcia, Spain). Parents were fully informed and gave their written informed consent before the participation of all students.

Urolithin metabotypes

Children and adolescents consumed 25 g peeled raw walnuts daily or 250 mL of pomegranate juice daily (in the case of individuals allergic to nuts) for three days. Packs of peeled walnuts were kindly provided by Borges International Group, S.L. (Reus, Tarragona, Spain) and pomegranate juice by the AMC Group (Espinardo, Murcia, Spain). In the morning of the fourth day, a sample of urine was provided for its analysis by high-performance liquid chromatography with diode array detection coupled to electrospray ionisation and ion-trap tandem mass spectrometry (HPLC-DAD-ESI-IT-MS/MS), and ultra-high performance liquid chromatography coupled with electrospray ionization-quadrupole-time-of-flight-mass spectrometry (UPLC-ESI-QTOF-MS) as described elsewhere[27]. This allowed the stratification of the participants according to their different capacity to metabolise ellagic acid derivatives into urolithins, i.e., urolithin metabotypes UM-A, UM-B or UM-0 as previously described[26].

Anthropometric measurements and validated questionnaires

The determinations of height, weight, and waist and hip circumference were performed always by the same research staff, using the same equipment in all cases, and in the presence of teachers from the educational centres. The child growth standards from the World Health Organization (WHO) were used to define the BMI (kg/m2) cut-offs for underweight, normoweight, overweight and obese individuals as a function of sex and age[5]. The students were asked to record possible incidences (medication, protocol compliance, etc.), and also their physical activity level[28], which took into account the two hours of physical activity a week in their schools (low activity) and the practice of additional extracurricular sports at least three days a week (high activity). Besides, a validated questionnaire to assess the adherence to the Mediterranean diet in children (KIDMED) was used[29]. The score in this questionnaire (ranging from 1 to 13, from very poor to optimum adhesion, respectively) was grouped as ‘Low’ (score from 1 to 4), ‘Medium’ (from 5 to 8), and ‘Good’ (from 9 to 13).

Selection of SNPs and genotyping

Candidate genes and polymorphisms were identified after browsing the Single Nucleotide Polymorphism Database (dbSNP) and examining the published literature regarding each known gene and variant (favourable and unfavourable) associations[30-32], and analysing their potential regulatory and biological functions with the Regulome DB and HaploReg v4.1 databases[33,34]. On the same day of the anthropometric evaluation, saliva samples were obtained by gently rubbing the inside part of the cheek with a sterile swab, free of human RNA and DNA (Deltalab, Barcelona, Spain). Children were asked to clean their mouths and avoid eating or drinking 60 min before collection of samples to prevent contaminations. Two samples were obtained per student. The swabs were immediately stored in refrigeration and further frozen at −80 °C until their processing. Genomic DNA extraction and genotyping were carried out at the GENYAL Platform (IMDEA-Food, Madrid, Spain) using the OpenArray™ AccuFill™ System (Life Technologies Inc. Carlsbad, CA, USA) as described elsewhere[35]. Data analysis was made by TaqMan Genotyper Software v1.3 (autocaller confidence level >90%).

Statistical analysis

A proportional-odds logistic ordinal regression was used to model the data with the software R version 3.5.1 (www.r-project.org). Nine subjects with >40% missing SNPs were removed, resulting in a final sample size of 406 students (Supplementary Fig. 2). The missing data was singly imputed using the missForest R package. Redundant predictors (rs9928094-FTO, rs9935401-FTO) were identified and removed using the Hmisc R package. SNPs with either favourable or unfavourable genotype frequencies below 5% were also removed (rs4994-ADRB3, rs7913948-ALOX5, rs7412-APOE, rs328-LPL, rs16139-NPY, rs6008259-PPARγ, rs2066826-PTGS2) (Supplementary Table 1). Using the ‘n/15 rule’[36], only 21 predictors could be used to get reliable estimates. Consequently, data reduction was applied to the SNPs by applying Multiple Correspondence Analysis (MCA) and using only the first 15 MCA dimensions. Ethnic groups representing less than 1% each were merged into the ‘Other’ category.

Wald tests for all the predictors in the model were generated, and they were further ranked by importance based on the χ2-degrees of freedom (-df) score. A simplified model was obtained by applying a ‘fast-backwards’ variable elimination approach[37] based on the Akaike’s Information Criterion (AIC)[38]. The approximate βs and ORs (and their 95% confidence intervals, CI) of the remaining variables were reported. The full model was validated through bootstrap to provide estimates of the performance of the model in new data in comparison with the training data, in the form of Sommers Dxy, R2, intercept, slope, Emax, and Briers B score. The significant contributions (coefficients of determination R2 with p values < 0.05) of SNPs to the essential MCA dimension were plotted to deconvolute it.

We used the proportional odds assumption in the model. ‘Physical activity’, ‘KIDMED’ and ‘Urolithin metabotype’ showed some deviation from this assumption. However, alternative extended continuation ratio models with this assumption relaxed for these variables did not result in better models according to the AIC, i.e., the higher complexity of the model was not compensated by the increase in the fit. The inclusion of transformations of some predictors, including a restricted cubic spline for both ‘Age’ and ‘KIDMED’, did not result in improved models, as judged by the AIC. Finally, a genetic score (computed as the sum of risk alleles) was also tested as a possible surrogate for the SNPs variables, but this did not result in a better replacement for the MCA dimensions. All tests were bilateral, with a significance level of 0.05.

Other statistical analyses were carried out using the SPSS software, v23.0 (SPSS Inc., Chicago, IL, USA). When more than two groups were compared, analyses of variance (ANOVA), followed by Bonferroni-corrected t-test (for post-hoc analysis) or the Kruskal–Wallis followed by Dunn’s test were used for normally and non-normally distributed data, respectively (KIDMED score vs FTO genotype TT, AT or AA, etc.). Comparison of non-normally distributed quantitative variables between two clusters was approached using the Mann-Whitney U-test (FTO TT genotype vs BMI or waist, etc.). Comparison of categorical variables was assessed using the Pearson’s χ2 test. Spearman’s rank or Pearson correlations were applied to explore possible associations between variables (BMI vs hip/height, etc.). Plots of data were performed using Sigma Plot 13.0 (Systat Software, San Jose, CA, USA).

Results

Characteristics of the cohort and associations with overweight-obesity

Table 1 shows characteristics of the cohort as a function of age and sex, including anthropometric values (hip, waist, weight and BMI), the distribution of urolithin metabotypes (A, B and 0), the KIDMED scores (grouped as low, medium and good adherence), physical activity and the percentage of normoweight, overweight and obesity. The participants were mainly Caucasian-Europeans (93.5%), with a small proportion of Arabs (2.9%) and Amerindians (2.2%), and a marginal presence of Black-Africans (0.96%), Asian-Chinese (0.22%), and Indo-Aryans (0.22%).

Table 1

Characteristics of the study population.

Age (y)	Sex	Hip (cm)*	Waist (cm)*	Weight (kg)*	Height (cm)*	BMI (kg/m2)*	BMI z-score*	NW (%)	OW (%)	OB (%)	KIDMED (%) (L/M/G)	Physical Activity (%) (H/L)	Metabotype (%) (A/B/0)
5	F (n = 16)	65.5 (58.0–81.0)	58.5 (48.0–76.0)	20.9 (15.9–31.1)	110.5 (103.5–119.0)	16.4 (13.6–25.7)	0.8 (−1.4–2.9)	50.0	18.8	31.2	28.6/64.3/7.1	76.9/23.1	81.3/0.0/18.7
	M (n = 14)	66.5 (58.0–83.0)	59.5 (51.0–71.5)	22.2 (17.4–34.2)	114.5 (110.0–124.0)	16.7 (13.7–24.4)	0.9 (−1.8–3.5)	50.0	21.4	28.6	57.1/28.6/14.3	84.6/15.4	64.3/0.0/35.7
6	F (n = 10)	69.5 (59.0–81.0)	59.3 (52.0–76.0)	23.7 (16.8–37.1)	118.9 (110.0–128.5)	15.9 (13.9–24.9)	0.4 (−1.1–2.6)	60.0	10.0	30.0	30.0/50.0/20.0	60.0/40.0	90.0/10.0/0.0
	M (n = 9)	68.0 (64.0–78.0)	59.0 (55.0–71.0)	22.7 (20.8–32.5)	121.5 (116.0–125.5)	16.3 (14.5–21.9)	0.7 (−0.8–2.5)	55.5	0.0	45.5	11.2/44.4/44.4	77.8/22.2	77.8/11.1/11.1
7	F (n = 9)	70.0 (66.0–82.0)	64.0 (55.0–76.0)	26.8 (24.5–40.5)	130.0 (123.5–132.9)	17.6 (15.0–24.0)	1.0 (−0.3–2.3)	55.5	11.1	34.4	22.2/44.4/33.3	88.8/11.2	55.5/33.3/11.2
	M (n = 16)	70.5 (65.0–85.0)	59.5 (55.0–76.0)	26.6 (21.1–44.5)	123.8 (116.0–140.6)	17.0 (13.9–24.6)	0.9 (−1.4–2.6)	56.3	12.5	31.2	33.4/53.3/13.3	93.8/6.2	68.8/18.8/12.4
8	F (n = 16)	74.5 (67.0–91.0)	60.5 (53.0–80.0)	32.1 (24.1–48.5)	132.5 (124.1–142.0)	18.3 (15.1–26.4)	1.0 (−0.4–2.4)	50.0	18.8	31.2	6.7/86.6/6.7	56.3/43.7	93.8/0.0/6.2
	M (n = 16)	75.0 (62.0–97.0)	63.0 (53.0–85.0)	28.4 (19.8–63.1)	128.2 (118.3–153.5)	17.6 (11.4–26.8)	0.9 (−5.7–2.5)	50.0	18.8	31.2	37.4/56.3/6.3	73.3/26.7	81.3/0.0/18.7
9	F (n = 12)	77.0 (68.0–103.0)	64.0 (53.0–82.0)	36.9 (27.9–63.2)	141.6 (134.9–149.8)	18.4 (14.5–30.0)	0.8 (−1.1–2.5)	58.3	16.7	25.0	41.7/50.0/8.3	100.0/0.0	100.0/0.0/0.0
	M (n = 12)	80.6 (71.0–97.0)	70.0 (58.0–90.5)	38.5 (27.9–59.6)	137.3 (125.7–152.2)	20.1 (16.1–28.3)	1.2 (0.0–2.4)	33.3	41.7	25.0	33.3/66.6/0.0	91.7/8.3	75.0/16.7/8.3
10	F (n = 14)	79.0 (74.0–99.0)	65.5 (57.0–82.2)	38.5 (30.4–61.7)	143.0 (137.8–159.5)	19.2 (15.6–27.3)	0.8 (−0.6–2.2)	50.0	28.6	21.4	35.7/50.0/14.3	57.1/42.9	78.6/14.3/7.1
	M (n = 15)	87.0 (73.0–98.0)	71.0 (60.0–90.0)	46.7 (30.5–66.4)	146.7 (132.6–152.5)	20.8 (16.1–30.7)	1.4 (−0.3–2.5)	26.7	33.3	40.0	40.0/53.3/6.7	93.3/6.7	86.7/13.3/0.0
11	F (n = 20)	86.5 (74.0–106.0)	65.0 (58.0–83.0)	46.4 (32.4–74.3)	153.0 (143.0–173.2)	19.3 (15.5–28.0)	0.6 (−1.0–2.1)	65.0	25.0	10.0	10.0/80.0/10.0	60.0/40.0	95.0/5.0/0.0
	M (n = 17)	88.0 (69.0–101.0)	73.0 (58.0–85.0)	47.0 (29.6–63.2)	149.0 (138.0–161.0)	20.8 (15.2–27.4)	1.2 (−1.2–2.1)	35.3	35.3	29.4	11.7/82.4/5.9	76.5/23.5	94.1/0.0/5.9
12	F (n = 21)	88.0 (79.0–108.0)	66.0 (60.0–84.0)	48.8 (36.7–72.6)	156.5 (146.0–165.2)	20.4 (16.5–29.5)	0.7 (−0.8–2.1)	66.7	19.0	14.3	33.3/52.4/14.3	71.4/28.6	61.9/23.8/14.3
	M (n = 20)	90.0 (74.0–107.0)	72.0 (62.0–93.0)	54.1 (33.9–77.5)	156.5 (142.9–170.5)	22.2 (16.1–28.8)	1.3 (−0.9–2.2)	25.0	50.0	25.0	30.0/60.0/10.0	80.0/20.0	85.0/10.0/5.0
13	F (n = 14)	89.0 (75.0–123.0)	67.5 (58.0–100.0)	51.2 (35.0–92.9)	154.6 (145.3–167.1)	20.8 (14.6–38.4)	0.6 (−2.2–2.6)	71.4	7.1	21.5	21.4/78.6/0.0	50.0/50.0	57.1/35.8/7.1
	M (n = 17)	91.0 (81.0–128.0)	71.0 (62.0–100.0)	59.1 (42.2–112.6)	165.6 (152.1–180.2)	20.4 (15.1–38.5)	0.7 (−2.0–2.3)	58.8	5.9	35.3	29.4/70.6/0.0	81.3/18.7	82.4/17.6/0.0
14	F (n = 19)	90.0 (81.5–106.0)	69.5 (58.0–90.0)	53.2 (43.4–73.5)	160.5 (151.7–171.5)	19.3 (16.1–29.6)	0.0 (−1.5–1.9)	68.4	26.3	5.3	36.8/63.2/0.0	36.8/63.2	78.9/15.8/5.3
	M (n = 11)	94.0 (71.0–123.0)	72.0 (67.0–101.0)	61.7 (46.0–106.9)	171.5 (162.5–189.5)	20.9 (17.1–36.1)	0.6 (−1.0–2.5)	63.6	27.3	9.1	9.1/90.9/0.0	81.8/18.2	81.8/0.0/18.2
15	F (n = 15)	88.0 (82.0–99.0)	68.0 (59.0–81.0)	52.5 (43.5–63.4)	158.5 (156.0–175.5)	20.3 (16.9–24.5)	0.1 (−1.3–1.1)	93.3	6.7	0.0	40.0/60.0/0.0	42.8/57.2	86.8/6.6/6.6
	M (n = 23)	94.0 (80.0–112.0)	75.5 (60.0–96.0)	62.6 (43.5–93.3)	173.1 (156.0–184.5)	21.2 (17.9–30.3)	0.5 (−0.9–2.1)	78.3	8.7	13.0	52.2/39.1/8.7	52.2/47.8	82.6/8.7/8.7
16	F (n = 23)	94.0 (82.0–138.5)	70.0 (61.0–137.0)	57.1 (41.0–146.1)	157.8 (143.5–172.2)	23.1 (17.5–49.3)	0.7 (−1.3–2.7)	73.9	21.7	4.4	34.8/52.2/13.0	39.1/60.9	65.2/21.7/13.1
	M (n = 22)	97.0 (83.0–121.0)	76.0 (64.0–100.0)	71.1 (49.4–110.8)	176.0 (164.2–186.5)	22.8 (18.2–35.8)	0.7 (−1.1–2.5)	63.6	27.3	9.1	22.7/68.2/9.1	77.3/22.7	77.3/4.5/18.2
17	F (n = 19)	91.5 (84.0–110.0)	72.0 (59.0–87.0)	60.6 (50.7–76.1)	163.0 (149.5–170.0)	23.0 (17.9–28.0)	0.6 (−1.3–1.5)	73.7	26.3	0.0	0.0/84.2/15.8	21.1/78.9	84.2/10.5/5.3
	M (n = 15)	92.0 (85.0–111.0)	73.5 (64.0–91.0)	66.5 (50.5–96.0)	176.2 (160.5–186.8)	22.2 (17.8–30.0)	0.3 (−1.2–1.9)	80.0	13.3	6.7	46.7/40.0/13.3	73.3/26.7	66.7/26.7/6.6

*Results are expressed as median and (range); F, Female; M, Male; KIDMED (L/A/G): Low quality, Medium quality and Good quality diet, respectively. NW, normoweight; OW, overweight; OB, obesity; Physical activity: H, high; L, low.

The hip-to-height ratio was the best anthropometric index associated with BMI (r = 0.78, p = 1.3 × 10–86) vs the waist-to-hip ratio (r = 0.13, p = 0.007) and the waist-to-height ratio (r = 0.32, p = 0.001) (Supplementary Fig. 3). As expected in growing children, BMI values increased on average from 5 to 17 years (Fig. 1A). The percentage of overweight-obesity decreased from 5 to 17 years (from 50% to 25%, respectively) with the exemption of boys from 9 to 12 years old (n = 64) that reached the highest prevalence of overweight-obesity (∼70%) (Fig. 1B). This categorized representation (Fig. 1B) was coincident with the quantitative evolution of BMI z-scores for girls (Fig. 1C) and boys (Fig. 1D) as a function of age.

Figure 1

Distribution (%) of (A) BMI, (B) overweight–obesity (OW–OB), (C,D) BMI z-scores, and (E) urolithin metabotypes in the cohort (n = 415) from 5 to 17 years. (F) Distribution of urolithin metabotypes (UM-A, UM-B and UM-0) in normoweight (NW) and OW + OB children. The shadow area in C and D designates overweight-obesity according to WHO (https://www.who.int/growthref/who2007_bmi_for_age/en/). ★, median of the group for each age.

Regarding the SNPs analysed, after correcting for multiple tests, three SNPs did not satisfy the Hardy-Weinberg equilibrium (HWE) (rs1801253, rs5082, rs11868035) (Supplementary Table 1). The rest of the SNPs were in equilibrium and were close to European frequencies. In the case of these three SNPs, we can speculate that the reason could be the association established by the presence of certain consanguinity (several sibling groups in the cohort). Nevertheless, this did not affect our results, since we aimed to compare variables from different domains and not only in estimating the particular effect of a single SNP. Overall, it is not absolutely necessary to have HWE in our approach, i.e., to estimate odds ratios of an MCA dimension and rank the predictors.

The distribution of urolithin metabotypes in this cohort (Fig. 1E) revealed that both a lower prevalence of UM-A and a higher occurrence of UM-0 were associated with an increased percentage of overweight-obesity after bivariate analysis (p = 0.015) (Fig. 1F).

No significant association between physical activity and overweight-obesity distribution was found (results not shown). On the contrary, there were many overweight-obese children with high physical activity, which should be explained as a consequence of their overweight-obesity status (results not shown). Regarding the diet, the KIDMED scores showed mean values of 6.9 ± 2.1 for the entire cohort (Supplementary Table 2) and ranged from the lowest value of 5.7 ± 1.4 in 15-year-old girls to the highest value of 8.8 ± 2.2 in 6-year-old boys (Fig. 2A). However, no significant differences were found between boys and girls as well as through the range of age (Fig. 2A). Besides, no significant association was found between overweight-obesity and KIDMED in this group (results not shown).

Figure 2

(A) KIDMED scores in girls (•) and boys (▴) from 5 to 17 years. Results are expressed as mean ± SD. (B) KIDMED scores (box plots) as a function of the genotypes (TT, AT and AA) of the rs9939609 SNP in the FTO gene. Significant differences are shown after the Kruskal Wallis test. Bar charts show the proportion of normoweight (NW) and overweight-obesity (OW + OB) in the cohort depending on the rs9939609-FTO genotype of the individuals.

We next explored the potential association of both adherence to the Mediterranean diet and percentage of overweight-obesity with the well-known obesity risk allele A of the rs9939609 SNP in the FTO gene. Figure 2B shows the KIDMED scores and the percentage of normoweight and overweight-obesity in children depending on their genotype TT, AT or AA. Remarkably, the adherence to the Mediterranean diet was significantly lower in children with the risk-associated genotypes AT and AA vs the TT genotype (Fig. 2B). Although the proportion of overweight-obesity was the highest for the AA genotype (45.5%), however, the difference vs that of the TT genotype (38.7%) did not reach statistical significance (Fig. 2B). Once again, many exemptions prevented the usefulness of the combination of KIDMED scores and rs9939609 genotypes as unique predictors of overweight-obesity in this cohort.

Therefore, although both the urolithin metabotypes and rs9939609 SNP-FTO could partially contribute to the overweight-obesity distribution as independent predictor variables in this cohort of children and adolescents; however, all the possible SNP-SNP interactions together with the rest of variables had not been taken into account. Therefore, we next developed an ordinal logistic model to identify the consortium of variables that could estimate the odds ratios of the overweight-obesity distribution in this cohort.

An ordinal logistic model to identify the consortium of variables associated with overweight-obesity

In this model, we used as the ordinal response the normoweight, overweight and obesity classification for children, based on sex and age to estimate the odds ratios (and the corresponding 95% confidence intervals) for different predictor variables, as well as their relative importance in the prediction of the response. In this holistic approach, we used as predictors the ‘Ethnicity’, ‘Urolithin metabotypes’, ‘KIDMED score’, ‘Physical activity’, and genetic polymorphisms (44 SNPs were finally included and compressed into 15 MCA dimensions, hereafter termed ‘SNP.Dim.’), together with sex and age. Although the distribution of normoweight, overweight and obesity WHO-based categories was apparently adjusted by sex and age, we still observed a trend for decreasing average age when moving from normoweight, overweight and obese children, as well as enrichment in boys in the same order. Therefore, we also included these variables in the model.

Figure 3 displays the predictors used in the ordinal logistic model ranked by their apparent importance, as measured by their χ2-df score. The model was highly significant (p < 0.0001), yielding the components ‘Age’, and ‘SNP.Dim.14’ as the two most important predictors. Therefore, in the prediction of overweight-obesity in this population, apart from ‘Age’ and ‘Sex’, the most critical contributing variables were SNPs (through the variable SNP.Dim.14), followed by ‘Ethnicity’, ‘Urolithin metabotype’ and ‘KIDMED’. The ‘Physical activity’ seemed to be irrelevant in this sample.

Figure 3

Apparent predictor importance in the ordinal logistic model. Variables are ranked by relevance (based on their χ2-df score). The χ2 of the Wald test for each predictor is also shown.

We next applied a fast-backwards approach based on AIC to obtain a reduced model and estimate the corresponding βs and ORs (as well as their 95% CI). The variables kept were ‘Sex’, ‘Age’, ‘Urolithin metabotype’, ‘KIDMED’, and the genetic components ‘SNP.Dim.3’, ‘SNP.Dim.11’, and ‘SNP.Dim.14’ (Table 2). The variable ‘Ethnicity’ was not retained in the simplified model, probably due to its high complexity (4 levels). Our results reveal that being a boy either with UM-B or UM-0 and having a higher SNP.Dimension.14, all increased the chances of overweight-obesity in this study population. On the contrary, ageing, better adherence to the Mediterranean diet (higher KIDMED score), and being UM-A was associated with lower probabilities of being overweighed-obese. Remarkably, all these variables operated additively to build the final probability of overweight-obesity for each subject.

Table 2

Beta estimates (βs), odds ratios (ORs), and their corresponding confident intervals (CIs).

Variables*	βs (CIs)	ORs (CIs)
Sex (boy)	0.51 (0.10, 0.918)	1.67 (1.11, 2.5)
Age	−0.12 (−0.17, −0.0664)	0.88 (0.83, 0.93)
UM-B	0.21 (−0.40, 0.837)	1.24 (0.66, 2.31)
UM-0	0.68 (0.02, 1.34)	1.98 (1.03, 3.81)
KIDMED	−0.06 (−0.16, 0.0256)	0.93 (0.85, 1.03)
SNP.Dim.3	1.05 (0.03, 2.07)	2.86 (1.03, 7.94)
SNP.Dim.11	1.18 (−0.01, 2.37)	3.25 (0.98, 10.7)
SNP.Dim.14	2.32 (1.08, 3.55)	10.1 (2.94, 35)

*SNP.Dim, dimensions (consortium) of SNPs obtained after Multiple Component Analysis; UM-B and UM-0, urolithin metabotype B and 0, respectively.

Figure 4 shows the R2 values of the 24 SNPs that significantly contributed to SNP.Dim.14 (p < 0.05). The most contributing SNP to this consortium was rs1801253-ADRB1, as well as other obesity risk-associated SNPs (UCP2, ADIPOQ, LEP, MC4R, etc.). However, there were other SNPs with less known involvement in obesity. Table 3 shows all the SNPs contributing to SNP.Dim.14 with their definitions and main processes in which they are involved.

Figure 4

Description of components of SNP Dimension-14. Significant R2 (p < 0.05) for each SNP within the SNP Dimension-14 are shown.

Table 3

The consortium of significant single-nucleotide polymorphisms (SNPs) integrated into the MCA dimension SNP.Dim.14.

R2*	SNP	Gene	Location	(Positive) contributing genotypes	Name	Gene-related process	Gene-related disease
0.294	rs1801253	ADRB1	Coding	C/C	Adrenoceptor β1	Lipolysis and fat oxidation	CVD#
0.096	rs4343	ACE	Coding	A/A A/G	Angiotensin I converting enzyme	Signalling, metabolism	CVD, MetS, Alzheimer’s
0.091	rs8061518	FTO	Intronic	G/G A/G	Fat mass and obesity-associated	Metabolism	Obesity, CVD, MetS, diabetes
0.076	rs1130864	CRP	Intronic	G/G	C-Reactive Protein	Inflammation	CVD
0.067	rs659366	UCP2	Intronic	C/C	Uncoupling protein 2	Metabolism	Obesity
0.058	rs6131	SELP	Coding	C/C	Selectin P	Cell adhesion	CVD
0.057	rs12535708	LEP	Intronic	C/C A/C	Leptin	Signalling, metabolism	Obesity
0.041	rs1501299	ADIPOQ	Intronic	T/T	Adiponectin	Metabolism	MetS, diabetes, cancer
0.036	rs708272	CETP	Intronic	A/A	Cholesteryl ester transfer protein, plasma	Metabolism	CVD
0.033	rs2241766	ADIPOQ	Coding	T/T	Adiponectin	Metabolism	MetS, diabetes, cancer
0.033	rs696217	GHRL	Intronic	T/T G/T	Ghrelin/obestatin prepropeptide	Signalling, metabolism	Obesity, cancer
0.031	rs11868035	SREBF1	Intronic	G/G	Sterol regulatory element-binding transcription factor 1	Metabolism	Diabetes
0.027	rs17782313	MC4R	Intronic	T/T	Melanocortin 4 receptor	Signalling, metabolism	Obesity, diabetes, cancer
0.023	rs1801282	PPARγ	Intronic	C/C	Peroxisome proliferator-activated receptor-γ	Signalling, metabolism	Diabetes, cancer
0.020	rs5443	GNB3	Coding	C/C	(G protein), β-polypeptide 3	Signalling, metabolism	Obesity, CVD, diabetes
0.019	rs3758538 (3944 A > C)	RBP4	Intronic	G/G G/T	Retinol binding protein 4	Inflammation	Obesity
0.018	rs1801133	MTHFR	Coding	G/G	Methylenetetrahydrofolate reductase	Metabolism	CVD, cancer
0.017	rs1143634	IL1B	Coding	A/A A/G	Interleukin 1β	Inflammation	Diabetes, periodontitis, cancer
0.016	rs693	APOB	Coding	G/G	Apolipoprotein B	Metabolism	MetS, diabetes, cancer
0.014	rs894160	PLIN1	Intronic	C/C	Perilipin 1	Metabolism	Obesity
0.013	rs1799883	FABP2	Coding	C/C	Fatty acid-binding protein 2	Metabolism	Diabetes, MetS
0.013	rs662799	APOA5	Intronic	G/G A/G	Apolipoprotein A5	Triglycerides	CVD
0.013	rs429358	APOE	Coding	C/C C/T	Apolipoprotein E	Metabolism	Alzheimer’s, CVD
0.010	rs9930333	FTO	Intronic	T/T	Fat mass and obesity-associated	Metabolism	Obesity, CVD, MetS, diabetes

*R2, Contribution of the SNPs to Dim.14 according to Fig. 3. #CVD, cardiovascular disease; MetS, metabolic syndrome.

Model validation

We used bootstrap to validate the model (Table 4). The model showed some degree of overfitting, reflected in a decrease of the indexes after correction for optimism. The discriminative capacity of the model was modest, although not negligible according to the Dxy, R2 indexes (it must be taken into account that the typical values of R2 in ordinal models are much lower than those observed in linear regression models), and the Brier score of 0.23. The calibration showed some degree of shrinkage, as seen from the deviation of the intercept and slope from 0 and 1, respectively. Emax, the maximum calibration error in predicting p (Y > normal), was 0.11, showing some degree of miscalibration. Overall, we obtained a significant model that still showed predictive power in external data.

Table 4

Validation of the model using Bootstrap.

	Index-orig	Training	Test	Optimism	Index-corrected
Dxy*	0.4136	0.4688	0.3518	0.1170	0.2967
R2	0.1824	0.2418	0.1323	0.1095	0.0729
Intercept	0.0000	0.0000	−0.1099	0.1099	−0.1099
Slope	1.0000	1.0000	0.6633	0.3367	0.6633
Emax	0.0000	0.0000	0.1110	0.1110	0.1110
B	0.2097	0.1978	0.2222	−0.0243	0.2340

*Dxy = bias-corrected Somers rank correlation coefficient that goes from −1 to 1; R2 = Nagelkerke R2 that goes from 0 to 1; Intercept and Slope of a logistic calibration equation (should be 0 and 1, respectively, for a perfect fit); Emax = maximum calibration error for p (Y = 0) based on the linear-logistic recalibration; B = Brier’s quadratic probability score, which goes from 0 (the best score) to 1 (the worst score).

Discussion

The present study shows that the probabilities of childhood overweight-obesity cannot be explained by specific isolated variables (that is, one or a few SNPs, or just physical activity, or diet, etc.), but by complex, multifactorial associations of environmental and genetic components. Besides, a genetic score, widely used in the literature[39], and understood as the theoretical equal-additive contribution of every single SNP could not explain the distribution of overweight-obesity in our cohort. In this regard, our analysis adds significance to previous studies that have associated the genetic background of adult individuals with the adherence to the Mediterranean diet and some markers related to obesity and metabolic syndrome[39-41]. However, to the best of our knowledge, there are no previous studies that explore the occurrence of overweight-obesity in children and adolescents by estimating the odds ratios for the predictor variables ‘urolithin metabotypes’ as gut microbiota biomarkers, ‘age’, ‘sex’, ‘adherence to the Mediterranean diet’, and an identified consortium of 24 SNPs from 22 genes, mainly related to obesity and cardiometabolic diseases.

The purpose of the ordinal logistic regression model was to estimate the odds ratios of the different explanatory variables and rank them by importance, rather than a tool to predict the probability of overweight and obesity of new individuals. There are in the literature several previous reports that model the probabilities of overweight and obesity in children through ordinal regression models[42-46]. However, to the best of our knowledge, this is the first time that this type of model is used with such a wide set of predictors comprising different putatively influencing domains (genetic, diet, exercise, microbiome urolithin metabotype, and ethnicity, in addition to sex and age) and ranks them according to their relative importance. The χ2-df score returned the following importance ranking: age >SNP.Dim.14 >ethnicity >sex >urolithin metabotype >KIDMED; while physical activity seemed to have a negligible importance. The use of an AIC-based fast-backwards reduction of the model removed ‘Ethnicity’ from the final model, probably because of the large number of categories in this predictor[42], that would be highly penalised by the AIC criterion. Besides, our sample was mainly of Caucasian-European origin (93.5%), and it is expected that in a sample with a more balanced distribution of ethnicity, this predictor would have high importance. Furthermore, the modest predictive capability of the model suggested the need for additional predictors, which is not unexpected, given the multifaceted aetiology of obesity. Still, the model, derived for estimating purposes, allowed us to rank predictors of different domains by their importance, as well as to estimate odds ratios for them.

The classification of overweight and obesity in adults is rather simple and independent of age and sex (i.e. 25 ≤ BMI < 30 for overweight, and BMI ≥ 30 for obesity). However, the WHO-based definition for overweight-obesity in children and teens uses age- and sex-based percentiles[5], and a higher BMI does not necessarily correspond with overweight-obesity in growing children. Therefore, using overweight and obesity as response variables instead of BMI values is especially useful in children and adolescents. This is due to the difficulty in comparing them at different ages and of the two sexes, given the remarkable, sex-dependent change in weight and height, during this period of development (i.e., while our model predicted that overweight-obesity decreases upon ageing, the use of BMI yielded the opposite result, masking the real result). Still, sex and age were essential predictors in our model, with estimated ORs of 1.67 for boys, and 0.88 for a one-year increase, respectively, which means an increased percentage of overweighed-obese boys, and a reduction of overweight-obesity upon ageing.

The gut microbiota contributes to the pathophysiology of obesity[3], and a recent report shows the role of gut microbial metabolites in the expression of the microRNA-181 family, which regulates white adipose tissue inflammation and obesity in children[14]. However, the gut microbiota, as a predictor of childhood obesity, has been scarcely approached and not usually considered together with SNPs and other variables. In a targeted approach, among a few analysed microbial groups, Scheepers et al. reported that the Bacteroides fragilis group was associated with childhood weight development[47]. Later, Rampelli et al. connected the onset of obesity in children to an increase of Proteobacteria and a decrease of Clostridiaceae and Ruminococcaceae[48]. In the same line, Bai et al. described that the Proteobacteria phylum had significantly enriched OTUs for higher BMI levels in a cohort of children from the American Gut Project[49]. This is relevant since the microbiota-associated with UM-B is enriched in Proteobacteria[25]. Interestingly, Nirmalkar et al. reported that the Coriobacteriaceae family was 3-fold more abundant in obese children and adolescents than normoweight[50]. Specifically, the genus Collinsella was more abundant in obese adolescents, a microbial group that could be related to the endothelial dysfunction[50,51]. Overall, all these results are remarkable since we recently reported that UM-B was enriched in the Coriobacteriaceae family, which was positively correlated with blood total-cholesterol, LDL-cholesterol, and BMI in adults[25]. Although we had previously observed a trend between UM-B occurrence and overweight-obesity in adult individuals, however, we could not establish a definitive link probably due to the lack of other interacting factors such as those included in the present study (SNPs, diet, etc.)[26,52]. Therefore, the Coriobacteriaceae family, and probably the Proteobacteria phylum, more abundant in obese children as well as in UM-B, could be the link between UM-B occurrence and overweight-obesity, which suggests that the microbiota associated with urolithin metabotypes could contribute in the prediction of the probability of being overweighed or obese. Regarding UM-0, its occurrence in the population is approximately constant (~10%), although there is both a higher occurrence and variability in the childhood[26]. The microbiota associated with UM-0 has been reported to show lower diversity than UM-B and UM-A, which could be indicative of an obesity-prone microbiota[4]. Nevertheless, we cannot exclude a possible shift of metabotype in children from UM-0 to either UM-A or UM-B determined by ageing[26] or after ellagitannin-rich diets as previously described in adults[21].

Regarding the possible health implications of belonging to one or another metabotype, we have recently reported that urolithin metabotypes determined the different restoration capacity of the gut microbiota and the anthropometric values (weight, waist and hip) of healthy women up to 12 months after delivery[53]. The gut microbiota of pregnant women is in dysbiosis, which persists at least 1 month after delivery. We observed that the gut microbiota associated with UM-B was more resilient than that of UM-A, which would have negative implications in the dysbiotic-prone UM-B. In contrast, the gut microbiota of UM-A women progressively became normal during the year after childbirth. Therefore, we suggested that the determination of urolithin metabotypes in pregnant and lactating women could be a useful tool to predict their predisposition to the recovery of the gut microbiota and anthropometric values, significantly altered during pregnancy and after childbirth[53].

Many studies describe the association, or lack of association, of specific SNPs with obesity. The rs9939609 SNP-FTO has been reported to confer a predisposition to obesity by regulating the control of food intake and food choice, suggesting a link to a hyperphagic phenotype or a preference for energy-dense foods in Scottish children[54]. This agrees with the connection between the low adherence to the Mediterranean diet and the risk allele A of rs9939609 SNP-FTO in our cohort. However, we did not observe a clear association of this SNP with overweight-obesity, but only a trend with many exemptions, in agreement with other studies where some SNPs, previously reported to be involved in obesity, such as the rs17782313-MC4R and rs9939609-FTO, exerted weak effects and very scarce contribution to obesity in 773 pre-pubertal Portuguese children[55]. In the present study, a consortium of 24 SNPs was identified as the second contributing predictor to overweight-obesity in our cohort. The rs1801253 SNP (also called Arg389Gly), located in the gen ADRB1, was the most contributing SNP within this consortium. It is known the involvement of the ADRB1 gene polymorphisms in cardiovascular diseases[56] but also in obesity[57,58].

The ADRB1 gene codifies a G-coupled protein (the β1-adrenergic receptor) that binds the catecholamines epinephrine and norepinephrine and controls sympathetic responses in the heart, kidney and adipocytes. Interestingly, Dionne et al. reported that the rs1801253 SNP-ADRB1 was associated with higher body weight and BMI in a cohort of Caucasian women (n = 931)[57]. In the same line, Aradillas-García et al. reported that the rs1801253 SNP-ADRB1, but not the Trp64Arg ADRB3, was associated with obesity in Mexican children[59]. The connection between ADBR1 SNPs and obesity could rely on the catecholamines, which are considered significant lipolysis regulators[60] and affect differentiation and proliferation of adipocytes[61]. In this regard, Lee et al.[62] established the association between impaired urinary epinephrine and norepinephrine excretion and obesity, insulin resistance, and metabolic syndrome in a cohort of 577 Chinese subjects.

Overall, all the above highlights again the need to consider SNPs consortia, interacting with other variables, instead of few SNPs in those studies aimed to associate SNPs with obesity. We are aware that the present study is an exploratory validation for a proof-of-concept, i.e., an ordinal logistic model that associates child overweight-obesity with a consortium of SNPs potentially interacting with the urolithin metabotypes-associated microbiota, adherence to the Mediterranean diet, age, and sex. Although we claim for the rationale of our approach and its potential usefulness, however, our results should be confirmed with additional research. We also acknowledge some limitations that should be considered in further studies, which also could improve its prediction capability. For example, it would be interesting to include other possible variables, such as the detailed composition and functionality of the individuals’ gut microbiomes, and dietary interventions to evaluate not only associations but also individuals’ responses. The latter would be even better than the use of validated questionnaires. Besides, a higher number of SNPs (or many SNPs associated with a specific gene) should be explored, especially in children from other geographical origins and ethnicities. Also, the inclusion of serobiochemical variables and traits related to obesity and its comorbidities (blood lipid profile, blood pressure, glucose homeostasis, etc.) could yield relevant information. Finally, a validation cohort (i.e., a parallel-group with all the children either normoweight or obese) should also be considered in further studies to confirm our model fully.

Conclusions

The present research highlights the need for a holistic approach to unravel the predictors of overweight-obesity in children. Our results confirm, in agreement with the multifaceted aetiology of obesity, the link of childhood overweight-obesity to multifactorial associations of environmental and genetic components. The ordinal logistic model revealed that child overweight-obesity prevalence was related to being a young boy with either UM-B or UM-0, low KIDMED score and high contribution of a consortium of 24 SNPs, being rs1801253-ADRB1, rs4343-ACE, rs8061518-FTO, rs1130864-CRP, rs659366-UCP2, rs6131-SELP, rs12535708-LEP, rs1501299-ADIPOQ, rs708272-CETP and rs2241766-ADIPOQ the top-ten contributing SNPs. Therefore, it is of particular relevance the evaluation of interactive SNPs consortia along with the stratification of the children according to their urolithin metabotypes, which could be early biomarkers, in the case of UM-B and UM-0, of a dysbiotic-prone obesity-associated microbiota.

Supplementary information.