HLA-DRB1 and HLA-DQB1 genetic diversity modulates response to lithium in bipolar affective disorders

Bipolar affective disorder (BD) is a severe psychiatric illness, for which lithium (Li) is the gold standard for acute and maintenance therapies. The therapeutic response to Li in BD is heterogeneous and reliable biomarkers allowing patients stratification are still needed. A GWAS performed by the International Consortium on Lithium Genetics (ConLiGen) has recently identified genetic markers associated with treatment responses to Li in the human leukocyte antigens (HLA) region. To better understand the molecular mechanisms underlying this association, we have genetically imputed the classical alleles of the HLA region in the European patients of the ConLiGen cohort. We found our best signal for amino-acid variants belonging to the HLA-DRB1*11:01 classical allele, associated with a better response to Li (p < 1 × 10−3; FDR < 0.09 in the recessive model). Alanine or Leucine at position 74 of the HLA-DRB1 heavy chain was associated with a good response while Arginine or Glutamic acid with a poor response. As these variants have been implicated in common inflammatory/autoimmune processes, our findings strongly suggest that HLA-mediated low inflammatory background may contribute to the efficient response to Li in BD patients, while an inflammatory status overriding Li anti-inflammatory properties would favor a weak response.

Bipolar affective disorder (BD) is a severe and often chronic psychiatric illness, characterized by recurrent dysregulation of mood with alternating episodes of mania and depression. BD affects around 48.8 million people 1 accounting for 9.9 million years of life lived with disability worldwide 1 . All-cause mortality and risk of suicide 2 are substantially increased in people with the disorder. Both genetic and environmental factors have been identified to contribute to the pathogenesis of BD 3 but the underlying molecular processes remain poorly understood. Nevertheless, well codified mood stabilizer-based treatments have totally changed BD patient management.
In this context, Lithium (Li) has gained a status as the 'gold standard' amongst mood-stabilizers used for both acute-and maintenance therapy in BD [4][5][6] . Li display strong anti-manic properties 7,8 . Accordingly, it is the only mood stabilizer with protective effects against further episodes of both manic and depressive polarity 6 , making it the most efficient therapeutic option for preventing relapses 9 . Even if Li is recommended as a first-line option for anti-manic and maintenance treatment [10][11][12][13] , its routinely use is hampered by several caveats. Firstly, the therapeutic response to Li in BD is highly heterogeneous. In acute mania, about 65% of patients respond at least partially to Li monotherapy while 35% are refractory 14,15 . In maintenance treatment, an efficient long-term response is reported only for approximately 30% of patients 16 , whereas 30% has intermediate outcomes and 30% poor response. Secondly, Li is toxic at high doses and interindividual variation in metabolism requires a stringent monitoring of Li plasma levels 17 . And thirdly, Li has long-term adverse effects including alteration of thyroid and parathyroid functions and an elevated risk of renal failure 5 .
Amongst biological Li modulator markers, genetic variants (i.e., single nucleotide polymorphisms, SNPs 16,18 ) and genetic traits (i.e., polygenic risk scores, PRS) have been identified that differentiate Li response groups. We recently reported that BD patients with high genetic load for schizophrenia and major depressive disorders (MDD) are less likely to have favorable Li treatment outcomes than those with lower genetic load of these traits 19,20 . Moreover, we showed in cross-trait meta-GWAS and pathway analyses that genetic signals tagging the major histocompatibility complex (MHC) and pro-inflammatory cytokines including tumor necrosis factor alpha (TNF-α), interleukin 4 (IL-4) and interferon gamma (IFNγ) could have a biological role in Li treatment response in BD 19 . More specifically, the most significant finding of the cross-trait GWAS and the SNP from the post-GWAS functional analyses suggested that the human leukocyte antigen (HLA) cluster which is hosted by the major histocompatibility complex (MHC) region could be implicated in the genetic susceptibility to schizophrenia (SCZ) and Li treatment response 19 . The MHC/HLA region is the most robust immunogenetic finding in SCZ 21 , and could be marking a schizophrenia-type pathogenesis associated with Li non-response. The extensive linkage disequilibrium (LD) in the HLA region, the presence of non-HLA genes surrounding the cluster and pitfalls in study designs may explain previous inconsistent data concerning HLA association with Li responsiveness in www.nature.com/scientificreports/ BD [22][23][24] . Nevertheless, the above-mentioned tags of both HLA and pro-inflammatory cytokines reasonably focus attention towards the possible implication of immune processes in responses to Li. Accordingly and given the role of HLA alleles in antigen presentation to downstream effector immune cell subsets, we took advantage of the above-mentioned GWAS 19 for an in-depth analysis of genetic variants of the HLA region in relation to Li treatment response. To this end, we imputed the HLA class-I and class-II genetic variants corresponding to HLA amino-acids and classical alleles and investigated their distribution in responders and non-responders to Li treatment. Table 1. Phenotypic characteristics of individuals used for the analyses. Data are n, n (%), or mean (SD). Alda scale refers to the Retrospective Criteria of Long-Term Treatment Response in Research Subjects with Bipolar Disorder; Non-resp = Li non responders (ALDA subscale A less than or equal to 3); Resp = Li responders (ALDA subscale A greater than or equal to 7).

Results
Genotyping data from 2139 BD patients who had long-term Li treatment were obtained from the International Consortium ConLiGen for HLA imputation (847 in GWAS1 and 1292 in GWAS2). The descriptive statistics of the phenotypes of the participants in each analysis are presented in Table 1. The main goal of the present study was Table 3. Results of the meta-analysis for continuous and dichotomous response outcome in additive model. AF: allele frequency; Q-value: the q-value gives the expected positive false discovery rate; Non-resp = Li non responders (ALDA subscale A less than or equal to 3); Resp = Li responders (ALDA subscale A greater than or equal to 7). www.nature.com/scientificreports/ to uncover HLA amino-acids or classical HLA alleles potentially associated with responses to Li in BD. We first performed association analyses with HLA variants separately for the GWAS1 and GWAS2, and then performed a meta-analysis of the two studies.
For the association analyses of both GWAS1 and GWAS2, no HLA variant reached the Bonferroni significance (threshold of 5.2 × 10 −5 ) in the additive model, neither for the continuous nor for the extreme phenotypes used to measure Li response with the Alda scale (comparison of extreme phenotypes, see "Materials and methods"). Nevertheless, several suggestive associations exhibiting a FDR < 0.25 were observed and are presented in Supplementary Table 2.
In the analysis of GWAS1 for the continuous trait (A score, individuals with B score > 4 were removed), the top ranked variants (p-values < 1 × 10 −3 FDR = 0.12) were located within the HLA-DRB1 and HLA-DQB1 genes. For HLA-DRB1, Tyrosine/Leucine at position 37, Arginine at position 71, and Phenylalanine at position 67 were associated with a better response to Li (Table 2). Interestingly, these amino-acids notably belong to the HLA-DRB1 heavy chain encoding the HLA-DRB1*11:01 allele. Whereas weak response to Li was found to be associated with Leucine at position 26 of the HLA-DQB1 locus ( Table 2). None of these signals were replicated in the GWAS2. The amino-acid frequencies in the GWAS1 and GWAS2 were in agreement with those published in the reference cohorts from Allele frequency net database 25 (Supplementary Table 3).
In a second step, the meta-analysis of GWAS1 and GWAS2 (in the additive model) identified two signals almost reaching the 0.25-FDR threshold for suggestive associations (Table 3). Importantly, the associations were of similar magnitude and direction of effect for both the continuous and dichotomous traits. Once again, the two types of signals were related to the HLA-DRB1 and HLA-DQB1 loci and the frequency of their variants were also comparable to those of the reference cohorts from Allele frequency net database 25 (Supplementary Table 3). The first signal in the HLA-DQB1 locus which consisted of several amino-acids belonging to the HLA-DQB1*02 heavy chain was associated with no response to Li ( Table 3). The second signal tags the amino-acid 74 of the HLA-DRB1 heavy chain. While Alanine or Leucine at position 74 are related to a better response to Li, Arginine or Glutamic acid at position 74 were associated to a weaker one (Table 3). For the analysis of the dichotomous phenotypes (defined by the Alda score ≥ 7 for the responders and ≤ 3 for non-responders), the frequency of the amino-acids Alanine or Leucine at position 74 was higher among the responders as compared to the non responders [72% vs 65% (GWAS1) and 73% vs 67% (GWAS2) in responders vs non-responder respectively]. Since the distribution of HLA alleles is known to exhibit South-North and/or East-West gradients even in European countries, we have specifically analysed the distribution of the alleles (A/L vs others) at position 74 in each country separately. For each country, we thus compared the allelic frequency under a dominant (patients carrying at least a A or a L allele) or a recessive (patients carrying only A or L alleles) model of the dichotomous phenotype of responders and non-responders (when N ≥ 15 in each group). Surprisingly, we observed that, in nine countries out of nine with sufficient number of patients in each dichotomous phenotype (Alda score ≥ 7 and Alda score ≤ 3), the frequency of patients carrying only A or L allele in the recessive model is systematically lower in non-responders than in responders (Table 4), suggesting a recessive effect. The probability of getting such a systematic difference in 9 countries out of 9 by chance is statistically very weak (p = 0.004). When comparing extreme groups (Alda score ≥ 9 and Alda score ≤ 3, see "Materials and methods") in GWAS1 and GWAS2 for the allele A/L at position 74 in the recessive model, we obtained a FDR of 0.09 (Supplementary Tables 1 and 5), stronger than the one obtained for the additive model and confirming the recessive effect. Finally, it is worth mentioning that the HLA-DRB1*11:01 allele bears an Alanine 74 which confirms the above-mentioned favourable effect of such variant on responses to Li. Additionally, it is important to note that the HLA-DRB1 position 74 was previously reported to be in LD with the positions 37, 67, and 71 that came out strongly in GWAS1 (Table 2) 26 .
In order to test dependence between the observed HLA-DRB1 and HLA-DQB1 signals, we evaluated the LD between their variants. We found a high level of correlation (r or D' > 0.9) between HLA-DQB1*02 and two variants of HLA-DRB1 (Supplementary Table 4). The r 2 was found to be the higher (r 2 > 0.5), between Alanine or Leucine at position 74 of HLA-DRB1 and HLA-DQB1*02 (Supplementary Table 4), which suggests a putative link between the HLA-DRB1 and HLA-DQB1. We used Haploview 27 to determine the haplotype between Alanine or Leucine at position 74 of HLA-DRB1 and HLA-DQB1*02 and we observed that all individuals carrying Alanine or Leucine at position 74 of HLA-DRB1 did not carry HLA-DQB1*02 and reciprocally. Given that Alanine and Leucine 74 have been found to be associated with good responders, while HLA-DQB1*02 with non-responders, it was not possible to determine whether the signal was due to the HLA-DRB1 or to the HLA-DQB1, or both of them.

Discussion
Understanding the mechanisms underlying inter-individual responses for a given therapeutic compound is a leading step towards precision medicine as, whatever the treatment efficacy, a subset of patients does not respond. Regarding BD, Li constitute the cornerstone of long-term therapy and management of the disease, but only around one-third of Li treated patients respond positively 28 . A great deal of research has thus been done to identify predictors of the degree of prevention against recurrences of mania and depression during long-term treatment. While several studies showed that Li response was polygenic, a large GWAS-based survey performed by the ConLiGen consortium in 2017, demonstrated that responses to Li were influenced both by the allostatic load of inter-individual SCZ polygenic scores and by genetic markers tagging the HLA antigen presentation pathway and genes encoding pro-inflammatory cytokines 19 .
Li has been demonstrated to exhibit potent anti-inflammatory properties mediated at least by the wellknown inhibitory effects exerted on the Glycogen synthetase kinase-3-beta (GSK-3β) enzyme and downstream on the NF-κB pathway with consequent reduction of the pro-inflammatory cytokines TNF-α and IL-6 levels among others 29,30 . In addition, it has a potent anti-inflammatory action mediated both by the inhibition of IL-1β www.nature.com/scientificreports/ production and by the reduction of cyclooxygenase-2 expression while enhancing expression of IL-2 and IL-10 anti-inflammatory mediators. Finally, Li is also thought to modulate glial cellular inflammatory processes 31,32 . Given that the HLA system is central for the control of inflammatory processes, we imputed putative functional HLA variants within individuals from the above-mentioned GWAS and analyzed their potential association with response to Li. Although no HLA variant per se reached the stringent Bonferroni significance threshold, we observed that several amino-acid positions had a low FDR score. Interestingly, all belonged to the HLA class-II cluster, suggesting the involvement of immune humoral processes in response to Li especially as some of them are located within the 3rd hypervariable (3HVR) sequences (amino-acids 67-74). Of interest, the 3HVR is the primary T-cell recognition site which conditions the immune adaptive processes mediated by the HLA class-II alleles 33 . Accordingly, the 3HVR is the core-associated region with immune processes including chronic inflammation and autoimmunity. It is thus important to mention that the HLA signal observed in the ConLiGen study 19 tags the HLA-DM gene which encodes molecules pivotal for the HLA class-II antigen processing pathway by facilitating antigenic peptide loading to HLA class-II alleles 34 .
More precisely, in the analysis of patients belonging to GWAS1, we observed that within the HLA-DRB1 heavy chain, Tyrosine or Leucine at position 37, Arginine at position 71 and Phenylalanine at position 67, amino-acids all belonging to the HLA-DRB1*11:01 classical allele, are associated to a better response to Li. Such a finding could be discussed either at amino-acid or allele level. Indeed, Arginine 71 has previously been demonstrated to modulate disease severity in inflammatory/autoimmune settings. For example, in the context of inflammatory polyarthritis, a study evaluating the relationship between amino-acids at position 11, 71 and 74 and clinically defined inflammation showed that Arginine 71 belongs to haplotypes conferring protection against severe inflammation. In this study, while Valine at position 11 is the primary risk factor for severe clinical inflammation, a Serine at this position confers protection against clinically-defined inflammation and disease progression 35 . Of interest, HLA-DRB1*11:01 bears the reported protective haplotype against inflammation namely, Serine 11, Arginine 71 and Leucine 74. In addition, another study evaluating whether clinical severity in rheumatoid arthritis (RA) correlated with HLA genetic heterogeneity, demonstrated that among patients having the prototypic at risk shared epitope, consisting of specific amino-acids combinations at positions 70 to 74 of the 3HVR, Arginine 71 characterized patients who were non-producers of the rheumatoid factor autoantibody (RF), a subset of patients usually characterized with less aggressive and destructive disease as compared to those who produce RF 36 .
Our meta-analysis showed that patients having Alanine or Leucine at position 74 were more represented in the group of responders. Of note, HLA-DRB1*11:01 allele which bears Leucine at position 74 was also likely associated with good response. Patients with Arginine or Glutamic acid at position 74 belong to the non-responder subgroup. These associations were also observed in the analysis of extreme groups where homozygous subjects for the A/L alleles were always found in smaller numbers in the extreme non-responders as compared to responders (Table 4, Supplementary Table 1), altogether suggesting a recessive effect. The beneficial effect exerted by the HLA-DRB1*11:01 allele can be discussed regarding its involvement in immune disease settings. A Brazilian study of the genetic influence on multiple sclerosis, a chronic autoimmune condition, showed that beside the canonical at risk HLA-DRB1*15 allele, the HLA-DRB1*11 variant exerts a strong protective effect against the disorder 37 . Such protective effect against multiple sclerosis was also reported in ethnically distant populationgroups 38 . Moreover, a recent study allowed us to establish a link between our finding and the signal identified in the primary ConLiGen GWAS study concerning the association with the HLA-DM region 19 . Indeed, the authors, analyzing the mechanism underlying the protective status conferred by HLA-DRB1*01:01 and HLA-DRB1*11:01 alleles against multiple sclerosis, demonstrated the involvement of the HLA-DM-mediated lower efficiency to bind antigenic peptide to the antigen peptide groove of the protective allele with consequent failure of peptide presentation at the cell surface and possible alleviation of inflammatory/autoimmune processes 39 . In another immune condition, namely chronic Hepatitis C Virus infection, the HLA-DRB1*11 was demonstrated to mediate protection against liver damage 40 . Altogether, these data suggest that good Li responders are associated with HLA variants conveying low inflammatory characteristics.
Concerning the non-responder patient group, an association was found with Glutamic acid or Arginine at position 74 in HLA-DRB1. The latter has been shown to drive the risk of developing Grave's disease, a common inflammatory/auto-immune disorder either in patients bearing the usually associated HLA-DRB1*03 specificity or in affected non-DR3 Graves patients 41 . In addition, we also observed an association between weak response to Li and multiple amino-acids scattered along the heavy chain of the HLA-DQB1 gene, almost all corresponding to the HLA-DQB1*02 allele. This association may suggest the involvement of the well-known HLA-8.1 ancestral haplotype (A*01 ~ B*08 ~ DRB1*03 ~ DQB1*02), which is characterized by a steady state pro-inflammatory status in healthy subjects and is the most associated HLA haplotype with inflammatory/autoimmune disorders. Another possibility explaining the association with HLA-DQB1*02 could be merely a LD with another HLA variant, not uncovered in the present study such as HLA-DRB1*03 known to be in LD with HLA-DQB1*02 and to bear the non-response associated Arginine 74. This hypothesis needs to be confirmed in larger cohorts of patients as it may have practical therapeutic implications such as add-on treatments of anti-inflammatory drugs in non-responders to Li.
Our study was performed on patients coming from multiple countries, which inevitably introduce heterogeneity in terms of data collection, methods of diagnosis/evaluation of the phenotype and information on patients such as comorbidities, overall constitutes possible confounding effects. Within this line, the fact that the distribution of Alda scores in the GWAS1 group was slightly more variable than that observed in the GWAS2 group and that patients belonging to the GWAS1 were slightly older than that of the GWAS2 group (Table 1), may also contribute to differences in the analysis of the two groups. In spite of these limits, the association between the HLA-DRB1 position 74 A/L allele and responders was consistently observed under the recessive model with extreme phenotypes in all the countries tested in the GWAS1 and in the GWAS2. We have already mentioned www.nature.com/scientificreports/ that the recessive effect for HLA-DRB1 cannot be dissociated from a dominant effect for the HLA-DQB1*02 allele in terms of genetic association. Such recessive-based observation may reflect the characteristics of the HLA system (and other innate immune sensors) in terms of structure, function and evolution where the heterozygosity advantage was selected to better fight pathogen diversity, a recessive status may prevent overwhelmed inflammatory processes [42][43][44] . Altogether, beside the well-documented anti-inflammatory effect exerted by Li on the innate GSK-3β pathways and consequently on circulating pro-and anti-inflammatory cytokines, our observations indicate that response to Li in BD patients may be underpinned by adaptive HLA-mediated inflammatory processes, fitting with the known anti-inflammatory effect driven by Li that may contribute to its therapeutic efficiency 31 . Our genetic findings showing an association between HLA-mediated lower inflammation in patients with efficient response to Li, may suggest that non-responders among BD patients could have inflammatory backgrounds that override Li mediated anti-inflammatory properties.

Materials and methods
Study population. This study includes GWAS genotype and phenotype data of Li-treated BD patients from 22 participating sites belonging to the International Consortium on Lithium Genetics (ConLiGen) 16 . The whole cohort consisted of 2588 patients resulting from two GWAS of 1187 and 1401 individuals named GWAS1 and GWAS2 in the original publication 16 . Samples were collected and genotyped in two distinct steps, determining the GWAS1 and GWAS2 groups [detailed rationale is available in the appendix of Hou et al. 16 ]. The various cohorts were initially developed separately and various chips were used for genotyping 16 . We kept the same distribution of the samples and the names GWAS1 and GWAS2. After removal of Japanese, Taiwanese and Sardinian patients given the poor applicability of HLA imputation in these population-groups on the one hand, and the difference in HLA frequencies among the Sardinian population on the other, 862 and 1309 patients belonging respectively to the GWAS1 and GWAS2 were subjected to analysis. We performed a quality control of the genotypes for each patient. As cryptic relatedness between individuals can introduce a bias in association analysis and in the following meta-analysis, we computed identity-by-descent (IBD) between all pairs of individuals for the two groups and removed one individual of the related pair if a second degree related can be detected (N = 14 for GWAS1 and N = 11 for the GWAS2). We also removed one patient in the GWAS1 and six patients in the GWAS2 because of an excess of heterozygosity, suggesting a poor DNA quality. After these two steps, 847 patients from GWAS1 and 1292 from GWAS2 remained for analysis.
The Conligen consortium used data obtained through an international collaboration. The University of Heidelberg Ethics Committee provided the central ethics approval for the consortium. Written informed consent was obtained from each patient (they were all aged > 18) according to the study protocols of the participating cohorts and their respective institutions. All procedures were performed in accordance with the guidelines of the institutional research committees and of the Helsinki Declaration.
Study cohort phenotypes. From the discovery GWAS participants 16 , long-term treatment responses to Li were evaluated using the Alda scale. Briefly, the Alda scale used criterion A to rate the degree of response (activity of illness while on adequate Li treatment) on a 10-point scale and criterion B which establishes the relationship between improvement and treatment (i.e. number of episodes off treatment, frequency of episode off treatment, duration of treatment, compliance, use of additional treatment). The total score is obtained by subtracting score B from score A and negative scores are set to 0, so that the total score ranges from 0 to 10.
In the present study, we performed 3 types of phenotype analyses in order to dissect the genetic associations with Li response: (1) A continuous phenotype (range 0-10) using the A score, excluding all individuals with a total B score greater than 4 for a dimensional assessment of Li response. In this study, there were 772 patients studied in GWAS1 and 1080 patients studied in GWAS2 by linear regression with the continuous phenotype in the additive model.
(2) A dichotomous phenotype comparing good vs non responders to Li, defined by an A Alda score ≥ 7 for responders and ≤ 3 for non-responders. In this comparative study (in the additive model), there were 186 nonresponders compared with 418 responders from GWAS1, and 201 non-responders compared with 685 responders from GWAS2. To take into account the variations of distribution of HLA alleles according to countries, we also checked country by country the frequency of patients carrying the HLA-DRB1 A or L alleles at position 74 in the non-responder and responder groups, when these groups had enough patients (N ≥ 15) to compute a reliable frequency. We added the Sardinian population (comprising 52 non-responders compared with 91 good responders, Table 4), even if it was excluded from the statistical analysis.
(3) We also performed a country-by-country analysis in extreme phenotypes defined by an A Alda score ≥ 9 for responders and ≤ 3 for non-responders: there were 186 non-responders compared with 275 responders from GWAS1, and 201 non-responders compared with 362 responders from GWAS2. In Sardinian population, there were 52 non-responders compared with 39 good responders (Supplementary Table 1). We also computed a FDR value comparing the non-responders with responders in the recessive model with this extreme phenotype (Supplementary Table 5).
Imputation of HLA variants. The SNP2HLA software 45 has been developed to impute directly HLA class-I and class II-classical alleles at the 4-digit level from the chromosome 6 genotypic data of each patient. In addition, the software computes directly the amino-acid variants from the imputed 4-digit HLA alleles, allowing to test associations also at the level of the HLA protein amino-acid variants (including the ones constituting the peptide antigen presenting groove) thus facilitating the biological interpretation. A total of 1708 imputed vari- www.nature.com/scientificreports/ ants (433 classical HLA alleles and 1275 amino-acid variations) were generated. The variants with minor allele frequency < 0.05 and/or p-value for Hardy-Weinberg departure test < 1 × 10 −6 were removed. The resulting 85 classical HLA alleles and 880 amino-acid variations were subjected to the statistical analysis pipeline. To validate the accuracy of the imputation, we compared the frequency of the variants with those listed in the Allele frequency net database 25 . Given that the studied populations were from European descent, we choose available representative populations from France (France, Rennes, N = 200; France, west, N = 100) and from Germany (Germany, Essen, N = 174) in the Allele Frequency Net Database with alleles given in at least four digits. The best variants obtained from our association analysis had a rather strong MAF, and since the MAF found in the control populations were similar it was sufficient to confirm the quality of the HLA imputation for thee variants, even if these control populations had a relatively modest size. For the amino-acid frequencies we used the tool of translation from Allele frequency net database 25 .
Statistical analysis. For the 965 HLA imputed variants having passed the quality controls, associations between HLA amino-acids or alleles and response to Li were tested. The statistical analysis was performed by a multivariate linear or logistic regression using PLINK software 46 under additive model for response to Li (continuous or binary phenotypes). To correct for possible population stratification, we performed principal component analysis (PCA) on the genotypes and used the top two principal components as covariates in statistical analyses. To avoid spurious association, we added age and gender of patients as covariates. Finally, a classical meta-analysis was performed between associations resulting from GWAS1 and GWAS2 for two phenotypes (continuous (1) and dichotomous (2)) using PLINK. We used the Bonferroni threshold (for 965 tests − 85 classical HLA alleles + 880 amino-acid-, p-value < 5.2 × 10 −5 ) to determine the significance of the associations and we also took the false discovery rate (FDR) < 0.25 for the suggestive associations. www.nature.com/scientificreports/