Epigenetic Regulation of Processes Related to High Level of Fibroblast Growth Factor 21 in Obese Subjects

We hypothesised that epigenetics may play an important role in mediating fibroblast growth factor 21 (FGF21) resistance in obesity. We aimed to evaluate DNA methylation changes and miRNA pattern in obese subjects associated with high serum FGF21 levels. The study included 136 participants with BMI 27–45 kg/m2. Fasting FGF21, glucose, insulin, GIP, lipids, adipokines, miokines and cytokines were measured and compared in high serum FGF21 (n = 68) group to low FGF21 (n = 68) group. Human DNA Methylation Microarrays were analysed in leukocytes from each group (n = 16). Expression of miRNAs was evaluated using quantitative PCR-TLDA. The study identified differentially methylated genes in pathways related to glucose transport, insulin secretion and signalling, lipid transport and cellular metabolism, response to nutrient levels, thermogenesis, browning of adipose tissue and bone mineralisation. Additionally, it detected transcription factor genes regulating FGF21 and fibroblast growth factor receptor and vascular endothelial growth factor receptor pathways regulation. Increased expression of hsa-miR-875-5p and decreased expression of hsa-miR-133a-3p, hsa-miR-185-5p and hsa-miR-200c-3p were found in the group with high serum FGF21. These changes were associated with high FGF21, VEGF and low adiponectin serum levels. Our results point to a significant role of the epigenetic regulation of genes involved in metabolic pathways related to FGF21 action.


Introduction
Fibroblast growth factor 21 (FGF21) acts as an endocrine factor that regulates glucose and lipids homeostasis in humans and animals although its function is still debated due to its different sites of production and actions [1][2][3]. FGF21 is primarily expressed in the liver but also in the skeletal muscles, adipose tissue, gut and pancreas [4]. FGF21 is regulated by nutritional status (fasting state, high carbohydrate and high protein ingestion) as well as physical activity [3,[5][6][7]. In humans, FGF21 expression in the liver was also found to be up-regulated by glucagon and by glucocorticoids, and repressed by insulin [3]. Circulating FGF21 concentrations exhibit a characteristic diurnal rhythm in humans [8]. Peroxisome proliferator-activated receptor α (PPARα) induces FGF21 transcription in the liver in the fasting state in humans and mice via binding to its promoter [3,9]. Nevertheless, other nuclear receptors such as retinoic acid receptor-related orphan receptor α (RORα), carbohydrate response element-binding protein (ChREBP), farnesoid X receptor (FXR) and liver X receptor (LXR) are also involved in the regulation of hepatic FGF21 expression [10][11][12]. Recently, the stimulatory effect of the fasting state (through KLF15 and BTG2) on hepatic FGF21 expression have been demonstrated in mice [13]. In contrast to the liver, FGF21 expression in brown and white adipose tissues (WAT) is induced mainly by cold exposure and adrenergic signalling as well as PPARγ, which directly mediates the expression of FGF21 in WAT [3,9]. In muscles, FGF21 expression is controlled by ATF4 [14].
FGF21 signalling to target cells is mediated by binding to and activation of four receptors: FGFR1, FGFR2, FGFR3 and FGFR4 in the presence of co-receptor β-klotho [15]. FGFR receptors are ubiquitously expressed, whereas β-klotho, an essential co-receptor, is expressed only in the liver, adipose tissue, pancreas and brain, suggesting that FGF21 functions selectively in these tissues [4]. Although therapeutic administrating of FGF21 prevents diet-induced obesity and related insulin resistance in mice and humans. Paradoxically metabolic disorders like obesity, diabetes are characterised by elevated serum FGF21. This highlights the important role of receptors. A recent study has shown the reduced expression of FGFR1c, FGFR3c and FGFR4 receptors in skeletal muscles and increased expression of FGFR1c in WAT of individuals with type 2 diabetes [16]. Signal activation via FGFR1 leads to weight loss, brown fat thermogenesis, insulin sensitisation and hepatic lipids reduction [17]. FGF21 is widely known to exert multiple metabolic effects: increases insulin sensitivity, fatty acid oxidation and ketogenesis but inhibits lipogenesis in liver; increases fatty acid uptake, lipolysis, browning, thermogenesis and mitochondrial activity in white adipocytes and enhances glucose uptake and thermogenesis in brown adipocytes [3,12].
Serum level of circulating FGF21 increases with age as well as in obesity and insulin resistance state [18][19][20]. High serum FGF21 level is associated with MAFLD (metabolic associated fatty liver disease) and the higher risk of developing metabolic syndrome and type 2 diabetes [21,22]. Elevated FGF21 levels are observed in patients with atherosclerosis, pancreatitis and reduced bone strength [3]. Moreover, recent findings have indicated that high versus low serum FGF21 levels correlate with mortality in the elderly [12]. FGF21 is a part of a complicated inter-organ cross talk, so other factors (e.g., cytokines, metabolites) may modulate FGF21 action in target tissues and consequently elicit a beneficial (healthy) or detrimental (unhealthy) effect.
It is hypothesised that proper action of FGF21 prevents the development of dietinduced obesity [21]. According to a recent hypothesis, FGF21 is a hormonal mediator of the human "thrifty" metabolic phenotype [23]. Efficient regulation, with large changes in FGF21 secretion in response to the supply of nutrients or nutrient deprivation, is one of the main factors determining the metabolic phenotype less susceptible to the development of obesity and its complications [24]. As FGF21 also plays a role in brown adipose tissue activation, it seems, to some extent, to affect the individual propensity to weight gain. However, in obesity the state of FGF21 resistance is widely observed, which is suggested as compensatory up-regulation response to chronically elevated serum levels of FGF21 [19]. Mechanisms mediating the interaction between environmental factors and the genome, such as epigenetics, may be of particular importance in the pathogenesis of FGF21-resistance. In recent study on mice it was demonstrated that Fgf21 gene methylation status can be modulated in a PPARα-dependent manner [25]. Therefore, expanding knowledge of the epigenetic basis associated with blunted FGF21 response to diet, is important for identifying individuals at risk and preventing the metabolic complications of obesity.
The aim of this study was to find differentially methylated genes and profile of miRNAs related to high FGF21 levels in obese subjects. For this purpose, a genome wide DNA methylation and miRNA expression patterns were determined in leukocytes. We chose blood samples, as easily accessible to collect, for screening of epigenetic biomarkers associated with high FGF21 levels. As a result, a blood epigenetic profile characteristic for high serum FGF21 in obesity was established.

Patients
The study group (from the BIOCLAIMS cohort) included adult overweight and obese individuals (n = 136; 100 women and 36 men) with BMI ranging from 27 to 45 (overweight: n = 30, obese: n = 106). Exclusion criteria were cardiovascular diseases, diabetes mellitus, kidney or liver failure, endocrine disorders and chronic inflammation as well as smoking or excessive use of alcohol, pregnancy or lactation. Subjects did not use any drugs or dietary supplements and led inactive lifestyle.
The study was conducted in accordance with the Code of Ethics of the World Medical Association (Declaration of Helsinki) and with the Good Clinical Practice guidelines. The study protocol was approved by the Bioethics Committee of the Jagiellonian University Patients were divided into 2 groups based on the fasting serum FGF21 levels (FGF21 < 213 pg/mL and ≥ 213 pg/mL). The cut-off point was the median serum FGF21 in the cohort.

Anthropometry Measurements
All subjects had measurements of body weight, height, waist and hip circumferences and blood pressure. Body fat percentage was estimated by bioelectrical impedance method using the Segmental Body Composition 116 Analyzer TANITA BC 418 MA (Tanita Corporation, Tokyo, Japan).

Sample Collection and Biochemical Tests
Fasting venous blood samples were always between 8 to 9 am collected and then centrifuged (1000 × g for 10 min at 4 • C) within 30 min of collection for serum and plasma separation. For biochemical tests (except free fatty acids (NEFAs), serum and plasma samples were immediately frozen and stored at −80 • C until assayed. For analysis of DNA methylation, fasting venous blood was collected into K3-EDTA tubes. At the same time, a second sample of venous blood was collected into PAX gene Blood RNA Tubes (Becton Dickinson, Bedford, MA, USA) for miRNAs expression analysis. All samples were frozen and stored at −80 • C until analysed.
Serum was used for biochemical tests: glucose, insulin, non-esterified fatty acids (NEFA), total cholesterol, HDL-cholesterol and triglycerides (TG) concentrations. Glucose, total cholesterol, HDL cholesterol, and TGs were assayed using enzymatic colorimetric methods on the MaxMat analyzer (MaxMat S.A., Montpeliere, France). LDL cholesterol concentration was calculated using the Friedewald formula. Insulin was measured by an immunoradiometric method (INS IRMA Diasource ImmunoAssays, Neuve Belgium). Homeostasis model of assessment (HOMA-IR) was calculated as a ratio according to Matsuda M and DeFronzo RA [26]. Plasma GIP was measured by ELISA (Human GIP (Total) 96-Well Plate Assay (EMD Millipore, St Charles, MO, USA)). Total plasma fatty acids content and composition of saturated, monounsaturated and polyunsaturated fatty acids were measured by gas-liquid chromatography and flame-ionization detector after direct in situ transesterification, according to Glaser et al. [27]. NEFAs concentrations were measured immediately in fresh plasma by an enzymatic quantitative colorimetric method (Roche Diagnostics GmbH, Mannheim, Germany). Activities of ALT and GGT were measured by the routine method using an automated analyser (Hitachi cobas c 701/702, Roche Diagnostics GmbH, Mannheim, Germany).

Genome-wide DNA Methylation Analysis
DNA methylation microarrays were analyzed in 16 samples selected randomly (matched to sex and BMI) from the two studied groups (8 representative for each group-High FGF21 and Low FGF21). The subgroups reflected the biochemical characteristics of the group they represented. Genomic DNA was isolated from peripheral blood leukocytes using the High Pure PCR Template Preparation Kit (Roche Diagnostics, Mannheim, Germany). The assessment of DNA quantity and quality was performed by spectrophotometry using the NanoDrop ND1000 UV-VIS spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). First, sonication of 5 µg of DNA was performed. Sonication efficiency was assessed by electrophoresis on a 2.0% agarose gel. The 4/5 of sonicated DNA was taken for immunoprecipitation, the fifth part was used as a reference input fraction. The immunoprecipitation of DNA containing 5methylcytosines (5-mC) were performed using monoclonal antibodies against 5-methylcytidine (Monoclonal Antibody to 5-Methyl Cytosine / 5-MeC Purified from Acris Antibodies, Inc, USA). Immunoprecipitated and reference DNA aliquots were labelled with fluorescent dyes Cyanine-3 and Cyanine-5, respectively. Competitive hybridisation of input material and methylated DNA enriched fraction was performed to oligonucleotide microarrays-Human DNA Methylation Microarray: High-definition 244K arrays contained 27,627 probes for annotated human CpG islands and 5081 undermethylated regions (G4495A, Design ID: 023795) from Agilent Technologies (Santa Clara, CA, USA). The selected arrays includes 3 types of probes: probes for known CpG islands, probes spanning 85 nucleotides proximal to each CpG island and probes for undermethylated regions representing 237,227 unique biological probes. Microarrays were hybridized for 40 h at 65 • C. Samples amplification and labelling, hybridisation, slides washing and scanning procedures and image extraction using Agilent Features Extraction software v 10.10.1.1 were performed according to the manufacturer's instructions.

Determination of miRNAs Expression by Real-Time PCR -TLDA (TaqMan Low Density Arrays)
For relative quantification of miRNA, 16 blood samples intended for epigenetic analysis were drawn from both group (8 representative for each group-High FGF21 and Low FGF21).
MiRNA was extracted using the PAXgene 96 Blood RNA Kit and GeneMATRIX Universal RNA/miRNA Purification Kit (EUR X , Gdansk, Poland). Finally, after isolation we obtained miRNAs from peripheral blood leukocytes. MiRNA quality was assessed in an Agilent Bioanalyzer 2100 using the RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA) and quantified by spectrophotometry using the NanoDrop ND-1000.

Real-time PCR Analysis
Total RNA was reverse transcribed using a reverse transcription kit (High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Carlsbad, CA, USA)) with random primers (1 µg). Quantitave real-time polymerase chain reaction (qPCR) was performed with the QuantiTect SYBR Green PCR kit. The specific TaqMan Gene Expression Assays with a pair of specific PCR primers and a TaqMan probe with a FAM (Applied Biosystem, Carlsbad, CA, USA) was used. Amplification was performed using the continuous fluorescence detection system 7900 HT Fast Real Time PCR system (Applied Biosystem, Carlsbad, CA, USA). The expression ratio of target mRNA was normalised to the level of 18 s RNA and compared to control group. Data was analysed using the ∆∆CT method. The statistically significant results were recognized, with a value of p < 0.05.

Statistical Analyses
For biochemical and anthropometric parameters, data distribution was analysed by Shapiro-Wilk test. Differences between groups were calculated by unpaired t-test for normally distributed data and Mann-Whitney U-test for non-normally distributed data (BMI, adipose tissue mass, WHR, glucose, HOMA-IR, VEGF, GIP, leptin, adiponectin, FGF21, myostatin, IL-6 and hs-CRP). The chi-squared test was used for nominal variables. Normally distributed data are presented as mean ± standard deviation (SD), otherwise as median (Q2) and interquartile range (Q1; Q3). All analyses were performed with the Statistica 13 software (StatSoft Polska, Krakow, Poland). The p-value < 0.05 was considered statistically significant.
Microarray DNA methylation data analysis was performed using the Feature extraction software version 10.10.1.1 (Agilent Technologies, Santa Clara, CA, USA), the BRB-ArrayTools software version 4.6 (National Institutes of Health, Bethesda, MD, USA) and the R software version 3.6.1 (Development Core Team, Vienna, Austria and Bioconductor). Feature extraction software was used to assess background subtracted intensity values for the two fluorescence dyes on each individual array feature and calculated as the ratio (Cy3/Cy5). Loess normalisation method for dual channel raw hybridisation signals were computed in BRB-ArrayTools. M values (log intensity ratio M = log 2 (Cy3/Cy5)) were used for downstream statistical analyses. A linear regression model for comparing two groups with high and low FG21 was adjusted for potential confounders including blood cell-type proportions and age. To account for potential differences in the proportions of blood cells, we estimated the proportion of lymphocytes, monocytes, eosinophils, basophils and neutrophils to adjust raw data as previously described by Houseman et al. in R [28]. We used the functions lmFit and eBayes in "limma" package to build linear models and compute moderated t-statistics (p-value was calculated). Then Benjamini-Hochberg multiple test correction FDR was performed and shown as q-value. Regarding methylation sites, probes not including CpG or not located in CpG islands, probes containing CNVs in the CpG site and probes located on X and Y chromosomes were removed from the analyses. Fold of change (FC) in relation to obese with high FGF21 versus obese with low FGF21 was used as the methylation level. Loci corresponding to q-value < 0.01 and FC in the range of > 1.4 or < −1.4 were classified as differentially methylated. Methylated regions had log ratios significantly above zero while less methylated regions had log ratios significantly below zero. Highly methylated regions are with FC ≥ 2.0 and moderately methylated with FC between 1.4 to 2.0.
Pathway analyses: We used the BiNGO plugin in Cytoscape software (version 3.7.2) to assess the involvement of selected genes in biological processes and molecular function [29]. The parameters as overrepresentation after correction (using statistical test as hypergeometric test with multiple testing correction as Benjamini-Hochberg False Discovery Rate (FDR) correction) were selected. Results with corrected p-value < 0.05 are presented in this manuscript.

Biochemical and Anthropometric Characteristics of the Study Groups
The two investigated groups were comparable in terms of age, sex, BMI, adipose tissue mass, blood pressure, fasting glucose and cholesterol, regardless of serum FGF21 levels. Nevertheless, the group with high circulating FGF21 was characterised by higher WHR, elevated fasting insulin and HOMA-IR, NEFAs as well as triglycerides levels ( Table 1). This group included 42 (62%) individuals who met the criteria for abdominal obesity, compared to the low FGF21 group, n = 27 (40%). Plasma GIP, an incretin hormone important for regulating carbohydrate as well as lipid metabolism, was also found to be elevated in the high FGF21 group. Liver damage markers (ALT, GGT) also differed between groups, indicating patients with high FGF21 are more prone to liver injury. Obese patients with high FGF21 had elevated VEGF, MCP1 and decreased adiponectin levels whereas other adipokines-leptin, resistin, visfatin, myokines-irisin, myostatin, and IL-6 levels did not differ between these groups (Table 1). Other inflammatory markers studied-CRP, sVCAM, sPECAM, sE-selectin-showed comparable levels regardless of circulating FGF21 ( Table 1). The patients in the study group presented metabolic profile related to elevated FGF21. Characteristics of subjects selected for epigenetic study are presented in Table S1 (in supplementary files).

Results of DNA Methylation Analysis
The analysis of the data obtained from DNA methylation arrays revealed 11 198 CpG probes differentially methylated (p < 0.01, representing approximately 4.6% of probes located on the array), of which 5425 were significantly hypomethylated and 5774 were significantly hypermethylated in the high FGF21 group compared to the low FGF21 group. We selected 375 CpG probes for further analysis of biological significance in pathways related to potential FGF21 action. We focused on those genes that were engaged in metabolism of lipids, glucose and adipokines and other pathways known or potentially related to FGF21 activity. Detailed results of each probe described in the manuscript are presented in Table  S2 enclosed as a supplementary file.
Analyses in Cytoscape software revealed biological processes affected by differentially methylated genes (results presented in Table 2). The most interesting finding is the involvement of these genes in the response to glucose and insulin and other nutrient levels. Subsequently, the differential methylation status of lipid transport and fatty acid metabolism, fat cell differentiation (including white and brown fat cell differentiation), generation of precursor metabolites and energy, response to hypoxia and regulation of ossification, were demonstrated (Table 2). Importantly, in the high FGF21 group, changes in methylation status of genes involved in browning program of adipose tissue: BMP4 (hypomethylated in promoter), FGF9 (hypermethylated in promoter) as well as PRDM16 (hypermethylated in promoter) and CPT1C (hypomethylated in promoter) were detected. Genes: MRAP (shown as hypomethylated in the high FGF21 group) and TBX1 (hypermethylated in the high FGF21 group) are the hallmark proteins, for fatty tissue (highly enriched in differentiated adipocytes) and beige adipocytes, respectively.
Interestingly, five genes coding for transcriptional factors: PPARA, RORA, ATF4, NR3C1 and KLF15, known to regulate FGF21 expression, were differentially methylated ( Table 3). Analysis of DNA methylation also revealed regulation of fibroblast growth factor receptor signalling (hypomethylated genes: FGF3, KLB, FGFR1 and FGFR3) and vascular endothelial growth factor receptor signalling (hypomethylated genes: VEGFA and FLT1) pathways (Table 2). High FGF21 individuals were also characterised by the LIPG gene, as well as ABCA1 and VLDLR genes (all hypermethylated) related to lipoprotein metabolism and transport. We found differential methylation in promoters in genes: IGF1R -receptor for IGF1 and IGFBP1an FGF21-induced pro-osteoclastogenic liver hormone (all with hypermethylated status in promoters in the high FGF21 group) ( Table 2). Interestingly, the promoter of GIPR (receptor for GIP) gene was hypermethylated in this group. We showed at the protein level that serum GIP and VEGF levels were significantly higher in obese patients with high FGF21 than in subjects with low FGF21.   Further analysis of differentially methylated genes in ClueGO plugin in CytoScape software allowed to create graphic representation of biological processes and pathways ( Figure S1). Figure S1 is enclosed as a supplementary file.
DNA microarray results was verified by real-time PCR in 17 samples (8 from high FGF21 group and 9 from low FGF21 group) for selected ATF4 gene. We demonstrated a statistically significant 2.75-fold decrease in expression of the ATF4 gene in leukocytes (p < 0.05). This result confirms that detected DNA hypermethylation in the promoter of ATF4 gene regulates its expression in patients with high serum FGF21 level.

Results of microRNA Expression Level in Leukocytes:
We demonstrated a statistically significant difference in expression of four miRNAs in peripheral blood leukocytes related to high FGF21 serum levels. Expression of the following miRNAs: hsa-miR-133a-3p, hsa-miR-185-5p and hsa-miR-200c-3p was downregulated and hsa-miR-875-5p was up-regulated (results are presented in Table 4). In Table 4, we also selected targets for miRNAs in relation to high FGF21 (among others the FGFR1 gene, receptor for FGF21) as well as related to high serum VEGF (VEGFA and its binding receptors genes) and low adiponectin (the ADIPOQ gene), found to accompany high FGF21 serum levels. Presented in Table 4, target genes for hsa-miR-133a-3p, hsa-miR-185-5p and hsa-miR-200c-3p were validated in previous manuscripts and/ or information based on selected online datasets. For hsa-miR-875-5p, we chose some putative target genes located in 3'UTR from miRWalK and TargetScan databases.

Discussion
In the current study, we have demonstrated for the first time genome-wide DNA methylation profile and miRNAs expression analysis in white blood cells associated with high serum FGF21 levels in human obese non-diabetic individuals. In our cohort, wellcharacterised in terms of many circulating active mediators, such as myokines, adipokines, incretins, cytokines and fatty acids, we observed increased parameters of insulin-and GIP resistance, lipids disorders, elevated VEGF, MCP1 and lower adiponectin serum levels suggesting disturbed nutrients signalling and disregulated inter-organ cross-talk in the high FGF21 group compared to the low FGF21 group.
It is well known that the altered DNA methylation of CpG sites in the genes' promoters, as well as distal regulatory sites, changes gene expression patterns by modifying the interaction of histones, thereby affecting the binding of transcription factors or recruitment of methyl-CpG binding proteins (MBPs) [39,40]. The detected and presented, differentially methylated genes may be responsible for aberrant expression of genes related to high FGF21 level in obese people.
Regarding lipid metabolism, differentially methylated genes: PPARA, LPIN1, CPT1A, CPT1B, CPT1C, ANGPTL4 and FADS2 are key regulators of liver, adipose tissue, muscle and total lipid homeostatis. PPARA (nuclear transcription factor) together with LPIN1 (acting as a nuclear transcriptional coactivator for PPARGC1A and PPARA) regulate various lipid metabolic pathways [42,43]. In the liver PPARA induces the expression of genes involved in fatty acid oxidation, binding, activation, elongation and desaturation, synthesis and hydrolysis of triglycerides and lipid droplets, lipoprotein metabolism, gluconeogenesis and bile acid metabolism pathways [44]. It has been demonstrated on mice models that enhanced lipin expression increases lipogenic gene expression in adipose tissue and promoted fat storage in obesity, whereas in skeletal muscle it participated in energy expenditure and fat utilisation [45]. CPT1C gene plays a role in the central control of feeding behaviour and whole-body energy homeostasis [46]. Hypomethylation of leukocytic CPT1C gene was previously thought to be associated with increased risk of developing MAFLD [47].
Among five detected TFs known to regulate FGF21 expression, the majority (except 1 probe in NR3C1 gene) were hypermethylated in patients with high FGF21. PPARA was activated by endogenous fatty acids and other PPARα agonists and, in turn, activated expression of FGF21 and genes involved in fatty acids beta-oxidation and ketogenesis pathways [33,48].
Research on mice showed that PPARα modulated the methylation of liver Fgf21 and represents a form of epigenetic memory that persists from the postnatal period into adulthood that influence the risk of obesity in later life [25]. On the other hand, we did not observe any differences in methylation status of PPARG, regulating FGF21 expression in WAT, between the groups. It has been reported previously that glucocorticoid receptor decreases hepatic PPARα expression and, in consequence, reduces liver FGF21 expression and plasma levels [10] and suppression of nuclear receptor RORA causing a decrease in FGF21 expression [11]. Observed hypermethylation of TFs (KLF15, NR3C1 (glucocorticoid receptor), RORA (all regulating hepatic FGF21 expression) and ATF4 (regulating muscle FGF21 expression)) suggests disturbed regulation of FGF21 in tissues. It was identified that glucose induced Fgf21 mRNA expression through ChREBP activation, as the ChoRE is located in the Fgf21 promoter. The feedback systems of ChREBP regulation are composed of two factors: PPARα and Fibroblast growth factor-21. Negative feedback is achieved by FGF21, which interacts with PPARα to repress their own transcription [49]. Although PPARα is the most important regulator of hepatic FGF21 expression, it is possible that other TFs may be down-regulated by FGF21, which is suggested by our DNA methylation results.
In our study, besides adipomyokine FGF21, two other myokines-myostatin and irisinwere measured in the serum. Since circulating myostatin and irisin did not differ between the groups, we supposed that FGF21 was mainly of hepatic origin and not skeletal muscle. Thus, we postulated that high FGF21 levels in our studied group were associated with a metabolic disorder, not induced by exercise. Receptors mediate FGF21 action in target tissues. FGFR1 and FGFR2 are ubiquitously expressed growth factor receptors that mediate most biological functions of the FGF family [4]. Regarding liver FGF21 receptors, hepatic mRNA expression is the highest for FGFR4 and FGFR2, then for FGFR3 and the lowest for FGFR1 [50][51][52]. We detected differential methylation in FGFR1 and FGFR3 receptor genes and their co-factors: FGFRL1 and KLB genes. This is in line with previous studies where it was demonstrated that hepatic levels for β-Klotho, FGFR1 and FGFR3 transcripts were significantly increased in patients with obesity [53]. Interestingly, patients in the high FGF21 group presented simultaneously high level of serum VEGF protein. We detected differential DNA methylation in pathway: positive regulation of vascular endothelial growth factor receptor signalling (with affected genes: FLT1, VEGFA and FGF9).
As non-coding microRNAs are post-transcriptional regulators of genes translation, we analysed the microRNA profile in patients with obesity and high FGF21 levels. Analysis of miRNAs expression in leukocytes revealed a statistically significant decrease in expression of has-miR-133a-3hashsa-miR-185-5hasnd hsa-miR-200c-3p in the group with high FGF21. Previously, it was demonstrated that miR-133 directly targeted and negatively regulated Prdm16 in mice [54]. In humans, it was shown that circulating miR-133 is high in patients with periprocedural myocardial injury and, as a consequence, decreased fibroblast growth factor (FGFR1) [55]. As FGF21 is known to be activated by stress conditions, low hsa-miR-133a-3p expression in our studied group may be a kind of response to metabolic stress. According to a recent hypothesis, FGF21 is a hormonal mediator of the human "thrifty" metabolic phenotype; thus, it may be one of the mechanisms of metabolism towards browning of adipose tissue. In studies on human cancer cells inverse correlation of miR-185 and miR-200c to modulate VEGFA translation was found [56,57]. This profile of miRNAs is in line with elevated serum VEGF level detected in the group with high FGF21 from our cohort. The present study revealed the increased expression of has-miR-875-5p. There is limited information regardihashsa-miR-875-5p. Potential involvemehasof hsa-miR-875-5p to regulate TNF, LEP and IRS1 genes in gestational diabetes mellitus was demonstrated in a recent study [58]. Previous studies have shown decreased lhasls of hsa-miR-875-5p in the adipose tissue of patients with non-alcoholic steatohepatitis (NASH) [59]. Searching the miRNAs databases, we found potential targets for hsa-miR-875-5p such as: ADIPOQ Simultaneously, we demonstrated statistically significant lower adiponectin protein level in serum in individuals with high FGF21; however, to confirm that hsa-miR-875-5p regulate the ADIPOQ gene are required further miRNA/mRNA interaction validations.
The main limitation of our study was the analysis of genome-wide DNA methylation and miRNAs expression in limited samples representative for both groups. Thus, results from our high-throughput method gives rise for further studies on targeted DNA methylation. We assessed DNA methylation in peripheral blood cells because they are easily accessible and the collection is acceptable by patients. Although epigenetic studies on different tissue samples are more informative, blood samples are generally used in most studies with non-surgery subjects. Nevertheless, previous studies demonstrated that levels of DNA methylation in blood tend to be broadly correlated with levels in other target tissues [60,61]. This suggests that DNA methylation and miRNAs analysis in leukocytes may reflect epigenetic changes in other tissues (liver, muscles, adipose tissue and bones) relevant for the pathogenesis of glucose and lipid disturbances related to high FGF21 in obesity. Another limitation is investigations among participants with metabolically healthy obesity, excluding general obese population with associated diabetes and cardiovascular diseases or other obesity complications.

Conclusions
The major finding of this preliminary study is that individuals with high FGF21 presents specific DNA methylation and miRNAs profiles in blood leukocytes. Altered DNA methylation in promoters of transcription factor genes (PPARA, KLF15, NR3C1, RORA and ATF4) known to regulate FGF21 expression and its binding receptors and co-factors (FGFR1, FGFR3 and FGFRL1 and KLB) has been revealed in relation to the elevated levels of circulating FGF21 in obesity. The mostly regulated processes are insulin secretion and signalling, lipid transport and homeostaais maintenance, thermogenesis and browning of adipose tissue and regulation of ossification and bone mineralization. Additionally, we demonstrated differentially expressed miRNAs known to target and inversely regulate the FGFR1 and VEGFA genes. As DNA methylation is reversible and depends on environmental factors, there is the potential to influence the methylation status of key genes by nutrition and a healthy lifestyle to prevent obesity-related complications [62][63][64][65]. These findings give rise for further studies on detailed dietary factors and targeted DNA methylation editing therapies that may regulate revealed genes.
Supplementary Materials: The following are available online at https://www.mdpi.com/2073-4425/ 12/2/307/s1, Figure S1: Graph of a biological network of identified as differentially methylated genes created by the ClueGO plugin in Cytoscape software, Table S1: Characteristics of subjects selected for epigenetics study, Table S2: Detail results of identified differentially methylated CpGs probes.
Author Contributions: T.P. performed the study, analysed and interpreted the data and drafted the manuscript. A.P. participated in conception and design of the study, analysed the data and revised article critically for important intellectual content. J.G. interpreted the data and revised article critically for important intellectual content. U.R. performed the study and analysed the data.  Data Availability Statement: Data from the BIOCLAIMS study are deposited in Nutritional Phenotype Database http://www.dbnp.org/nutritional-instance/ (accessed on 19 January 2021)). The data of differentially methylated CpGs probes presented in this study are available in Supplementary materials.