To evaluate sex-specific effects of genetic variants within the genes of interest (CKB, CKMT1B, and GATM), data from a sex-stratified meta-analysis of GWAS for BMI (Pulit et al., 2019) were utilized (https://zenodo.org/record/1251813#.YbsPznyZNaR; downloaded on 8 June 2022). A total of 434,794 female and 374,756 male samples as well as 27, 381, 302 SNPs were included (Pulit et al., 2019). Information regarding the genomic position, rsID, minor allele frequency (MAF), standard error (SE), beta and p-values (p) of the SNPs located within the genes and a 500 kb region adjacent to the genes were extracted from the datasets. Given the spatial proximity of GATM and CKMT1B, these genes were analyzed collectively within one region (see Table 1). After assessing whether the effect alleles for males and females were homogeneously assigned, palindromic SNPs with a MAF between 0.48 and 0.52 were excluded. The SNPs’ effect directions were subsequently harmonized to eliminate strand issues by inverting the beta values of SNPs tested for their major allele. Next, z-scores were calculated for each SNP (see Eq. 1; Khramtsova et al., 2019; Zheng et al., 2022). Furthermore, additional data from a GWAS for AN (Watson et al., 2019) were obtained (https://www.med.unc.edu/pgc/download-results/; downloaded on 21 October 2021). However, as no information on MAF was available in this dataset, effect size harmonization was omitted. SNPs were considered trait-associated if they exhibited a GWAS-derived p < 5*10^–8. z − s c o r e = b e t a f e m a l e − b e t a m a l e S E 2 f e m a l e + S E 2 m a l e (1)

SNPs were defined to have a sexually dimorphic effect if they exhibited an absolute z-score of at least 3 (|z-score| ≥ 3), as described previously (Khramtsova et al., 2019).

We sequenced the DNA of a screening group consisting of 192 children and adolescents with severe obesity (age: 13.88 ± 2.50 years; BMI: 32.36 ± 6.04 kg/m²; all had a BMI percentile >90th; 182 had a BMI >97th percentile; see Table 2; Hinney et al., 2007; Scherag et al., 2010), 192 patients with anorexia nervosa (acute or recovered; age: 21.85 ± 10.60 years; BMI: 16.38 ± 3.12 kg/m²; see Table 2; Boraska et al., 2012; Watson et al., 2019) and 192 healthy-lean controls (age: 25.55 ± 3.91 years; BMI: 18.38 ± 1.10 kg/m²; see Table 2; Hinney et al., 2007). Additional genotyping was performed in a confirmation study group of 781 families comprising a child or adolescent with severe obesity and both biological parents (referred to as trios; age of offspring: 13.45 ± 3.11 years; BMI of offspring: 32.15 ± 5.92 kg/m²; see Table 2) from a family-based obesity GWAS, 320 children and adolescents with severe obesity (age: 14.60 ± 4.24 years; BMI: 33.77 ± 7.00 kg/m², see Table 2) and 253 healthy-lean controls (age: 26.40 ± 6.84 years; BMI: 17.85 ± 1.11 kg/m²; see Table 2) from a case-control GWAS for obesity (Hinney et al., 2007; Scherag et al., 2010). The 781 trios were composed of a child or adolescent with severe obesity (BMI percentile 90th) and both biological parents. Over 90% of the offspring (723; 92.57%) had a BMI at or above the 97th BMI percentile. The mean BMI of the 1,562 parents was 30.40 kg/m² (standard deviation: 6.3 kg/m²). Of these, 724 (46.35%) had a BMI above 30 kg/m², while 822 (52.62%) showed a BMI lower or equal 30 kg/m². For 16 parents, no BMI data was available. Written informed consent was given by all participants of the screening and confirmation study groups and, in the case of minors, by their parents. This study was further approved by the Ethics committees of the respective universities and institutions and was performed in accordance with the Declaration of Helsinki.

Furthermore, the European, non-Finnish study group extracted from the Genome Aggregation Database (gnomAD; v2.1.1; GRCh37; Karczewski et al., 2020; https://gnomad.broadinstitute.org/; see Supplementary Table S3) was used as an additional control population. Data from either sex were considered as controls for the study group with severe obesity, while data from females served as controls for the screening group of exclusively female patients with AN. Phenotypic data regarding the gnomAD controls and genotype data pertaining to specific variants were unavailable. Additionally, the number of samples varied for each variant.

For gene expression and correlation analyses (based on RNA-sequencing data), paired omental visceral and abdominal subcutaneous adipose tissues were collected from 1,479 individuals of the Leipzig Obesity BioBank (LOBB). Individuals were either non-obese (n = 31; 51.6% female; age: 55.8 ± 13.4 years old; BMI: 25.7 ± 2.7 kg/m²) or had obesity (n = 1,448; 71.2% female; age: 46.9 ± 11.7 years old; BMI: 49.2 ± 8.3 kg/m²). Tissue samples were collected during elective laparoscopic abdominal surgery as described before (Langhardt et al., 2018). Laboratory measurements of body composition and metabolic parameters were obtained as detailed previously (Kloting et al., 2010). The study was approved by the Ethics Committee of the University of Leipzig (approval no: 159-12–21052012) and performed in accordance with the Declaration of Helsinki. All individuals gave written informed consent before participating in the study.

The coding region of the genes CKB (chr14: 103,519,667–103,552,833; GRCh38; ENST00000348956.7), CKMT1B (chr15: 43,593,054–43,604,901; GRCh38; ENST00000441322.6), and GATM (chr15: 45,361,124–45,402,327; GRCh38; ENST00000396659.8) were Sanger sequenced in 192 children and adolescents with severe obesity and 192 healthy-lean controls (see Table 2). The CDS of CKB was further analyzed in 192 female patients with AN (see Table 2). The respective genomic fragments were amplified by performing polymerase chain reactions (PCR; Veriti96-well Thermal Cycler, Applied Biosystems, Foster City, CA, United States) with coding region-specific primers (see Supplementary Table S1). For the amplification of the CDS of CKMT1B, a nested PCR approach was conducted (see Supplementary Table S1). Confirmation of the desired fragments was achieved by 2.5% agarose gel electrophoresis. Subsequent Sanger sequencing was performed by MicroSynth Seqlab GmbH in Göttingen, Germany. Sequence analysis and genotype assignment were conducted independently by at least two experienced scientists using the SeqMan Pro software (version: 11.0.0; DNAstar, Inc., Madison, WI, United States). Discrepancies were solved by either reaching consensus or by re-sequencing.

The CKB- and CKMT1B-located non-synonymous variants identified in patients with severe obesity were additionally genotyped in a larger independent confirmation study group, either by a TaqMan assay (ThermoFisher Scientific, Inc., Waltham, MA, United States) or a restricted fragment length polymorphism (RFLP) analysis (see Figure 1). Therefore, 781 trios with severe obesity as well as 320 children and adolescents with severe obesity and 253 healthy-lean controls from a previous GWAS (Hinney et al., 2007; Scherag et al., 2010) were analyzed (see Table 2). Custom TaqMan assays were designed for the variants rs758572075 (assay ID: ANXGZUF) and rs1230355611 (assay ID: ANYMVED). For the SNP rs146047573, a pre-designed assay was ordered (assay ID: C_162127713_10). All TaqMan assays were performed on the StepOnePlus Real-Time PCR system by Applied Biosystems (Foster City, CA, United States). Genotyping of SNP rs149544188 was performed by subjecting the samples to the restriction enzyme BsaHI (R0556S; New England Biolabs GmbH; Ipswich, MA, United States). Genotypes were assigned by at least two experienced scientists, while discrepancies were solved by either reaching consensus or by repeating the experiments.

To exclude putative genotyping errors, the Hardy-Weinberg-Equilibrium (HWE) was checked. For two variants, a significant deviation of the HWE was detected. Consequently, the respective genotyping was re-checked and repeated. Once genotyping errors were excluded, the corresponding variants remained in the analyses. Furthermore, the transmission disequilibrium test (TDT) was conducted in the 781 trios with severe obesity for one variant analyzed in the confirmation study group (see Table 2; Spielman et al., 1993). Next, putative disease associations with either obesity or AN were assessed by a two-sided chi-square test or a two-sided Fisher’s exact test using IBM SPSS Statistics software (version 28.0.0.0). Associations were considered nominally associated if p < 0.05.

All identified variants were subsequently analyzed by various in silico tools regarding their pathogenic implications, conservation (MutationTaster2021; Steinhaus et al., 2021), and potential to alter splice sites (ESEfinder, Cartegni et al., 2003; SpliceMan; Lim and Fairbrother, 2012). The impact of non-synonymous variants on the protein structure was examined with the tool PANTHER-PSEP (position-specific evolutionary preservation; Tang and Thomas, 2016) and the Cologne University Protein Stability Analysis Tool (CUPSAT, Parthiban et al., 2006). The following protein data bank (PDB) reference protein structures were utilized: CKB - 3DRE (Bong et al., 2008), CKMT1B- 1QK1 (Eder et al., 2000), and GATM - 2JDW (Humm et al., 1997).

RNA-sequencing (RNA-seq; rRNA-depleted, single-end) data from the LOBB cohort were generated based on a SMARTseq protocol (Picelli et al., 2014). In brief, RNA enrichment and reverse transcription were conducted with oligo (dT) and TSO primers. cDNA amplification was carried out with ISPCR primers, and cDNA was processed with Tn5 using the Nextera DNA Flex kit. RNA libraries were sequenced on a Novaseq 6,000 instrument at the Functional Genomics Center Zurich (FGCZ), Switzerland.

The raw reads were adapter trimmed and checked for sufficient quality by applying Fastp (v0.20.0; Chen et al., 2018), permitting a minimum read length of 18 nts and a quality cut-off of 20. Aligned to the human reference genome (GRCh38.p13), gene level expression quantification (gene model definition from GENCODE release 32) was carried out using Kallisto (v0.46; Bray et al., 2016). Samples with more than 20 million read counts were downsampled to 20 million read counts using the R package ezRun (v3.14.1; https://github.com/uzh/ezRun). Data were homoscedastically normalized with respect to library size using the variance-stabilizing transformation implemented by DESeq2 (v1.32.0; Love et al., 2014). To effectively neutralize the effects of in vitro RNA degradation, normalized counts were calibrated with transcript integrity numbers (TINs; Wang et al., 2016). TINs were estimated using the R package RSeQC (v4.0.0; Wang et al., 2012). Finally, the normalized data were adjusted for the gender and age batches.

Correlation analyses were performed by determining Pearson’s correlation coefficient and a confidence interval of 0.95. Between-subject comparisons for the genes were conducted with a non-parametric one-way Kruskal–Wallis ANOVA and pairwise Dunn’s test utilizing the R package ggstatsplot (v0.9.1; Patil, 2021). p-values were corrected for multiple inference using the Hommel’s method (Hommel, 1988). Adjusted two-sided p-values (p_adj.) < 0.05 were considered significant. Analyses were performed under R version 4.2.1.

First, we analyzed whether GWAS-identified SNPs for BMI in CKB, CKMT1B, and GATM show sexually dimorphic effects. In the investigated GWAS (Pulit et al., 2019), a total of 24 SNPs were located in the genomic region of CKB, while 35 and 256 were encompassed in the regions of CKMT1B and GATM, respectively (see Supplementary Table S2). Notably, only CKB contained three SNPs (rs1803283, rs8017780, and rs1136165) that were associated with BMI (GWAS-derived p < 5*10^–8) when both sexes were analyzed concomitantly. Two of these variants (rs1803283 and rs1136165) showed a significant BMI association in females (see Figure 2A), with their minor alleles predisposing to a lower BMI. No BMI-associated SNPs were detected in the region of CKMT1B and GATM (see Figure 2B; Supplementary Table S2). Calculating the z-scores with harmonized data of BMI-associated SNPs in CKB revealed that one variant (rs1136165; 4.17% of all CKB-located SNPs included in the GWAS) showed a more than 3-fold standard deviation aberration of effect sizes between females and males (absolute z-score = |3.196|; p = 0.0014; see Figure 3; Supplementary Table S2). As the tested allele at rs1136165 was solely associated with BMI in females (GWAS-derived p = 7.51*10^–14 for tested allele T; GWAS-derived beta = 0.02 for tested allele T; Pulit et al., 2019) and not in males (GWAS-derived p = 0.013 for tested allele T; GWAS-derived beta = 0.008 for tested allele T; Pulit et al., 2019), this sexual dimorphism indicates a stronger BMI-altering effect in females than males. In the 500 kb region surrounding CKB, 225 SNPs (2.05% of all SNPs in this ±500 kb region) showed a significant sex-dimorphic effect (z-score > |3|). The highest z-score in this adjacent region was detected for rs72712857 (z-score = |4.27|; p = 0.00002; see Supplementary Table S2), which is exclusively BMI-associated in females (GWAS derived p = 1.44*10^–9) (Pulit et al., 2019).

Given the lack of BMI associations within CKMT1B and GATM and adjacent regions (±500 kb), we solely calculated the z-score for one SNP that was associated with BMI in the sex combined analysis (rs765787; GWAS-derived p = 1.80*10^–9) located in the 500 kb region surrounding the two genes of interest (CKMT1B and GATM; see Figure 2B; Supplementary Table S2). However, no sex-specific effect was detected (z-score = |0.38|, p = 0.71). Likewise, in the GWAS data for AN (Watson et al., 2019), we did not find SNPs associated with AN within the gene regions of CKB, CKMT1B, and GATM or in adjacent areas (see Figures 2C,D; Supplementary Table S2). Thus, we did not further analyze putative sex-specific effects within this dataset.

Our sex-stratified analyses suggested a sex-dimorphic effect of the variant rs1136165 (p.Arg152=) located within CKB (see Figures 2, 3; Supplementary Table S2). Considering this GWAS hit and the observation that 1) AN and BMI share certain loci (Bulik-Sullivan et al., 2015; Hinney et al., 2017; Zheng et al., 2022), 2) AN mainly affects women (Hübel et al., 2019; Mitchell and Peterson, 2020), and 3) our GWAS-related sex-specific analyses may have been affected by insufficient power due to the limited sample size of the AN GWAS, we performed a mutation screen of CKB’s CDS in 192 female patients with AN (acute or recovered) as well as in 192 children and adolescents with severe obesity and 192 healthy-lean controls (see Table 2).

In total, we observed five variants (rs762206402, rs1136165, rs200890896, rs146047573, and rs1803283) in the CDS of CKB, one of which is non-synonymous (rs146047573, p.Tyr269Cys; see Table 3). All variants fulfilled the HWE.

While the synonymous variant rs1136165 (p.Arg152=) was shown to have a sex-dimorphic effect on BMI (see Figures 2, 3; Supplementary Table S2), similar frequencies were found in our cases with severe obesity or AN as well as in healthy-lean controls (see Table 3). Hence, we were unable to detect associations of the T-allele (based on the forward strand) at rs1136165 with obesity (p = 0.72) or AN (p = 0.63) in our screening group (see Supplementary Table S4 and Supplementary Table S5). Nevertheless, in silico tools predicted a putative splice site-altering effect (see Table 4).

The rare, non-synonymous variant (rs146047573, p. Tyr269Cys) was detected heterozygously in a young girl from our screening group with severe obesity (BMI: 25.53 kg/m²; BMI-SDS: 2.23; 96th BMI percentile; age: 12.54 years) and neither in patients with AN nor in healthy-lean controls (see Table 3). It was predicted to be conserved and might alter splice sites, while the overall pathogenicity was indicated to be benign (see Tables 4 and 5). However, it putatively destabilizes the protein (see Table 5). Consequently, we genotyped rs146047573 in additional independent confirmation samples of 781 families with severe obesity (child or adolescent with severe obesity and both biological parents), 320 children and adolescents with severe obesity, and 253 healthy-lean controls (see Table 2). Thereby, we detected one additional female with severe obesity carrying the mutation heterozygously (BMI: 35.85 kg/m²; BMI-SDS: 2.96; 100th BMI percentile; age: 16.60 years) and no further carriers among the controls (see Table 3). No obesity association for the alternative allele (C) of rs146047573 was detected in our study groups (p = 1; see Supplementary Table S4). However, a nominal association of this variant’s alternative allele with obesity was evident upon comparison of an alternative non-Finnish European control group extracted from gnomAD with the screening group of 192 children or adolescents with severe obesity (p = 0.03) and our larger study group, including our screening and confirmation group (total of 1,293 children or adolescents with severe obesity; p = 0.02; see Table 2; Supplementary Table S4).

The rare synonymous variant rs200890896 (p.Ala175=) was exclusively detected in one heterozygous female with lifetime AN (BMI: 19.53 kg/m²; BMI-SDS: 1.84; age: 45.68 years see Table 3). It was predicted to be benign while still possessing the potential to alter splicing products (see Table 4). Within our study groups, no association with AN could be ascertained (p = 1; see Supplementary Table S5). However, a nominal AN association of the variant’s G allele was found (p = 0.02) when consulting an exclusively female control group extracted from gnomAD (see Supplementary Table S3 and Supplementary Table S5).

An additional rare synonymous variant (rs762206402; p.Pro112=) indicated to exert non-pathogenic effects (see Table 4) was detected in one heterozygous adolescent with severe obesity (BMI: 39.12 kg/m²; BMI-SDS: 3.08; 100th BMI percentile; age: 15.22 years see Table 3) and not in any other group. The infrequent C-allele of rs762206402 was nominally associated with obesity (p = 0.006) when analyzing the alternative gnomAD control group (see Supplementary Tables S3, S4). Additionally, one synonymous variant (rs1803283, p.Glu364 =) was frequently found (MAF ∼40%) in all three investigated screening groups (see Table 3). Non-pathogenic but putatively splice site-altering effects were determined by in silico tools (see Table 4). Previously, a BMI association of rs1803283 was found in females (p = 1.71*10^–12) and both sexes combined in a GWAS (p = 1.48*10^–16; Pulit et al., 2019).

Nine variants in the protein-coding sequence of CKMT1B were detected in the mutation screen of 192 children and adolescents with severe obesity and 192 lean control individuals (see Table 3). Notably, deviations from the HWE were detected for the novel variant (p.Thr79=) in our screening group with severe obesity and rs13234 (p.Ala362b =) in control individuals. Accordingly, we re-checked and repeated our sequencing in the corresponding samples. The results did not differ from the initial analyses. For instance, for the novel variant p.Thr79=, we detected one heterozygous and one homozygous carrier (see Table 3), leading to deviations from the HWE. Nevertheless, as we were able to exclude genotyping errors, it is feasible that the deviation occurred by chance. Therefore, we retained these variants in the subsequent analyses. All other detected variants fulfilled the HWE.

In total, we observed six synonymous (novel p.Thr79=, novel p.Asp120=, novel p.Thr142=, rs9571, rs13234, and rs144820945), one non-sense (rs1230355611) and two non-synonymous variants (rs758572075 and rs149544188) in the CDS of CKMT1B. The non-synonymous variant rs758572075 leading to a non-conservative amino acid exchange (p.Val191Ala) was found in one individual with obesity (BMI: 35.79 kg/m²; BMI-SDS: 2.96; 100th BMI percentile; age: 16.72 years) and two healthy-lean controls, while the other non-synonymous variant, rs149544188 (p.Arg399His) was identified twice among the children and adolescents with severe obesity and once in the lean control group of the screening samples (see Table 3). Both non-synonymous variants are likely benign but potentially alter splice sites (see Table 4). Furthermore, rs758572075 was predicted to have a stabilizing effect on CKMT1B as well as a favorable impact on its torsion (see Table 5). The non-sense variant rs1230355611 (p.Arg184Stop) was exclusively detected in one lean control individual analyzed in our screening study group (see Table 3). The in silico analyses indicated a deleterious effect of this non-sense mutation (see Table 4).

To evaluate a potential relevance for weight regulation of the protein-altering variants, we genotyped the non-sense (rs1230355611) and both non-synonymous (rs758572075 and rs149544188) variants in additional independent confirmation study groups. Here, within the 781 families with severe obesity, 320 independent children and adolescents with severe obesity, and 253 additional controls (see Table 2), we detected the non-synonymous variant rs758572075 (p.Val191Ala) in one additional boy with obesity (BMI: 32.85 kg/m²; BMI-SDS: 2.57; 100th BMI percentile, age: 14.36 years). The boy was heterozygous for this variant, which the mother who had obesity herself transmitted (BMI: 38.03 kg/m²; BMI-SDS: 1.95, age: 42.94 years), while the lean father (BMI: 22.63 kg/m²; BMI-SDS: 1.54, age: 48.28 years) did not harbor the variant. No association of the C-allele of this SNP with obesity was detected within our screening (p = 1) and confirmation study groups (p = 1; both groups combined: p = 0.27) (see Supplementary Table S4).

Moreover, the non-synonymous variant rs149544188 (p.Arg399His) was detected heterozygously in eleven additional patients with obesity and five healthy-lean controls (see Table 3). Of the eleven cases with obesity, nine were derived from the index cases of the 781 trios with severe obesity. In five cases, the mother transmitted the mutant allele, while the paternal allele was inherited in another four. In an additional 16 families, the child or adolescent with obesity was not a carrier of rs149544188, but one parent harbored the variant heterozygously (not shown). Six fathers and ten mothers were mutation carriers without transmitting the respective allele to their offspring. Consequently, the transmission rate for the mutant A-allele was 36%, indicating a transmission disequilibrium with a reduced transmission rate of the infrequent A-allele (p_chi-square = 0.16).

Furthermore, the non-sense variant rs1230355611 (p.Arg184Stop) was observed heterozygously in one girl with severe obesity (BMI: 31.34 kg/m²; BMI-SDS: 2.52; 99.fourth BMI percentile, age: 14.62 years; see Table 3). It was transmitted from the normal-weight mother (BMI: 24.98 kg/m²; BMI-SDS: 0.11; age: 45.12 years) to the offspring with severe obesity. There was no evidence of an obesity-predisposing effect of this variant (screening group: p = 1; confirmation group: p = 1; both groups: p = 0.45; see Supplementary Table S4).

Evidence for further phenotype-related associations was lacking (see Supplementary Table S4). Only when considering an additional control study group extracted from gnomAD (see Supplementary Table S3), the G-allele of the synonymous rs9571 (p.Lys352 =) was nominally associated (p = 0.00001) with obesity (see Supplementary Table S4). This variant was detected with the same frequency in our 192 patients with severe obesity and healthy-lean controls (MAF = 0.01; see Table 3).

Sequencing of the coding region of GATM in 192 children and adolescents with severe obesity and 192 healthy-lean controls revealed one synonymous (rs1145086; p.Leu418=) and four non-synonymous variants (rs1288775, p.Gln110His; rs146057680, p.Asp234Val; rs773358289, p.Pro287Ser; and rs747005297, p.His292Arg; see Table 3). All identified variants fulfilled the HWE and were detected heterozygously. The missense variant rs1288775 and the synonymous SNP rs1145086 were frequently detected in children and adolescents with severe obesity and healthy-lean controls (MAF >20%; see Table 3). Additionally, two rare non-synonymous variants (rs146057680, p.Asp234Val, and rs773358289, p.Pro287Ser) were observed exclusively among controls. One non-synonymous variant, rs747005297 (p.His292Arg), was exclusively identified in a boy with severe obesity from our study group (BMI: 29.55 kg/m²; BMI-SDS: 2.90; 100th BMI percentile; age: 9.92 years). In silico analyses predicted the potential to affect splice sites and exhibit an overall deleterious effect (see Table 4). All non-synonymous SNPs were classified as either protein destabilizing or unfavorable regarding GATM’s torsion (see Table 5). Except for rs1288775 (p.Gln110His), all other non-synonymous variants were predicted to be deleterious (see Table 4). However, we were unable to find evidence for nominal associations of the variants localized in GATM with obesity (see Supplementary Table S4). But, when investigating the additional control cohort of Europeans extracted from gnomAD (see Supplementary Table S3), a nominal obesity association of the alternative C-allele of rs747005297 (p = 0.009) was observed (see Supplementary Table S4).

To complement our results with transcriptional data, we used RNA-seq data from 1,479 individuals of the LOBB. Based on these data, we analyzed correlations between CKB, CKMT1B, and GATM expression levels and specific metabolic and anthropometric parameters in both VAT and SAT of participants with and without obesity.

The unadjusted correlation analyses (p < 0.05) suggested relationships between several anthropometric and metabolic parameters and our genes of interest. In VAT, CKB correlated positively with CKMT1B (n = 1,479; ρ = 0.43; p < 0.001) and GATM (n = 1,479; ρ = 0.12; p < 0.001) as well as with body fat (n = 669; ρ = 0.11; p = 0.05), waist (n = 244; ρ = 0.16; p = 0.01), and non-esterified fatty acids (NEFA; n = 44; ρ = 0.30; p = 0.05; see Supplementary Figure S1). In the SAT, inverse correlations were detected. CKB was negatively correlated with CKMT1B (n = 1,479; ρ = −0.26; p < 0.001) and GATM (n = 1,479; ρ = −0.37; p < 0.001) as well as BMI (n = 1,479; ρ = −0.06; p = 0.01), homeostatic model assessment for insulin resistance (HOMA-IR; n = 416; ρ = −0.11; p = 0.02), NEFA (n = 44; ρ = −0.35; p = 0.02), and C-reactive protein (n = 1,435; ρ = −0.07; p = 0.005; see Supplementary Figure S1). CKB correlated negatively with hemoglobin A1C (hbA1C) in VAT (n = 790; ρ = −0.07; p = 0.04) and SAT (n = 790; ρ = −0.08; p = 0.03).

However, after correction for multiple inferences, only the inter-gene correlations were significant in both adipose tissues (see Figure 4 and Supplementary Figures S2, S3). As already seen in the data not corrected for multiple testing, all genes were positively correlated with CKB in VAT (CKMT1B: n = 1,479; ρ = 0.43; p_adj < 0.001; GATM: n = 1,479; ρ = 0.12; p_adj < 0.001; see Figure 4A and Supplementary Figure S3), while CKB was negatively correlated with CKMT1B (n = 1,479; ρ = - 0.26; p_adj. < 0.001) and GATM (n = 1,479; ρ = −0.37; p_adj. < 0.001) in SAT (see Figure 4B and Supplementary Figure S2). CKMT1B and GATM exhibited positive correlations in both VAT (n = 1,479; ρ = 0.21; p_adj. < 0.001 see Figure 4 and Supplementary Figure S3) and SAT (n = 1,479; ρ = 0.12; p_adj. < 0.001; see Figure 4 and Supplementary Figure S2).

Obesity and body weight regulation are influenced by sex hormones (e.g., Brettle et al., 2022; Vigil et al., 2022). As the mean age of our study groups was around 50 years (non-obese: 55.8 years old; participants with obesity: 46.9 years old) and most women experience their menopause by the age of 49 years (e.g., Davis et al., 2015), we subsequently excluded all participants above the age of 50 years (n = 465) to circumvent a putative menopause-based bias. Despite the removal of putatively post-menopausal women, differences in correlations were marginal (see Supplementary Figure S4). The inter-gene correlations in the visceral adipose tissue remained unchanged (see Supplementary Figure S4) from the analyses with all participants (see Figure 4A). Yet, the positive correlation between CKMT1B and GATM in the subcutaneous adipose tissue was not observed when participants above 50 years of age were excluded (see Supplementary Figure S4).

Considering the inverse correlations of CKB with CKMT1B and GATM in visceral and subcutaneous adipose tissue (see Figure 4), we subsequently performed between-subject comparisons of non-obese participants and individuals with obesity from the LOBB for all three genes studied. Increased gene expression levels in the VAT of individuals with obesity in comparison to the SAT of this study group were found for CKB (p_adj. < 0.001 see Figure 5A), CKMT1B (p_adj. < 0.001; see Figure 5B), and GATM (p_adj. < 0.001; see Figure 5C). Deviations between participants with and without obesity in the VAT were exclusively observed for CKMT1B (p_adj. < 0.001, see Figure 5B).

To again prevent bias due to the menopause, we excluded females with an age above 50 years. Again, solely minor differences were detected (see Supplementary Figure S5). Thus, the VAT of participants with and without obesity showed a higher expression of CKB than the SAT of patients with obesity (see Supplementary Figure S5).

Considering these findings, we first explored the putative sex-stratified effects of SNPs located within CKB, CKMT1B, and GATM based on data from a BMI GWAS (Pulit et al., 2019). This analysis revealed one common SNP (rs1136165, p.Arg152 =) located in the genomic region of CKB with a sex-dimorphic effect. The BMI-increasing effect (beta-value derived from GWAS = 0.02; Pulit et al., 2019) was shown to be stronger in females than in males. In the 500 kb adjacent region of CKB, more than 4% of the GWAS SNPs revealed a BMI association. Given that enhancers may have a distance of more than 1 million base pairs to the gene’s promoter they affect (Nizovtseva et al., 2017), it is feasible to assume that some of these BMI-associated SNPs in the CKB adjacent region might disrupt enhancer sequences pertaining to the gene and eventually alter its transcription. In contrast, we did not detect evidence for SNPs exerting any effects, including sex-specific effects on CKMT1B and GATM. Moreover, by analyzing data from the latest AN GWAS (Watson et al., 2019), we could not determine variants associated with AN in the regions of interest (CKB, CKMT1B, GATM, and adjacent regions).

The CDSs of CKB, CKMT1B, and GATM were subsequently screened for variants in 192 children and adolescents with severe obesity as well as 192 healthy-lean controls. Furthermore, considering the abovementioned sex-specific GWAS hit, the observation that (1) AN and BMI share certain loci (Bulik-Sullivan et al., 2015; Hinney et al., 2017; Zheng et al., 2022) and 2) AN mainly affects females (Hübel et al., 2019; Mitchell and Peterson, 2020), we additionally screened the coding regions of CKB as the most likely target regarding creatine metabolism in 192 female patients with AN.

We identified five genetic variants in the coding regions of CKB and GATM each, while the CDS of CKMT1B revealed nine variants. The presence of infrequent non-synonymous variants in CKB and CKMT1B (rs146047573, rs758572075, rs1230355611, and rs149544188) was confirmed in independent study groups comprising 781 families with severe obesity, 320 children and adolescents with severe obesity, and 253 healthy lean controls (see Table 2). No evidence for associations of these variants with either obesity or AN was found (see Supplementary Tables S4, S5). Moreover, three novel variants within CKMT1B were detected. As these variants have not previously been described in any other population or regarding other phenotypes, it is not unlikely that, even though no association was detected, these variants might have a putative relevance for obesity. Furthermore, we detected that the variant rs149544188 located in CKMT1B was nominally less frequently transmitted by a heterozygous parent to the child with obesity than expected (transmission rate of 36%, p_Chi-square = 0.16).

In our study groups, we were unable to determine associations of the detected variants with obesity or AN. Presumably, the comparatively small sample size (see Table 2) as well as the low minor allele frequencies (see Table 3) of the detected variants may account for the lack of associations. Thus, we used larger datasets from gnomAD’s non-Finnish European population as an alternative control group (see Supplementary Tables S3–S5). With this, the infrequent alleles of the variants rs76220640 (p.Pro112 =) and rs146047573 (p.Tyr269Cys) located in CKB, the CKMT1B variant rs9571's allele (p.Lys352 =) as well as the allele of the GATM-located variant rs747005297 (p.His292Arg) were nominally associated with obesity. The infrequent allele at rs200890896 (p.Ala175 =) exhibited a nominal association with AN.

Regarding the interpretation of findings in comparison to the gnomAD control group, it is important to note that no phenotypic data were available regarding the data in this database. The gnomAD database mainly includes genotype information from case‒control studies regarding conditions with an onset in adulthood, such as cardiovascular diseases (Karczewski et al., 2020). Since these diseases are comorbidities of an increased BMI and obesity (Williams et al., 2015; Blüher, 2019), our statistical results using this control cohort might be biased by the recruitment of participants with an elevated BMI. However, as our study group consisted of children and adolescents with severe obesity at the extremes of the BMI distribution in the general population, we can assume that the gnomAD control group included only a few individuals with a similar degree of obesity. Finally, we cannot exclude confounding with respect to ethnic background upon the comparison of our German study groups with the gnomAD control group.

Most of the identified non-synonymous variants in the genes of interest were predicted to impact the protein structure. Previously, alterations in the biochemical and biophysical properties of CKB caused by distinct genetic variations leading to amino acid substitutions have been analyzed (Li et al., 2012). Eight non-synonymous variants in CKB were studied. The results demonstrated that half of the variants’ amino acid substitutions affected the catalytic efficiency. For instance, the resultant amino acid exchange of rs13558 (p.Lys267GLu) increases the enzymatic activity by ∼ 30% compared to CKB containing the wild-type lysine (Li et al., 2012). Further amino acid substitution due to other SNP alleles conversely impaired the catalytic efficiency. In mice, some mutated alleles even implied reduced thermal stability below the physiological temperature of 37 °C (Li et al., 2012). Li et al. (2012) suspected that this impairment is related to thermoregulation in the double-KO mouse model (Streijger et al., 2005; Streijger et al., 2009; Li et al., 2012). Our in silico analyses predominantly pointed to benign effects, while the majority of variants exhibited the potential to alter splicing products. Only four non-synonymous variants were predicted to be deleterious overall (see Table 4). One of these variants (rs747005297) located in GATM, which was nominally associated with obesity (based on the alternative control group from gnomAD), was indicated to be protein stabilizing but unfavorable for GATM’s torsion (see Table 5). However, in silico analyses can merely indicate putative consequences of genetic variants and are no substitute for subsequent in vitro studies, which might reveal underlying pathological mechanisms. Particularly considering findings regarding the double-KO mouse model discussed above (Streijger et al., 2005), it is conceivable that polymorphisms in the genes studied here might have an influence on body weight in humans as well. Additionally, as mouse models already suggest compensatory processes between creatine kinases (Streijger et al., 2005), it is possible that mutations in either kinases or even all three genes would have a more severe impact on BMI.

By analyzing RNA-seq data of VAT and SAT of 1,479 individuals (1,448 with obesity and 31 non-obese) of the Leipzig Obesity Biobank, we detected that all three genes of interest (CKB, CKMT1B and GATM) showed positive correlations with each other in visceral fat. However, a negative correlation of CKB with CKMT1B and GATM was observed in SAT. Further correlations with anthropometric or metabolic parameters, such as BMI or body fat, were not significant after correction for multiple comparisons. For CKB, CKMT1B and GATM, subsequent analyses revealed gene expression differences between VAT and SAT in 1,448 adipose tissue donors with obesity. We found consistently higher gene expression levels of the three genes in VAT than in SAT. Additionally, VAT CKMT1B expression was higher in patients with obesity in comparison to individuals with a BMI lower than 30 kg/m².

As expected for genes belonging to the creatine metabolism pathway (Kazak et al., 2015; Kazak et al., 2017), we detected positive correlations between the genes in VAT. Yet, negative correlations among the three genes were found in SAT. Moreover, we found positive correlations between uncorrected gene expression of CKB, CKMT1B, and GATM with body fat in VAT and a negative correlation with BMI in SAT. After correcting for multiple testing, however, these correlations were not significant which indicate small effects that require even larger sample sizes.

CKB was negatively correlated with the other two genes in SAT, while positive correlations were observed in VAT. Differences between the two fat depots are known (e.g., Korsic et al., 2012; Spoto et al., 2014; Ronquillo et al., 2019), and CKB may reflect fat depot specific expression signatures. Expression differences of creatine metabolism pathway genes may therefore reflect functional differences between visceral and subcutaneous fat depots that have been previously reported (Gesta et al., 2006; Spoto et al., 2014; Mittal, 2019; Norreen-Thorsen et al., 2022). Furthermore, global transcriptomic differences in VAT and SAT of individuals with obesity and healthy controls revealed differential expression profiles in a recent study (Ronquillo et al., 2019).

Further, after the removal of 465 female participants older than 50 years, deviations in correlation and gene expression were minimal in comparison to the analysis including all individuals. However, we could depict differences between the inter-gene correlation of CKMT1B and GATM in the subcutaneous adipose tissue when excluding participants above the menopausal age. Based on our data, we can assume that CKB, CKMT1B and GATM are influenced only slightly by sex hormones. Yet, one proteome study in BAT of rats from 2011 revealed 48 proteins exhibiting a sex-specific regulation. Numerous of these proteins affect thermogenesis, such as the creatine kinase, and are more highly expressed in female than male rats (Choi et al., 2011). In order to gain unambiguous evidence whether sex hormones impact the expression of those genes, further and, above all, more targeted analyses are necessary. However, the results of our gene correlation analyses need to be interpreted with caution as the associations were of a small to moderate size (Schober et al., 2018). Nevertheless, we detected adipose tissue-specific expression and correlations of our genes of interest indicating that the alternative creatine-dependent mechanism regulating thermogenesis might be more active in VAT.