Suggestive Evidence for Causal Effect of Leptin Levels on Risk for Anorexia Nervosa: Results of a Mendelian Randomization Study

Genetic correlations suggest a coexisting genetic predisposition to both low leptin levels and risk for anorexia nervosa (AN). To investigate the causality and direction of these associations, we performed bidirectional two-sample Mendelian randomization (MR) analyses using data of the most recent genome-wide association study (GWAS) for AN and both a GWAS and an exome-wide-association-study (EWAS) for leptin levels. Most MR methods with genetic instruments from GWAS showed a causal effect of lower leptin levels on higher risk of AN (e.g. IVW b = −0.923, p = 1.5 × 10−4). Because most patients with AN are female, we additionally performed analyses using leptin GWAS data of females only. Again, there was a significant effect of leptin levels on the risk of AN (e.g. IVW b = −0.826, p = 1.1 × 10−04). MR with genetic instruments from EWAS showed no overall effect of leptin levels on the risk for AN. For the opposite direction, MR revealed no causal effect of AN on leptin levels. If our results are confirmed in extended GWAS data sets, a low endogenous leptin synthesis represents a risk factor for developing AN.

In this study, we focus on leptin levels in anorexia nervosa (AN). AN is a severe psychiatric disorder characterized by underweight and hypoleptinemia in addition to AN specific cognitions and emotions (Hebebrand et al., 1995;Grinspoon et al., 1996;Hebebrand et al., 1997;Hebebrand and Bulik, 2011;American Psychiatric Association, 2013;Misra and Klibanski, 2014). Our study was stimulated by recent case reports suggesting pronounced effects of recombinant human leptin on such cognitions and emotions (Milos et al., 2020;Antel et al., 2021).
Based on recent genetic findings, AN has been defined as a metabo-psychiatric disorder (Watson et al., 2019). This hypothesis is based on genetic correlations to a range of anthropometric and metabolic phenotypes. Thus, Watson et al. (2019) showed inverse genetic correlations between the risk for AN and levels of leptin, fasting insulin, body mass index (BMI), body fat percentage, risk of obesity, insulin resistance, and type 2 diabetes; a positive correlation was observed with levels of high density lipoproteins (HDL). However, not all genetic correlations, including those with leptin levels and type 2 diabetes, remained significant after adjusting for BMI (Watson et al., 2019). The negative genetic correlations between body fat percentage and AN were stronger in women than in men (Hübel et al., 2019).
To evaluate causality and direction of leptin levels for the risk of AN, Mendelian Randomization (MR) is the appropriate approach (Davey Smith and Ebrahim, 2003;Burgess et al., 2015;Zhao et al., 2021). Mendelian randomization uses genetic data and the assumption of a random transmission of genetic information at conception to the next generation with a random distribution of covariates. MR is therefore a kind of a randomized control study (RCT) with quasi-randomization Lawlor et al., 2008). The basic assumption of MR is that there are genetic variants that alter the level of a modifiable environmental exposure or mirror the biological effects of that environmental exposure and that this modifiable environmental exposure alters the risk of a specific disease. Polymorphisms can be markers for such variants with a biological function (Davey . Common genetic polymorphisms identified to be relevant for the environmental exposure of interest e.g. in genome-wide association studies (GWAS) can be used as proxies for the exposure whose measurement is impossible, inaccurate or confounded. In this way it is possible to investigate the effect of this exposure on disease risk excluding the risk of residual confounding. Since the genetic variants cannot be altered by the disease or the environment, reverse causality is ruled out in an MR study, as there is no evidence for selection bias in which participants enter into the genetic studies on which MR is based . Additional advantages of MR over RCTs or longitudinal observational studies are that practical issues (such as high costs and long study duration) or ethical problems can be neglected (König and Greco, 2018). In MR, genetic variants are equivalent to lifetime differences in exposure and indicate the long-term effects on disease (outcome). Accordingly, they generate more realistic estimates of causal effects, because they are free of measurement error or shortterm variations in exposure (Davey Smith and Ebrahim, 2005).
The aim of this study was to investigate whether there is a causal relationship between leptin levels and AN. We performed a two-sample MR which uses two different study samples for the risk factor and the outcome phenotype as data (Pierce and Burgess, 2013). Corresponding reverse MR analyses were conducted to test for opposite causality. We ran analyses based on both leptin levels unadjusted and adjusted for BMI and in data sets based on both sexes combined. However, we a priori considered the analysis based on leptin levels in females as the most valid. The reasons include our intention to attempt to eliminate effects of the well-known association between BMI and leptin levels and to focus on leptin secretion in females. Thus, both a higher percent body fat and increased synthesis of leptin per unit of fat mass potentially explain why females have higher circulating leptin levels than males (Hellström et al., 2000).

Leptin Levels
We performed the MR analyses based on the GWAS by Kilpeläinen et al. (Kilpeläinen et al., 2016), which included 32,161 individuals of European ancestry (stage 1). In stage 2, an additional 19,979 individuals with European ancestry were studied. The measured leptin level (µg/ml) was transformed logarithmically and adjusted for age, sex and study-specific covariates (e.g. genotype-derived principal components). Five loci were genome-wide associated (p < 5 × 10 −8 ) with leptin levels: LEP, SLC32A1, GCKR, CCNL1 and FTO. After adjustment of leptin levels for BMI, FTO did not remain significant. For genome-wide significant single nucleotide polymorphisms (SNP), results (non-adjusted and adjusted for BMI) were published for the total sample and separately for males and females. However, complete summary statistics for leptin levels with/without adjustment for BMI were only available for sexcombined datasets. SNP heritability (h 2 SNP ) was not reported. The most recent study on genetic determination of leptin levels was an exome-based analysis (EWAS) in individuals of diverse ancestries (Yaghootkar et al., 2020). This analysis was based on a total of 57,232 adults from 35 cohorts, of whom 50,321 were of European, 4,387 of African, 2,036 of East Asian, and 488 of Hispanic descent. In order to use the most recent currently available data, we performed MR using the SNPs from this study as an instrumental variable (IV). The studies included in the EWAS study computed residuals for leptin concentrations (in ng/ mL) using linear regression, adjusting for age, genome-wide principal components, and any study specific covariates (e.g. study center). The residuals were calculated with and without adjustments for BMI. The residuals were transformed using rank based inverse normal transformation to follow a distribution with a mean of 0 and a standard deviation of 1. Results were reported for all ancestries combined and separately for Europeans and for both sexes and separately for women and men. In addition, three different models (additive, recessive and dominant) were calculated. An array-wide Bonferroni-corrected threshold of p < 2 × 10 −7 for ∼250,000 variants was defined as significance level. Ten loci associated with leptin concentrations were found in the total sample, four of which had previously been detected by Kilpeläinen et al. (2016). Only nine variants were detected in Europeans [SNP rs17151919 (LEP) was not observed in Europeans], five and six of which attained significance after study-wide correction based on leptin levels unadjusted and adjusted for BMI, respectively. We used the effect sizes for unadjusted leptin levels determined in Europeans calculated with an additive model for our MR. How much of the total variance is explained by the SNPs (h 2 SNP ) was not reported (Yaghootkar et al., 2020).
We used the results of GWAS and EWAS on leptin levels adjusted for BMI to assess the effect of BMI on leptin levels. First, we calculated MR with all SNPs that associated significantly with unadjusted leptin levels. In the second step, we considered only those SNPs that showed a significant effect on leptin levels even after adjustment for BMI (p < 0.05). As there is no GWAS for BMI-adjusted AN, we continued to use effect sizes for nonadjusted leptin levels here.
In addition, we performed MR with leptin levels in females only as exposure.

Sample Overlap (GWAS/EWAS)
According to the Supplementary tables in Kilpeläinen et al. (2016), a total of 28 studies were included in Stage 1 and Stage 2. Of these studies, 11 studies overlap with the EWAS by Yaghootkar et al. (2020). Considering the number of samples in each study, the sample overlap is 39%. Conversely, in EWAS, the 11 overlapping studies result in a sample overlap of 28%.

Anorexia Nervosa (AN)
We used the latest GWAS for AN (Watson et al., 2019) for both directions of MR analyses, i.e. using the complete summary statistics for the analysis considering AN as outcome variable and for the reverse MR analysis with AN as exposure considering the eight genome-wide significant SNPs for the definition of the genetic instrument. This GWAS included 16,992 patients with AN and 55,525 controls with European ancestry. Case definitions included a lifetime diagnosis of AN based on hospital or registry records, structured clinical interviews or online questionnaires based on standardized criteria (DSM-III-R, DSM-IV, ICD-8, ICD-9 or ICD-10). In the UK Biobank, cases had self-reported a diagnosis of AN. Controls were matched to ancestry. In some control cohorts, individuals were screened for lifetime eating disorders and/or some or all psychiatric disorders. As AN is infrequent, large, unscreened control cohorts were considered suitable. About 88% of the cases are women, whereas the proportion of women and men in the controls was more balanced. Reported SNP heritability was 11-17% (Watson et al., 2019).
In case of unavailability of SNPs in GWAS for the outcome phenotype, we used proxy-SNPs (Hartwig et al., 2016). The SNPs with linkage disequilibrium (LD) of at least r 2 ≥ 0.80 (on the basis of GRCh37. p13, Ensemble version 87, 1,000 genomes: phase 3 version 5 for European ancestry) were extracted from the in silico tool SNIPA (Arnold et al., 2015) (http://www.snipa.org. Accessed between December 2020 and April 2021). Selection criteria for proxy-SNPs were defined: first highest r 2 , second smallest distance to the lead SNP, third no strand-ambiguous or palindromic alleles (Hartwig et al., 2016).

Sample Overlap (Leptin Levels/AN)
There does not appear to be any sample overlap between the GWAS or EWAS on leptin levels and GWAS on AN (Kilpeläinen et al., 2016;Watson et al., 2019;Yaghootkar et al., 2020).

Testing MR Assumptions and Statistical Analysis
To perform an MR study, three main assumptions have to be met (Lawlor et al., 2008;König and Greco, 2018): 1) The genetic instrument must have a true association with exposure. 2) The genetic instrument is independent of potential confounding factors in the relationship between the exposure and the outcome. 3) the outcome is associated with the genetic instrument only through the effect of the exposure (König and Greco, 2018). Therefore, the effect of the genetic instrument on the outcome has to be mediated only by the exposure and there should be no direct effects .
Only the first assumption can be directly tested . To fulfill this assumption, it is recommended to use independent genome-wide associated SNPs as instrumental variable (IV) for exposure (p < 5 × 10 −8 ). We deviated from this threshold in part because some variants used as instruments did not achieve genome wide significance in some subgroups (e.g. Europeans or women) in contrast to the significant findings for the overall study group. Among these SNPs, we selected those with p < 0.05 and calculated the F-statistics in the next step for each of these SNPs using formula F (beta/se) 2 . We excluded the SNPs with F < 10 from the MR (Swerdlow et al., 2016;Burgess et al., 2017b). Several approaches were applied to investigate assumptions two and three. We analysed horizontal pleiotropy by estimating the intercept of Egger's regression. If Eggers intercept is not significantly different from zero, it can be assumed that there is no horizontal pleiotropy (Burgess and Thompson, 2017). In order to detect potential pleiotropic instruments, we also conducted MR PRESSO, which applies a global distortion test, to evaluate whether the removal of the potentially pleiotropic instrument leads to a significant difference in the overall causal estimate. Simulation studies showed that MR PRESSO is more sensitive to horizontal pleiotropy than MR Eggers intercept (Verbanck et al., 2018).
Genetic polymorphisms are sometimes associated with multiple aspects or dimensions of a single trait . To test such heterogeneity of the IV, we used Cochran's Q-statistic (Burgess et al., 2017a). This test examines whether causal estimates of genetic variants (SNPs) are comparative (Bowden et al., 2015).
It is highly implausible that all genetic variants used as IV in MR fulfill the instrumental variable assumptions. The different methods vary in their robustness to violations of the assumptions for MR. There is no method that can provide an infallible test of causality. Therefore, it is recommended to perform different methods to assess whether a causal effect as determined via MR is appropriate (Burgess et al., 2017a;Burgess et al., 2019). We carried out the following statistical methods: 1) Inverse-Variance Weighted method (IVW) assumes that all ratio estimates provide independent evidence on the causal effect and there is no pleiotropic effect. IVW thus assumes, that all genetic variants are valid instruments. There is no intercept term in the regression model (Burgess and Thompson, 2017). 2) MR-Egger estimates the intercept as part of the analyses. The intercept term is interpreted as the average pleiotropic effect of a genetic variant included in the analyses. The pleiotropic effect is the effect of the genetic variant on the outcome that is not mediated via the exposure. A non-zero intercept from MR-Egger shows that there is directional pleiotropy, or that IV assumptions are violated, or both Burgess and Thompson, 2017). If the intercept term is exactly equal to zero, then the MR-Egger estimate will equal the IVW estimate. Alternatively, if the pleiotropic effects are independently distributed from the genetic associations with the risk factor (InSIDE assumption: INstrument Strength Independent of Direct Effect), then the MR-Egger estimate will be a consistent estimate of the causal effect as the sample size and number of genetic variants both increase (Burgess and Thompson, 2017;Hartwig et al., 2017). Comparison of IVW and MR-Egger estimates is helpful to interpret results of a MR analysis (Burgess and Thompson, 2017). 3) The mode-based estimate (MBE; simple mode, weighted mode) is robust to horizontal pleiotropy in a different manner to that of the IVW, MR-Egger or weighted median methods (Hartwig et al., 2017). MBE estimates the true causal effect consistently given the assumption that across all instruments, the most frequent value of bias through pleiotropy is zero [ZEro-modal pleiotropy assumption (ZEMPA) is met] (Hartwig et al., 2017). We calculated both simple as well as weighted MBE. Simple MBE is less precise than weighted MBE, but simple MBE is less prone to bias due to violations of the InSIDE assumption. Comparing both methods is perceived as being useful, but requires care caution, since the simple MBE may in some cases be imprecise (Hartwig et al., 2017). 4) Median based estimators are consistent even when up to 50% of the information comes from invalid instrumental variables. Weighted median estimator has an efficiency similar to that of IVW method, simple median estimator is more inefficient compared with IVW and weighted median methods ; penalised weighted median estimator is more robust in the case of heterogeneity of IV (Rees et al., 2019). 5) Robust Adjusted Profile Score (MR RAPS) provides an overall estimator, which is robust against systematic and idiosyncratic pleiotropy (some genetic instruments have a large effect on outcome) (Zhao et al., 2020). 6) MR PRESSO (see above) (Verbanck et al., 2018).
To investigate the relationship between study accuracy and effect size, we created a funnel plot . Asymmetry in the funnel plot indicates that assumptions for MR are not met (Katikireddi et al., 2018). To assess whether a single SNP had a large effect on the regression coefficients, we performed a leave-one-out approach. To conduct this analysis, we ran the IVW regression by omitting each genetic variant in turn (Burgess and Thompson, 2017).
Forest and scatter plots were used to visualize combined results of single and multi-SNP analyses. The scatter plots show the single SNP effects on the exposure against the single SNP effects on the outcome with corresponding standard deviations and estimated regression lines of the multi-SNP analyses.
The power analysis estimated whether the analysis, given sample size, proportion of cases in the study, and the proportion of variance explained, provides sufficient statistical power to detect a true casual effect (Brion et al., 2013). Because the proportion of variance explained (h 2 SNP ) was not published in the papers on GWAS and EWAS on leptin levels, we calculated these using the LDSC tool (Bulik-Sullivan et al., 2015).

RESULTS
In the analysis with leptin level as exposure (Kilpeläinen et al., 2016) and risk for AN as outcome including five SNPs (Table 1), the results were inconsistent depending on the method of analysis ( Table 2). Simple median, weighted median, penalized weighted median, MR PRESSO and MR RAPS revealed a significant causal effect of low leptin levels on the risk of AN; however, MR Egger, Inverse variance weighted (IVW), simple mode and weighted mode did not (Table 2; Figures 1A,B). Single SNP analyses showed that two SNPs were significant: rs10487505 (LEP) and rs6071166 (SLC32A1). The MR Eggers intercept showed no evidence of pleiotropy, but MR PRESSO identified rs8043757 (FTO) as an outlier and calculated outlier-corrected MR estimates: p 0.008. The heterogeneity tests were significant. The overall effect calculated with the penalised weighted median method, which is robust for heterogeneity, was significant ( Table 2). Since the SNP near LEP rs10487505:C:G is palindromic, we inspected the effect directions for SNPs in high LD (rs791600:G:A and rs2167289:G:T; r 2 0.98); the effect directions for leptin and AN were the same as for rs10487505. Funnel and leave-one-out plots are shown in Supplementary Figures 1,2. The h 2 SNP calculated with the LDSC tool for GWAS on leptin levels was 0.097 (genomic inflation factor λ 1.0574). Our analysis had a power of 80% to detect an OR of 1.08 for AN per 1 standard deviation decrease in leptin levels and 100% power to detect an OR of 1.14 (Supplementary Figure 3).
The GWAS on BMI-adjusted leptin levels showed that rs8043757 was no longer significantly associated with leptin levels (  6-9). In the single SNP analysis, rs10487505 and rs6071166 were again significant. After excluding the SNP rs8043757, the results were similar, except that the p-value from the simple mode method was no longer significant. However, the overall effect sizes were slightly higher (Figure 2, Supplementary Table 3, Supplementary Figures 10,11). In both analyses, there was no evidence of pleiotropy and heterogeneity.
We also used the SNPs from the exome study by (Yaghootkar et al., 2020) as genetic instruments (Table 3). Loci LEP, GCKR, CCNL1 and FTO were detected by both Kilpeläinen et al. (2016) and Yaghootkar et al. (2020), albeit with different SNPs. The loci KLHL31 (lead SNP rs3799260), ZNF800 (lead SNP rs62621812), ACTL9 (lead SNP rs234055) and KLF14 (lead SNP rs972283) had not been detected by (Kilpeläinen et al., 2016) and revealed novel genetic associations with leptin levels. However, for Europeans,  (Kilpeläinen et al., 2016) and their association with the risk for anorexia nervosa (AN) (OR is transformed to beta) (Watson et al., 2019).  SNPs rs3799260 (KLHL31) and rs234055 (ACTL9) do not seem to be particularly relevant, at least they appear to be weak IV (F < 10). MR analysis did not detect an overall effect of leptin levels on the risk for AN. In the single SNP analysis, SNP rs791600 (LEP) was significantly associated with AN (beta -0.905, se 0.316, p 0.004) ( Table 4, Supplementary Figures 12-15). There is no evidence of pleiotropy. Heterogeneity test based on IVW was significant (p 0.029). The h 2 SNP for EWAS on leptin levels calculated with the LDSC tool was 0.1913 (genomic inflation factor λ 1.1612). Power analysis showed that our MR with EWAS data had a power of 80% to detect an OR of 1.057 and a power of 100% to detect an OR of 1.10 for AN per 1 standard deviation decrease in leptin level (Supplementary Figure 16).

SNP
Because the SNP rs1121980 (FTO) was no longer significant in EWAS for leptin levels adjusted for BMI, we repeated the MR without this SNP (Table 3). Again, there was no overall effect of leptin on risk for AN (Table 4, Figures 3, Supplementary  Figures 17,18). Evidence of horizontal pleiotropy or heterogeneity was not found.
The analyses using effect sizes for SNPs associated with leptin levels identified in the EWAS were also performed in females only (Supplementary Table 4). Again, no significant finding emerged. The SNP rs791600 (LEP) was again significant in the single SNP In the reverse MR analyses, there was no causal effect of a higher risk of AN on leptin levels using the GWAS data ( Figures  4, Supplementary Tables 8,9, Supplementary Figures 27,28). There was no evidence of either horizontal pleiotropy or heterogeneity of genetic instruments.
The reverse analysis of the effect of risk of AN on leptin levels using EWAS could not be performed because the SNPs associated with AN (Watson et al., 2019) are not included in EWAS (Yaghootkar et al., 2020). We could only find proxies (LD ≥ 0.8) for two SNPs (rs10747478 and rs9821797), which would not have been sufficient for MR analyses.

Sensitivity Analyses
We performed two sensitivity analyses where we used the loci LEP, GCKR and CCNL1 which were detected in both GWAS and EWAS as IV. We did not include the SNPs rs8043757 and rs1121980 (FTO). In the first MR we used SNPs rs10487505 (LEP), rs780093 (GCKR), rs900400 (CCNL1) detected in GWAS by Kilpeläinen et al. (2016). In the second MR SNPs rs791600 (LEP), rs1260326 (GCKR) and rs900399 (CCNL1) reported by Yaghootkar et al. (2020) were included. These two analyses  For further sensitivity analyses we proceeded as follows: as IV we used SNPs that were reported as significant in EWAS (Yaghootkar et al., 2020), but for MR we used the effect sizes from GWAS (Kilpeläinen et al., 2016) (Yaghootkar et al., 2020) on the risk for anorexia nervosa (AN) (Watson et al., 2019)   Frontiers in Genetics | www.frontiersin.org September 2021 | Volume 12 | Article 733606 Table 12). In this MR analysis, only median-based methods showed a significant effect of leptin levels on the risk of AN (Supplementary Table 13, Supplementary Figures 37-40).

(Supplementary
Other MR methods were not significant. The reverse analysis (SNPs from GWAS on leptin levels with effect sizes from EWAS as IV for exposure on risk for AN) was not possible because only two SNPs (rs780093, rs900400) from GWAS were also reported in the EWAS. In the EWAS there were also no proxies for the other SNPs.

DISCUSSION
Our MR analyses suggest a causal effect of lower leptin levels for risk of AN using data from GWAS (Kilpeläinen et al., 2016;Watson et al., 2019). This also applied upon restriction to leptin SNP data based on females only. We performed our MR analysis with different MR methods to take into account respective strengths and weaknesses of the statistical methods against violation of the assumptions for MR (Burgess et al., 2017a). We also carried out different sensitivity analyses to take into account the specificities of the datasets. Not every MR method showed a significant effect. This is due to properties of the methods with respect to the characteristics of the data used.
In the analysis with leptin levels unadjusted for BMI as exposure and AN as outcome, IV heterogeneity and horizontal pleiotropy seemed to be an issue. The outlier-corrected overall estimator [MR-PRESSO (Verbanck et al., 2018)] was significant after excluding pleiotropic SNP rs8043757 (FTO). MR-RAPS, whose estimator is also more robust to pleiotropy (Zhao et al., 2019), showed a significant causal effect of leptin levels on the risk for AN. A significant effect was also observed upon use of the weighted median methods, of which the penalized weighted median method is more robust to heterogeneity of IV (Rees et al., 2019). These methods show that the weighted majority of genetic variants are associated with the outcome (Burgess and Thompson, 2017). Since pleiotropy is obviously an issue in this analysis, the assumptions for IVW were not met. The overall estimator from MR Egger was not significant, but the direction was positive, unlike the estimators from the analyses with other methods. This confirms the suspicion of the violation of the assumptions for IV, because an influential genetic variant changes the sign of the MR-Egger estimate (Burgess and Thompson, 2017). Upon consideration of all MR results giving the same direction of effect (except for MR-Egger), we interpret the findings using data from GWAS as indicative of an elevated risk for development of AN due to factors associated with low leptin levels.
After excluding pleiotropic SNP rs8043757 (FTO), the difficulties with horizontal pleiotropy and heterogeneity were no longer evident. The estimators of MBE and MR Egger were not significant, however IVW, RAPS, all median-based methods (single, weighted and penalised weighted median) showed that there is a significant causal effect of lower leptin levels on higher risk for AN. The analyses based on females only showed similar results [mode-based methods and MR Egger were not significant, IVW, RAPS, all median-based methods (single, weighted and penalised weighted median) were significant]. We again conclude that our results indicate that lower leptin levels may have a causal effect on higher risk of AN, and this is also independent of BMI.
We were unable to assess the effects of leptin levels on the risk of AN in women compared to men due to lack of data. Such a study would be possible with availability of sex-dependent analyses of AN-GWAS. Such a GWAS including only affected females (or males) and sex-matched healthy controls could help to clarify whether the substantially lower risk for AN in men can also in part be explained by differences in genetic factors. It should be mentioned that in GWAS on AN (Watson et al., 2019) 12% of the cases were men. The GWAS on leptin levels (Kilpeläinen et al., 2016) revealed that only SNP rs10487505 (LEP) resulted in different effect sizes and p-values in women and men (Kilpeläinen et al., 2016). Yaghootkar et al. (2020) reported two genes (CNTD1, DNAJC18), which may show sex-differences in associations with leptin levels. Upon use of the GWAS (Kilpeläinen et al., 2016), our MR study was able to show the effect of leptin levels based on females only on the risk for AN, despite lower power.
Single SNP analyses showed that rs10487505 and rs6071166 were consistently associated with risk for AN in all analyses. The common variant rs10487505 (7:128220110) is located 21 kb from LEP (Kilpeläinen et al., 2016) within a regulatory region. The functional implications of the association of SNP rs10487505 in LEP are not yet understood. MicroRNA 129-1 (7:128207872-128207943) is located upstream to rs10487505. This microRNA has been associated with e.g. glycerol levels (Raitoharju et al., 2014), type II diabetes mellitus (Hara et al., 2014), educational attainment and cognitive function (Lee et al., 2018) and tumorigenesis of retinoblastoma (Zhao et al., 2009). More likely than microRNA 129-1, the interplay of long noncoding RNA (IncOb), enhancer sequences (LE1 and LE2) and LEP promoter transcript, which are also located in this region, could be relevant for the expression of LEP (Dallner et al., 2019). Thus, higher expression of the IncOb was associated with higher leptin mRNA levels in fed, fasted, and obese mice (Dallner et al., 2019). The IncOb is located upstream (−28 kb) of the LE1, which in turn is upstream from LEP. The enhancer LE2 is located downstream from LEP. Both enhancers interact with the proximal leptin promoter, which is required for proper leptin expression (de la Brousse et al., 1996;Mason et al., 1998;Munzberg and Heymsfield, 2019). The proximal promoter alone is not sufficient for adipose-specific leptin expression: at least one enhancer (LE1 or LE2) is required to maintain leptin expression in adipose tissue. IncOb binds exclusively to the promoter. The absence of IncOb significantly decreases leptin gene expression, therefore IncOb seems to be required to exploit the full capacity of leptin expression (Dallner et al., 2019;Munzberg and Heymsfield, 2019). Expression of the analog human lncRNA (EST EL947753) is correlated with leptin mRNA expression in human adipose tissue. The variant rs10487505, which showed the lowest p-value and largest effect size with BMI-adjusted leptin levels, is located in this lncRNA.
The role of leptin levels for risk of AN was previously hypothesized in a GWAS on stringently phenotyped patients with AN that provided hints for a dysregulation of leptin signaling in this eating disorder (Li et al., 2017). However, the detected SNP (rs929626 in EBF1) did not reach genome-wide significance (p 2.04 × 10 −7 ) for AN.
The other consistently significant single SNP, SNP rs6071166 (20q11.23), is an intergenic variant ∼20 kb from SLC32A1 (Kilpeläinen et al., 2016). However, the knockdown of the Adig gene [Adipogenin; alias SMAF1 (Small Adipocyte Factor 1)] in mouse models decreased Lep expression and leptin secretion and thus provided evidence that ADIG controls leptin levels. The ADIG gene is located ∼116 kb from rs6071166 (Kilpeläinen et al., 2016) and plays a role in adipocyte differentiation and development (Kim et al., 2004).
We can only speculate why the analysis with genetic instruments from the EWAS did not show the overall effect of leptin levels on AN. The study by Yaghootkar et al. (2020) analysed the samples with the HumanExome BeadChip, which takes into account variants in gene-coding regions (some versions of the chip apparently also take into account SNPs outside of exons). This study also revealed that the LEP gene has a causal effect on the risk of AN: SNP rs791600, an intergenic variant (intron variant of predicted ncRNA LOC105375494) near the LEP gene, was significantly associated in the single SNP analyses. SNP rs791600 is in linkage disequilibrium (LD) (EUR r 2 0.70) with rs10487505 (see above). Since the EWAS contained SNPs from protein-coding regions and only a few selected intronic/intergenic SNPs, most intergenic SNPs are consequently missing from the study: both rs10487505 and rs6071166 (SLC32A1) were not included in the exome array. In addition to the LEP locus, the GCKR and CCNL1 loci were identified in both the GWAS by Kilpeläinen et al. (2016) and in the EWAS by Yaghootkar et al. (2020), albeit the lead SNPs being different. The EWAS reported four novel loci in Europeans: rs62621812 (p 8.2 × 10 −8 , after adjusting for BMI p 2.8 × 10 −12 , ZNF800, missense), rs379926 (p 3.8 × 10 −3 , after adjusting for BMI p 3.8 × 10 −6 , KLHL31, missense), rs972283 (p 1.1 × 10 −10 , after adjusting for BMI p 3.8 × 10 −18 , KLF14, intergenic), and rs2340550 (only nominally significant after adjusting for BMI p 2.6 × 10 −02 , ACTL9, missense; not included in MR). SNPs rs62621812 and rs972283 revealed non-significant positive associations between leptin level and risk of AN in single-SNP MR. In addition, as already mentioned, EWAS does not include the SNP rs6071166 (SLC32A1), which showed a significant negative effect on the risk of AN in single-SNP analyses. Sensitivity analyses for overlapping loci (LEP, GCKR and CCNL1) from GWAS and EWAS as IV support the suggestion that there is an effect of lower leptin levels on risk for AN. Further sensitivity analysis with SNPs from EWAS as IV, but with effect sizes from GWAS could not find a significant effect with most MR methods. Reverse analysis with SNPs from GWAS, but effect sizes from EWAS was not possible due to lack of data. More powerful GWAS are required to resolve the discrepancy between our results based on these two approaches and to figure out whether our MR results based on the GWAS by Kilpeläinen et al. (2016) are false positive.
For a SNP to be a valid IV, it does not necessarily have to be located in a protein-coding region. A SNP that is appropriate as an IV can also act as a cis, located close to the encoding gene influencing mRNA expression. However, trans-acting SNPs located outside the target gene or even on a different chromosome can also serve as suitable genetic instruments. However, these SNPs are more likely to be pleiotropic (Swerdlow et al., 2016). The more recent EWAS study can be advantageous as a data basis for IV, if they take less frequent variants into account (Swerdlow et al., 2016). For the GWAS by Kilpeläinen et al. (2016), different genotyping platforms included variants with at least MAF >1% or >5%, depending on the cohort. No information could be found on the different exome Frontiers in Genetics | www.frontiersin.org September 2021 | Volume 12 | Article 733606 genotyping arrays used for the genotyping of the cohorts in EWAS (Yaghootkar et al., 2020) that would indicate that SNPs with MAF <1% were taken into account. Thus, there is no evidence regarding this aspect that MR analysis with SNPs from EWAS is by definition superior to MR analysis with genetic tools from GWAS. With regard to the number of subjects, GWAS (Kilpeläinen et al., 2016) and EWAS (Yaghootkar et al., 2020) were similar in size. The power analysis showed that the IV from EWAS had slightly higher power to detect the true effect (EWAS: 100% power to detect OR 1.10 vs GWAS: OR 1.14). It should be pointed out that the genomic inflation rate for EWAS was elevated, which may lead to an overestimation of the h 2 SNP . One reason we did not find an effect with EWAS data, if it exists, could be due to the rank-based inverse normal transformation of the leptin level measurements. Due to this transformation, the beta estimator cannot accurately reflect the effect per allele on the trait.
Our analyses (GWAS and EWAS) showed that the leptin levelreducing alleles of SNPs that are close to LEP associate with a higher risk of AN. The detection of a causal role for AN of variants in linkage disequilibrium with an effect on leptin levels in both GWAS and EWAS is in itself an important finding linking the two phenotypes. Based on our MR analyses, we conclude that low leptin levels may well increase the risk for AN. These results would further support the labeling of AN as a metabo-psychiatric disorder (Watson et al., 2019).
The causal reason for the observed hypoleptinemia in patients with acute AN (Hebebrand et al., 1997;Balligand et al., 1998;Misra and Klibanski, 2014) is not the disease AN per se, but the associated starvation, which entails a decreased body fat. However, based on the results of this study, a low endogenous leptin synthesis may represent a risk factor for the development of AN. Potentially young females [and potentially males, too, albeit with a much lower a priori risk (Nicholls et al., 2011;Steinhausen and Jensen, 2015)] with low leptin levels are at a higher risk for AN due to their genetically lower leptin levels (Kilpeläinen et al., 2016) and/or dysregulated leptin signaling (Li et al., 2017) in addition to other risk factors [e.g. personality traits, emotion dysregulation, hormonal and metabolic effects (Wonderlich et al., 2005;Thompson-Brenner et al., 2008;Miller, 2011;Lavender et al., 2015)]. For such individuals, self-induced weight loss (experimenting with diets, high physical activity/sport) may entail an increased risk of the onset of hypoleptinemia, which could be a relevant factor for patients to become entrapped in AN (Hebebrand et al., 2019). Physical activity has been shown to lead to lower levels of leptin (Lichtenstein et al., 2015), which may also increase the risk of AN in predisposed persons (Davis et al., 1994;Vayalapalli et al., 2018). As AN evolves, leptin levels drop into a critical range as a consequence of loss of body fat due to a negative energy balance. In the state of starvation associated with low body fat and hypoleptinemia, a positive energy balance only slowly entails increments in leptin levels upon re-alimantation (Hebebrand et al., 1997;Balligand et al., 1998), patients remain entrapped in the disease for a prolonged period of time. This slow increase of leptin secretion may underlie the slow improvement of patients during realimentation. The finding that an increase in leptin levels through exogenous application of recombinant human leptin via off-label treatment with human recombinant leptin (metreleptin) allows an escape from the "cage" (Milos et al., 2020;Antel et al., 2021) supports this hypothesis. In conclusion, a genetically lower leptin level in patients may not only represent a risk factor for a more rapid entrenchment in the eating disorder due to the in comparison to non-predisposed individuals earlier development of hypoleptinemia and as a consequence the initiation of central effects of starvation. The genetic predisposition to low leptin levels may also contribute to the slow recovery among realimentation. Finally, this genetic predisposition may play a role in the frequent relapses observed in AN.

LIMITATIONS
When interpreting the results, it is important to note that the studies included in GWAS on leptin levels (Kilpeläinen et al., 2016) were mostly population-based or family-based. The three included case-control studies were not specific to leptin levels (NHS: type 2 diabetes and breast cancer; GEMS: dyslipidemia; Health 2000: the condition for cases not mentioned). The EWAS on leptin levels (Yaghootkar et al., 2020) included populationbased, family-based and cohort studies. Therefore, it is more reasonable to assume that very low leptin levels (and also very high ones) were not or only partially taken into account in the estimation of genetic associations. In MR, a linear relationship is assumed; non-linear relationships cannot be taken into account with these models. The GWAS on leptin levels and AN, as well as EWAS on leptin levels, do not focus on these phenotypes in childhood. Although some or several of the included studies included children or adolescents, it is only possible to apply our results to children and adolescents if the genetic factors are similar to those in adults. Since most patients with AN are female, GWAS on AN in males are lacking. It is possible that the effects of leptin levels on the risk of AN differ between the sexes. Two-sample MR makes the assumption of independent samples. Violation of this assumption leads to inflated type 1 error and biased effect estimates . In a simulation analysis for binary outcomes with sample overlap in the controls -which would be theoretically possible for our study -there were no detectable biases in the IV estimates, not even with extremely weak instruments, and no inflation of the type 1 error .

CONCLUSION
Our MR analysis provides support for a causal effect of lower leptin levels on a higher risk of AN. This holds up upon use of data based on leptin levels in females only. This conclusion is based on the analyses of GWAS data, the analyses with EWAS data did not implicate a causal effect. AN itself has no causal effect on leptin levels, although starvation induced hypoleptinemia characterizes patients with acute AN. Apart from a genetic predisposition to lower leptin levels representing a risk factor