Placenta and appetite genes GDF15 and IGFBP7 are associated with hyperemesis gravidarum

Hyperemesis gravidarum (HG), severe nausea and vomiting of pregnancy, occurs in 0.3–2% of pregnancies and is associated with maternal and fetal morbidity. The cause of HG remains unknown, but familial aggregation and results of twin studies suggest that understanding the genetic contribution is essential for comprehending the disease etiology. Here, we conduct a genome-wide association study (GWAS) for binary (HG) and ordinal (severity of nausea and vomiting) phenotypes of pregnancy complications. Two loci, chr19p13.11 and chr4q12, are genome-wide significant (p < 5 × 10−8) in both association scans and are replicated in an independent cohort. The genes implicated at these two loci are GDF15 and IGFBP7 respectively, both known to be involved in placentation, appetite, and cachexia. While proving the casual roles of GDF15 and IGFBP7 in nausea and vomiting of pregnancy requires further study, this GWAS provides insights into the genetic risk factors contributing to the disease.

N ausea and vomiting of pregnancy (NVP) affects 50-90% of pregnant women 1 and as many as 18% of pregnant women take medication to treat this condition 2 . Hyperemesis gravidarum (HG) is the most severe form and occurs in 0.3-2% of pregnancies 3 . Its clinical presentation includes severe intractable vomiting, often associated with dehydration, weight loss (>5% pre-pregnancy weight), ketonuria, nutritional deficiencies, and electrolyte disturbances 4 . The severity of the disease has been well documented, for instance with the death of author Charlotte Brontë 5 . To this day, HG remains the second leading cause of hospitalization during pregnancy 6 .
Despite the prevalence of NVP and the gravity of HG, decades of research have failed to identify the cause, and a clinically proven, safe, and effective treatment has yet to be found 7 . While the absence of NVP is associated with a higher risk of miscarriage 8 , having the most severe form of nausea and vomiting (HG) is also associated with poor fetal outcomes including preterm birth, neurodevelopmental delay, and vitamin K deficient embryopathy [9][10][11] . Several possible etiological factors have been investigated but the cause is unknown. Support for a genetic component to NVP come from twin studies in Finnish, Norwegian, and Spanish cohorts 12,13 . Heritability estimates for presence of NVP are as high as 73% 13 . Mothers and sisters of women affected by HG are at increased risk compared to mothers and sisters of unaffected women 14 . The risk of HG is increased by 17fold for sisters of HG patients, and there is a >27-fold increased risk of mother-daughter recurrence when a mother has HG with two daughters 15,16 .
The aim of this study is to use human genetics as a steppingstone towards elucidating the causes of HG. Herein, we report a genome-wide association study (GWAS) aimed at identifying variants that influence the risk of HG. We performed two genome-wide association scans, one for a binary HG phenotype, and one for an ordinal phenotype related to HG (severity of NVP). The key findings were replicated in two independent populations of women with clinically defined HG.
Genome-wide association scan of the ordinal phenotype. For SCAN2, we next conducted a genome-wide association scan of the 53,731 unrelated female 23andMe research participants of European descent, using NVP severity as an ordinal response: none (N = 14,988), slight (N = 14,292), moderate (N = 17,786), severe (N = 5445), and very severe NVP (N = 1220), coded on a scale of 0 to 4. Details of the definition of each response based on the survey data can be found in the Methods section. Cohort statistics are shown in Table 3. Manhattan and quantile-quantile plots are shown in Fig. 1b, and SNP-level QC information is shown in Supplementary Table 1c.
The two most strongly associated loci in SCAN2 were also genome-wide significant in SCAN1: chr19p13. 11 Figure 3) and there were eight additional association signals in SCAN2 with p < 10 −6 (Supplementary Table 2).
Replication. Two loci were genome-wide significant in both association scans and were selected for replication genotyping. Given that the replication cohorts contain information on HG as a binary trait and that HG is the more clinically relevant phenotype, three additional loci with p < 10 −6 in the binary trait scan (SCAN1) were also selected for replication. The first replication cohort (HG-IV) included 789 women with HG requiring IV fluid treatment and 606 controls reporting normal NVP. The second replication cohort (HG-TPN) included only women at the extreme ends of the clinical spectrum, 110 women requiring total parenteral nutrition (TPN) and 143 women reporting no NVP. Basic demographics of the replication cohorts are shown in Table 5. Over 90% of participants self-reported that they were of European descent.
We selected rs16982345 in the chr19p13.11 locus as a proxy for the lead SNP rs45543339 (r 2 = 0.98, D′ = 0.99) for replication in HG-IV (because rs16982345 had a commercially available high quality SNP assay using the Taqman genotyping system). rs16982345 was successfully genotyped in 789 individuals with clinically defined HG and 581 controls reporting normal NVP using the TaqMan genotyping platform ( We also successfully genotyped the variant in 103 women with HG who required TPN and 136 women who reported no NVP in at least two pregnancies (HG-TPN). The genotyping results had    (Table 6). We also performed replication of the second genome-wide significant locus. rs4865234 was used as a proxy for rs143409503 (r 2 = 0.95, D′ = 1.0) (because rs143409503 is an indel, which cannot be assayed using the TaqMan genotyping system). rs4865234 replicated in HG-IV (p = 3.5 × 10 −4 , OR = 1.35 [1.14, 1.59]) ( Table 6). The genotyping results of DNA from HG-TPN were also supportive (p = 2.8 × 10 −3 , OR = 1.81 [1.21, 2.73]).
Combined analysis. For the replicated loci, we combined the replication results with SCAN1 using a fixed-effect meta-analysis (Table 7). Briefly, this test computes a weighted average of oddsratio estimates across studies while accounting for standard errors. For rs16982345 on chr19p13.11, we estimated a meta-OR of 1.50 [1.38, 1.65] (p = 2.12 × 10 −19 ). Heterogeneity of individual effect sizes did not appear to influence the combined estimate (p het = 0.66). For rs4865234 on chr4q12, the meta-OR was 1.33 [1.23, 1.62] (p = 9.29 × 10 −12 ) with little evidence of effect size heterogeneity (p het = 0.727).
Conditional analyses. In order to identify additional association signals, we performed stepwise conditional analysis at all loci that reached genome-wide significance in either SCAN1 or SCAN2. Conditional analysis detected secondary associations in the chr19p13.11 locus in SCAN1 and in the chr19p13.11 and PGR/ TRPC6 loci in SCAN2 (Supplementary Table 3 and Supplementary Figures 4 and 5). LD between the primary and secondary signals is low (r 2 < 0.01), suggesting that they are independent signals corresponding to multiple risk variants within the locus.
The genomic context of the replicated loci chr19p13.11 and chr4q12 was analyzed using HaploReg 17 (Supplementary Data 1). Within those two regions, we analyzed A) the index SNPs for strongest association in SCAN1 (rs45543339 and rs143409503) and the proxy SNPs that were confirmed in B) the independent replication cohort (rs16982345 and rs4865234).
The lead SNP at chr19p13.11, rs45543339, was located in an intron of LRRC25, and the only notable functional annotation was an eQTL for KCNN1 in mucosa of the esophagus. rs16982345, which is within the 99% credible set at this locus and is in tight LD (r 2 = 0.98) with the lead SNP rs45543339, is located within an intron of GDF15 (Supplementary Data 1 and 2). It overlapped enhancer histone marks in three tissues, including placenta, and DNAse hypersensitivity sites in nine tissues, including ovary. A summary of the variants in LD with the lead SNP and associated with altered GDF15 expression is shown in Supplementary Data 2 [17][18][19][20][39][40][41][42]44 , with the full list of variants in LD with the lead SNP identified using HaploReg listed in Supplementary Data 1. One of the variants associated with HG in SCAN1, rs17725099 (p = 1.03 × 10 −13 ) and in LD with rs16982345 (r 2 = 0.77), was the top association signal (β = 0.16, Table 6 Replication results based on Fisher's exact test for the two most significantly associated loci in SCAN1 and SCAN2: rs16982345 (chr19p13.11), rs4865234 (chr4q12), and the three association signals in SCAN1 with p > 5 × 10 −8 and p < 1 × 10 −6 , rs5754397 (SYN3) rs78353059 (MMADHC/RND3), and rs56108151 (TMEM38B/ZNF462) The number (N) for each group is the total number of samples that were successfully genotyped. The call rate was >95%. OR estimates, 95% confidence intervals, and p-values were computed using 2 × 2 contingency tables in R. The effect allele is assumed to be the allele in the left-most (i.e., first) homozygous cell  18 . The SNP rs17725099 was also associated with circulating levels of GDF15 in the first study to report on the heritability of GDF15 plasma concentration 19 . In addition, rs1058587, a SNP in perfect LD (r 2 = 1) with the lead SNP and a part of the 99% credible set at this locus in SCAN1, has previously been associated with altered GDF15 expression. Expression levels of GDF15 were compared in an in vitro model between the cell line DU145 transfected with wild-type GDF15 (rs1058587, C allele) and GDF15 (rs1058587, G allele) 20 . The C allele of rs1058587 was associated with increased GDF15 levels in the in vitro model and was also the risk allele associated with HG, providing evidence that the direction of effect is toward higher levels of GDF15 for HG compared to controls. The lead indel at chr4q12, rs143409503, mapped to an intergenic region with the closest protein-coding gene being IGFBP7,~380 kbp away from the marker. No notable features were identified using HaploReg 17 , the GTEx portal, nor ExSNP 21 databases. Other phenotypes associated with SNPs in GDF15 and IGFBP7 are reported in Supplementary Note 1 [41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56] .

Discussion
This GWAS of HG identified two genome-wide significant signals that were subsequently replicated in the larger of two independent cohorts (HG-IV). The most significantly associated locus on chr19p13.11 contained genes GDF15 and LRRC25. While we cannot be certain which gene or genes are implicated, SNPs in LD with the associated variants at the chr19p13.11 locus were found to be associated with altered expression of GDF15. In addition, the gene encoding the receptor for GDF15, the brainstemrestricted receptor GFRAL 22 , was associated with NVP in SCAN2 (Table 4) adding further evidence supporting a previously unknown biological connection between GDF15 and HG. While the results of this study do not establish a causal link between GDF15 and HG, the association between this gene and HG is of particular importance because it highlights the possibility of a pathway involved in the etiology of the condition. GDF15 encodes a TGF-β superfamily member that is expressed at its highest levels in the trophoblast cells of the placenta 23 . The protein is found in maternal serum and increases significantly in the first two trimesters 23 . GDF15 is believed to suppress production of proinflammatory cytokines in order to facilitate placentation and maintain pregnancy 24 . In addition to its role in pregnancy, GDF15 has been shown to be a regulator of physiological body weight and appetite via activation of neurons in the hypothalamus and area postrema (vomiting center) of the brainstem 25,26 . It is also notable that abnormal overproduction of GDF15 in cancer was recently found to be the key driver of cancer anorexia and cachexia which, like HG, exhibits symptoms of chronic nausea and weight loss 27,28 . Of particular clinical interest, inhibition of GDF15 restored appetite and weight gain in a mouse model of cancer cachexia 28 , suggesting a therapeutic strategy that may be applicable to patients with HG, if GDF15 proves to be the implicated gene.
The other locus that reached genome-wide significance and was confirmed in the independent replication cohort, is chr4q12. Although there were no genes that can be directly related to the association signal, the protein-coding gene closest to the lead variant is IGFBP7. We have not found eQTLs that link the variants associated in this study to expression of IGFBP7, however, this can be explained at least in part by a lack of eQTL data in relevant tissues during pregnancy. The notably similar and functionally relevant roles of IGFBP7 and GDF15 (in placentation, cachexia, and feeding behavior) make a compelling argument for IGFBP7 despite the lack of evidence linking the association signal to its expression. Insulin-like growth factor binding protein 7 (IGFBP7) is a gene involved in implantation and decidualization of the pregnant uterus 29 . Like GDF15, IGFBP7 is significantly upregulated after implantation and highly expressed in the developing placenta 30 . Inhibition of IGFBP7 caused pregnancy loss in a mouse model by shifting uterine cytokines to Th1 type dominance and repressing uterine decidualization 30 . Thus IGFBP7 may play roles in both miscarriage and in the severity of NVP, providing a genetic mechanism for the protective effect of NVP. In addition, like GDF15, evidence points to IGFBP7 as a promising biomarker for cachexia associated with end-stage disease 29 . An intriguing finding is that the Drosophila homolog of IGFBP7 has been suggested to play a role in neuronal coordination between metabolic status and feeding behavior, potentially, like with HG, causing food aversion to normally attractive food, even when starving, and vice versa 31 .
The strengths of this study include the well-defined phenotype and large sample sizes in two scans and the HG-IV replication sample. The strong independently replicated association signals, and the fact that both candidate genes have roles in early pregnancy when HG sets in, support these genes as functional as well as positional candidate genes 30,32 . However, it is important to replicate the findings reported herein, both using larger replication sample sizes (for HG-TPN vs no NVP) and in studies of different ethnicities to determine whether the findings can be validated and generalized to other populations. Another caveat is that while the conditional analysis yielded additional association signals at chr19p13.11 that may act independently of the lead association signal, they were not replicated. Replication is necessary to determine if additional independent variants are associated at this locus. Finally, the lack of eQTL data and considerable distance between the chr4q12 association signal and the closest gene IGFBP7 is a noteworthy limitation.
Future research should focus on functional follow-up. Since GDF15 and IGFBP7 levels are upregulated during placentation and cachexia 23,24,28,29 , and downregulated prior to miscarriage 30,32 , serum concentrations of GDF15 and IGFBP7 should now be studied in pregnant women with and without HG. If GDF15 or IGFBP7 prove to be relevant, it is conceivable that drugs targeting these proteins may have clinical utility in treating HG. Such a study may also lead to new methods for prediction and diagnosis. In addition, the SCAN2 results suggest a great potential for discovering additional loci associated with HG or NVP. Finally, the findings herein suggest an answer to an age-old paradox. HG can lead to prolonged dehydration and undernutrition, which can be detrimental to maternal and fetal health and can decrease reproductive fitness. The dual roles of GDF15 and IGFBP7 in maintaining pregnancy and in increasing the risk of HG may provide a molecular explanation for why NVP still exists in nature.

Methods
Human subjects for GWAS. GWAS participants were customers of 23andMe, Inc. who consented to participate in research, and provided answers to the morning sickness-related questions. Two genome-wide association scans were performed, for (1) a binary HG phenotype and (2) an ordinal phenotype related to HG. Phenotype definitions are described below. Due to the nature of the phenotype definitions, some research participants were included in both association scans. All research participants included in the analysis provided informed consent and answered on-line surveys according to a human subjects protocol approved by Ethical & Independent Review Services, a private institutional review board.
Phenotype definitions. For the binary HG phenotype analyzed in SCAN1, we compared two ends of the clinical spectrum of NVP: 1306 research participants who reported via an on-line survey that they received IV therapy for NVP were classified as HG cases and 15,756 participants who reported no NVP served as controls. The phenotype definition was ascertained using the interview questions listed in Supplementary Note 2. For the ordinal phenotype related to HG analyzed in SCAN2, data were pulled from four questions also listed in Supplementary Note 2.
Genotyping and imputation and association testing. A total of 17,062 females in SCAN1 and 53,731 females in SCAN2 were genotyped on one of four custom Illumina genotyping arrays and additional genotypes were imputed using the September 2013 release of the 1000 Genomes Project Phase 1 33 reference haplotypes as described previously 34,35 . 16,165 individuals in SCAN1 were also included in SCAN2. Breakdown of the number of participants in specific phenotype categories in SCAN1 and SCAN2 are shown in Supplementary Table 1a. Fields that contain 5 or fewer individuals have been masked to protect the privacy of 23andMe customers. All females were filtered to select for European ancestry, and close relatives were removed. We performed logistic (SCAN1) and linear (SCAN2) regression assuming an additive model for allelic effects on NVP, using age and five principal components of genetic ancestry as covariates. Our previously published analysis of categorical phenotypes using ordinal and linear regression showed high concordance between resulting p-values 36 . Due to comparative ease of implementation and lower computational burden, we used linear regression in SCAN2. The SNP-level quality control information is shown in Supplementary Tables 1b  and 1c. Test statistics were further adjusted for the genomic control inflation factor of λ = 1.044 for SCAN1 and λ = 1.087 for SCAN2. The index SNPs were identified by choosing the SNPs with the strongest association in each associated region. Each region contained SNPs with p < 10 −5 that were grouped into intervals separated by a gap of a minimum of 250 kbp. The SNP with the smallest p-value within each interval was chosen as the lead SNP.
Credible sets. To account for the uncertainty in identifying the true causal variant within an association signal, we computed 99% credible sets at each locus. Under the assumptions that there is a single causal variant at a locus and that this variant was genotyped in the association study, 99% credible sets are defined as smallest sets of variants that are 99% likely, based on posterior probability, to contain the causal variant 37 . While these assumptions are not always true in practice, credible sets provide a helpful summary of the available evidence that one of the SNPs contained in them is causal, and may help in fine-mapping the association signal.
Replication cohorts. The first replication cohort (HG-IV) included 789 HG cases treated with IVs and 606 controls with NVP that did not require treatment. The second replication cohort (HG-TPN) included 110 cases requiring TPN due to HG and 143 controls reporting no NVP in any of at least two pregnancies. All participants gave informed consent. This study was approved by the UCLA Institutional Review Board.
Phenotype definition in the replication cohort. Affected individuals were included if they had an HG diagnosis and were treated with IV fluids and/or TPN. Participants affected by HG recruited acquaintances that were non-blood related and had at least two pregnancies lasting more than 27 weeks gestation. Controls were eligible to participate in the study if they reported normal or no NVP, did not lose weight due to nausea/vomiting, and did not receive medical attention due to NVP.
Recruitment of the replication cohort. For the replication study, the source population for HG cases included patients residing in the US. The majority of cases were recruited through a posting from 2007 to 2017 on the Hyperemesis Education and Research Foundation website (www.HelpHer.org). Minors (under 18 years) were excluded because the risk to benefit ratio to control minors would be difficult to justify and the study requires controls to have had two pregnancies, which is uncommon in minors. Women with pregnancies with chromosome abnormalities and multiple gestations were excluded due to possible distinct etiologies for NVP in these cases. Medical records and recruitment of an acquaintance with at least two pregnancies to serve as a control, were requested of each case. All participants were required to go over an information sheet by phone and return a signed information sheet with all elements of consent in order to enroll in the study.
Sample collection. Saliva samples were collected for DNA analysis from all cases and controls. DNA Genotek saliva kits (Oragene, Ottawa, Canada) were mailed to all cases and controls. The saliva collection kit is self-administered and comes with directions for submitting 2 ml of saliva into a collection vial and returning the sample to the study site via an addressed and postage-paid return envelope provided with the collection kit.
DNA extraction. DNA was extracted from the saliva samples according to manufacturer's instructions (Oragene, Ottawa Canada). Using the kit, we have successfully isolated, on average, 197 µg of DNA of high quality (260/280 1.84) from 2 ml of saliva. The low end of expected DNA quantity reported by the manufacturer is 30 µg/ml of saliva or 60 µg/sample. After the extraction, the DNA was stored at −20°C.
Quality control for replication data. The TaqMan genotyping platform was performed on 384-well plates with a minimum of two blank samples per plate and a minimum of two duplicate samples per plate. Once genotypes were determined from the first 384-well plate, a minimum of three positive controls or one positive control for each genotype was added to the remaining plates. The minimum call rate for each SNP was >95%.
Genotyping and statistical analysis in the replication study. rs16982345, the SNP with the strongest association to NVP in the chr19p13.11 locus was selected for genotyping in the replication sample HG-IV consisting of 789 cases and 606 controls. For the non-coding marker rs16982345 in chr19p13.11, genotyping a minimum of 780 cases and 580 controls, and assuming a proportion of the GG genotype similar to SCAN1 (0.53 for controls and 0.64 for cases), we estimated a power greater than 97% to detect between-group genotype differences using a twosample, two-sided test of proportions at significance level alpha = 0.05.
In addition, we genotyped an extreme cohort (HG-TPN) consisting of 110 women who reported TPN for HG and 143 controls with no NVP in two or more pregnancies. TaqMan genotyping primers for rs16982345 were available from Thermo Fisher Scientific and Applied Biosystems PRISM 7900HT Sequence Detection System (TaqMan) was used for large-scale screening. The call rate was >95%.
Genotypes were tested for significant association with HG using Fisher's exact test using 2 × 2 contingency tables counting the number effect vs non-effect alleles for cases and controls using R.
The SNP rs4865234 was genotyped as a proxy for the most significant association signal at chr4q12 because the most significant association signal, rs143409503, is an indel, which cannot be assayed using the TaqMan genotyping system. The SNP rs4865234 is in tight LD with rs143409503 (r 2 = 0.95, D′ = 1.0) and genotyped on DNA isolated from HG-IV and HG-TPN cohorts as stated previously for rs16982345.
The 3 additional loci that reached significance p < 10 −6 in the binary trait scan (SCAN1) were also genotyped in the two replication cohorts (HG-IV and HG-TPN) using the TaqMan genotyping platform. The SNP rs5754397 was used as a proxy for the lead SNP rs5998706 (r 2 = 0.55, D′ = 0.83), and the SNP rs78353059 was used as a proxy for the lead SNP rs114571265 (r 2 = 0.87, D′ = 1.0) because they were the closest SNPs listed in HaploReg 17 to the lead SNP with an available commercially validated assay from Thermo Fisher Scientific and they also reached significance p < 10 −6 in SCAN1. The lead SNP rs56108151 was commercially available and used for genotyping the two replication cohorts.
Meta-analysis of SCAN1 and replication. We performed a fixed-effect metaanalysis using the R language and software meta to compute the weighted average odds-ratio across the replication group and SCAN1 for sites rs16982345 and rs4865234. The fixed-effect meta-analysis assumes that there exists a true underlying effect in the population and combines individual study estimates (e.g., replication and SCAN1) to obtain a meta estimate. This assumption may be invalid in practice due to differences in linkage disequilibrium patterns across studies, or differences in covariates included. To test for hetereogeneity in effects across studies, or significant study-wide variance, we estimated Cochran's Q 38 for each SNP. This measures the variance across studies using the meta estimate as the true mean, which can be tested in a null-hypothesis framework under a χ 2 nÀ1 distribution. Here, the number of studies is n ¼ 2, so our heterogeneity p-values reflect a onetailed test under a χ 2 1 distribution.
Conditional analysis. The stepwise conditional analysis detects statistically independent signals in a GWAS locus. At each step, we assessed association between variants within 20 kbp around a GWAS locus and the HG or related ordinal phenotype, using top SNPs from preceding steps as additional covariates. This was repeated until no significant association was detected at p < 10 −5 . Finally, we fit a joint model with all of the primary and secondary signals and evaluated its goodness of fit through the likelihood ratio test (LRT) against the model with just the primary signal. The reported effect sizes of the secondary signals are from the joint model, and the p-values are from the likelihood ratio test comparing the full model with the leave-one-out models.
In silico functional follow-up. We used the HaploReg 17 v4.1 tool to identify SNPs in the haplotype blocks containing the two significant GWAS association signals, their LD r 2 and D′ scores, associated alleles, allele frequencies, genes, functional annotation, conservation, regulatory elements, protein-binding sites, motif changes, and eQTLs. We queried the two most significant loci identified, chr19p13.11 and chr4q12. Within those two regions, we analyzed SNPs with HaploReg, which uses LD information from the 1000 Genomes Project 33  We also queried the GTEx portal and ExSNP integrated eQTL databases 21 for the top candidate genes (based on loci identified in this study and similar function in early pregnancy reported previously 30,32 , GDF15 and IGFBP7) to identify eQTLs (17 and 28 studies, respectively) included in the databases. The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS. The data used for the analyses described in this manuscript were obtained from (https://gtexportal.org/home/bubbleHeatmapPage/) the GTEx Portal on 05 May 2017. The SNPs identified were cross checked with the SNPs identified in haplogroups by HaploReg.
Functional annotation of all SNPs (SNPs r 2 > 0.2 with SNPs in the two loci that were found to be in coding regions by HaploReg were queried for all predicted consequences using the 1000 Genomes Project 33 . Finally, published studies reporting SNPs linked to altered expression of the top candidate genes (GDF15 and IGFBP7) and SNPs in the top candidate genes reported in association studies in the GWAS catalog (http://www.ebi.ac.uk/gwas/ accessed on May 2017) and/or in published studies were cross checked with the SNPs identified in HaploReg to determine whether they are in LD with the SNPs linked to HG in our study.
Code availability. All relevant software, including version details, used for the primary GWAS and conditional analyses are described in the references included in the Methods section. All relevant codes with respect to the replication are available from the author M.S.F. (mfejzo@mednet.ucla.edu).
Data availability. Qualified researchers can contact apply.research@23andMe.com to gain access to full GWAS summary statistics following an agreement with 23andMe that protects 23andMe participant privacy. All relevant replication data are available from the author M.S.F. (mfejzo@mednet.ucla.edu).