A phenome-wide association study of 26 mendelian genes reveals phenotypic expressivity of common and rare variants within the general population

The clinical evaluation of a genetic syndrome relies upon recognition of a characteristic pattern of signs or symptoms to guide targeted genetic testing for confirmation of the diagnosis. However, individuals displaying a single phenotype of a complex syndrome may not meet criteria for clinical diagnosis or genetic testing. Here, we present a phenome-wide association study (PheWAS) approach to systematically explore the phenotypic expressivity of common and rare alleles in genes associated with four well-described syndromic diseases (Alagille (AS), Marfan (MS), DiGeorge (DS), and Noonan (NS) syndromes) in the general population. Using human phenotype ontology (HPO) terms, we systematically mapped 60 phenotypes related to AS, MS, DS and NS in 337,198 unrelated white British from the UK Biobank (UKBB) based on their hospital admission records, self-administrated questionnaires, and physiological measurements. We performed logistic regression adjusting for age, sex, and the first 5 genetic principal components, for each phenotype and each variant in the target genes (JAG1, NOTCH2 FBN1, PTPN1 and RAS-opathy genes, and genes in the 22q11.2 locus) and performed a gene burden test. Overall, we observed multiple phenotype-genotype correlations, such as the association between variation in JAG1, FBN1, PTPN11 and SOS2 with diastolic and systolic blood pressure; and pleiotropy among multiple variants in syndromic genes. For example, rs11066309 in PTPN11 was significantly associated with a lower body mass index, an increased risk of hypothyroidism and a smaller size for gestational age, all in concordance with NS-related phenotypes. Similarly, rs589668 in FBN1 was associated with an increase in body height and blood pressure, and a reduced body fat percentage as observed in Marfan syndrome. Our findings suggest that the spectrum of associations of common and rare variants in genes involved in syndromic diseases can be extended to individual phenotypes within the general population.


Introduction
Genetic syndromes are rare diseases defined by a specific and clinically recognizable set of phenotypes across multiple organ systems. The era of next-generation sequencing has enabled substantial progress in linking syndromic disease to specific genetic loci, coupled with public databases of genotype-phenotype relationships to facilitate the classification of genetic variants from "benign" to "pathogenic" for use in clinical decision making. Large population-scale databases of genetic variation without phenotypes, such as ExAC, have provided additional context for characterizing genotype-phenotype relationships in genetic disease [1]. For mutations previously thought to cause disease, population databases have often suggested lower estimates of penetrance than initially recognized [2,3].
The diagnosis or classification of an individual with genetic syndrome relies upon expert recognition of a characteristic pattern of signs or symptoms or a set of defined diagnostic criteria. However, individuals displaying single phenotypes of a complex syndrome may not meet criteria for clinical diagnosis or genetic testing; expanding a binary definition of syndromic phenotypes to phenotype scores can identify more individuals with Mendelian disease patterns [4]. Similarly, individuals with clearly pathogenic mutations may be affected with only a single component phenotype of a genetic syndrome [5,6]. However, the descriptions of allelic heterogeneity, penetrance, and expressivity in syndromic disease genes have focused almost exclusively upon rare or familial alleles [7,8] Recent studies have shown that rare and common variants in or near mendelian diseases genes are associated with complex traits in the general population [9][10][11]. Moreover, Freund et al [11] demonstrated an enrichment of signal from the summary statistics of Genome Wide Association Studies (GWAS) near syndromic disease genes for the related phenotypes. However, this work was based on the curation of available GWAS summary statistics.
Here, we present a phenome-wide association study (PheWAS) approach to systematically explore expressivity of common and rare alleles in genes associated with four well-described syndromic diseases in the general population. Using the UK Biobank, we linked individuallevel medical and morphometric data to the characteristic phenotypes of Alagille (AS), Marfan (MS), DiGeorge (DS), and Noonan (NS) syndromes. These data allow a survey of the association of common and rare alleles to single component phenotypes of each syndrome within the general (non-syndromic) population.

Results
Based on the Human Phenotype Ontology (HPO)-an ontology-based system developed using medical literature and other ontology-based systems [12]-we identified 196 HPO terms related to AS, MF, DS, and NS. Of these 196 HPO terms, 53 were shared between at least two syndromes, and seven terms were included in all four syndromes (S1 Table). After grouping the HPO terms into categories based on affected organ systems, there were 115 HPO terms of which 73 could be matched to 100 phenotypes available in the UKBB. We additionally included liver and renal serum biomarkers such as alanine aminotransferase, creatinine, and direct bilirubin to capture liver and renal dysfunction observed in some of the genetic syndrome. Most of the unmatched phenotypes were related to specific abnormalities of body structure or the musculoskeletal system, which were poorly represented in clinical and billing codes, or measurements such as impaired T-cell function, not available in the UKBB.

Genotype-phenotype associations are common across all syndromic genes
We tested the associations between all variants and all phenotypes included in our study. A total of 1,824,564 tests (84 phenotypes x 21,721 variants) were performed. Overall, we found significant association between 20 phenotypes, and multiple variants across JAG1, FBN1, PTPN11, SOS2, RIT1, RAF1, KAT6B, RASA2, MAP2K1, CBL, DGCR2 and COMT (Fig 1A and  S3 Table). Using stepwise conditional analysis implemented in GCTA, we identified a subset of 46 variants independently associated with those phenotypes (Fig 1B and S4 Table); among which 9 variants were associated with more 2 or more phenotypes (Fig 1B and Table 1). Among the phenotypes with significant associations, hypothyroidism (HP0000821), diastolic BP (HP0005117), systolic BP (HP0004421), standing/sitting height ratio (abnormality of body height; HP0000002), birth weight (small for gestational age; HP0001518), amount of subcutaneous adipose tissues or body fat percent (reduced subcutaneous adipose tissue; HP0003758), growth abnormality (HP0001507), body mass index (abnormality of body mass; HP0045081), hyperlipidemia (HP0003124), direct bilirubin level, creatinine level, aspartate amino transferase level, and alkaline phosphate phosphatase level were significantly associated with variants across multiples genes (Fig 1A and Fig 1B).
The set of independent variants associated with multiple phenotypes are represented in Table 1. Diastolic BP and systolic BP along with body mass index displayed a genetic association in each of the four syndromes. Birth weight, subcutaneous adipose tissue (body fat percent), and tall stature for MS or short stature for NS are common phenotypes; while hypothyroidism, and growth retardation are reported in both NS and DS. In order to replicate our findings, we looked in the GWAS catalogue, for sets of variants in genes in association with phenotypes or proxy-phenotypes included in our study. We excluded studies performed . The color indicates variant in each gene and the shape indicates each phenotype (ex: variants in PTPN11 are represented in blue, and the association with diastolic blood pressure indicated with the sign +). Correlation plot of association between phenotypes and the subset of independent variants within each gene is displayed in panel (b). The points represent the z-score and the color represent the direction of the association. The color varies from purple (inverse association) to red (positive association). The size of the point corresponds to the p-value (-log10(p)); the stars indicate variants associated with multiple phenotypes.
https://doi.org/10.1371/journal.pgen.1008802.g001 in non-Europeans or in the UK Biobank. We then extracted for each remaining phenotype, the association between candidate variants (defined variants significantly associated with the same phenotype in our study or in high LD (R2>0.8) with the associated variant) and the corresponding phenotype. The association between variants in FBN1, PTPN11, SOS2, and JAG1 and, diastolic and systolic BP, were reported in the GWAS catalogue with the direction of effect concordant with the observed effect in our study. Similarly, the association between variants in PTPN11 and hypothyroidism as well as the association between JAG1 and birth weight were reported (S5 Table). At a gene level, when using SKAT combined with weighted CADD score, 24 phenotypes including standing/sitting height ratio, blood pressure (systolic and diastolic), amount of subcutaneous fat, hypothyroidism, and hypercholesterolemia were significantly associated with several genes after multiple testing correction (p.fdr<0.05) ( Table 2, Fig 2A and Fig 2B).

Variation in syndromic genes are associated with component phenotypes
Marfan syndrome (MS) is a primary disorder of connective tissue with diagnostic criteria centered around cardiovascular, musculoskeletal, and ocular phenotypes linked to a single gene FBN1 which encodes an extracellular matrix protein. Several variants in FBN1 were significantly associated with increased standing/sitting height ratio and an elevated diastolic BP. An increased risk of aortic dissection and a lower percent of body fat (two major phenotypes in MS) were observed for several of these variants although the association was merely suggestive (Fig 3A and S3 Table). All variants in FBN1 displaying associations were located within the same LD block and were highly correlated with each other (Fig 3A and 3B). Using conditional regression analysis implemented in gcta on each phenotype, one independent signal was identified and rs589668 was tested with multiple phenotypes in FBN1. The variant rs589668 displays the top signal with high standing/sitting height ratio (p = 8x10 -28 , Table 1 and Fig 1B), an elevated diastolic BP (beta = 0.02, se = 0.002, p = 8x10 -09 ), and a lower percent of body fat (beta = -0.08, se = 0,02, p = 5x10 -05 , Table 1). The association observed with these phenotypes were as expected, in the same direction of effect as observed in Marfan syndrome. For instance, individuals with Marfan syndrome often have an elevated standing and sitting height ratio, and thin skin due to very small amounts of subcutaneous fat. At the gene level association using the SKAT test with variants weighted by CADD score (Combined Annotation Depletion Dependent), standing/sitting height ratio, systolic and diastolic BP, subcutaneous adipose tissue, and aortic dissection were significantly associated with FBN1 ( Fig 2B and Table 2). For NS and RAS-opathy related phenotypes, variants in PTPN11 were associated with increased risk of hypothyroidism, high diastolic and systolic BP, and high standing/sitting height ratio (Fig 1A, Fig 1B and S3 Table). variants in SOS2 were associated with lower systolic and diastolic BP, and lower percent of body fat (Fig 1A and S3 Table). variants in MAP2K1, and KAT8B were associated with lower body mass index, as well as lower level of cutaneous adipocyte tissues (Fig 1A, Fig 1B and S3 Table). variants in RASA2 were associated with body mass index, level of subcutaneous adipose tissues, a lower ratio of standing/sitting height, and growth abnormality. We also found a significant association between variants in SHOC2 and high level of Cystatin C, Creatinine, and Urea which reflect kidney dysfunction (Fig 1A, Fig 1B   Fig 2. Primary PheWAS results at a gene level. Plot of PheWAS results for all genes and phenotypes (a). The red line represents the level of significance after FDR correction. Genes are represented by color and phenotypes are indicated by shape. Correlation plot for the set of significant gene-phenotype pairs. (b). The color represents the p-value of association and varies from none associated (grey-white) to significant association (red-dark red). https://doi.org/10.1371/journal.pgen.1008802.g002

PLOS GENETICS
and S3 Table). In addition, direct bilirubin concentration was significantly associated with variants in SOS1, RIT1, CBL and SHOC2 (Fig 1A and S3 Table). variants in PTPN11 display a moderate to low correlation, while high correlation was observed between variants in RASA2, SOS2 and MAP2K1 indicating that the association observed within each gene represents a single signal (Fig 4B). To identify independent signals within each gene for each associated phenotype, we performed conditional regression using stepwise selection procedure implemented in GCTA. For each subset of associated SNPs-phenotypes pairs within each gene, we identify one independent signal (Fig 1B). Among SNPs associated with multiple phenotypes, rs11066309 in PTPN11 displays a strong association with increased risk for hypothyroidism (ALT freq = 0.40; OR [95% CI]: 1.19; [1.16-1.21]; p = 6x10 -59 ) along with five other phenotypes, including decreased body mass index (beta = -0.012, p = 1.13x10 -06 ) and birth weight (beta = -0.020, p = 2.95x10 -10 ) (Fig 1B and Table 1).
DiGeorge syndrome encompasses a recurrent microdeletion of multiple genes at the 22q11.2 locus due to the presence of segmental duplications, with affected individuals displaying neuropsychiatric, immunological, and cardiovascular phenotypes originating from defects in neural crest cell formation and migration. At the variant level, rs807747 in DGCR2 and rs71313931 in COMT were significantly associated with abnormal body height (p = 2x10 -07 ) and systolic BP (p = 5.6x10 -17 ), respectively (Fig 1A and 1B) while rs72646939 was associated with elevated creatinine levels (Fig 1A and 1B). At the gene level, COMT was associated with systolic BP, body mass index, and hypercholesterolemia while DGCR2 was associated with attention deficit disorder, urine level concentration, and DGCR8 associated with body mass index, systolic blood pressure, and standing/sitting height ratio (p.fdr<0.05) ( Table 2).

Evaluation of pleiotropy and epistasis between phenotypes and genotypes
After identifying multiple pairs of variant-phenotype associations among which, 43 SNPs were independent; 9 of the 43 SNPs were associated with two or more phenotypes (Fig 1B). To access whether association between a single variant and multiple phenotypes are independent or due to correlation between phenotypes, we performed a formal test of pleiotropy between each independent variant and the associated phenotypes. We found a significant association between all pleiotropic SNPs and the associated phenotypes (S7 Table). Similarly, to evaluate whether associations between a phenotype and multiple variants were explained by variantvariant interaction or epistasis we performed a stepwise linear or logistic regression with interaction terms between pairs of associated variant for each phenotype. After multiple testing  Table).

Discussion
Here, we systematically describe the association of variation in 26 mendelian genes linked to four syndromic diseases with the component phenotypes of the corresponding syndromic disease. We hypothesize that in the general population, common and rare alleles for syndromic diseases display pleiotropic effects with the phenotypes related to genetic syndromes. Using the UKBB, we linked individual-level data to the characteristic phenotypes of Alagille, Marfan, Noonan, and DiGeorge syndromes, showing clearly the association of common and rare alleles to single component phenotypes of each syndrome. Individual phenotypes that are modulated by variants within the loci syndromic disease genes appear to be present within the general population. These findings are consistent with the data reported by Bastarache et al [4] which suggested that scaling of component symptoms of rare disease into a continuous phenotyping score can improve the identification of individuals with rare diseases.
Within families of individuals affected by syndromic disease carrying the same pathogenic mutation, the expressivity of component phenotypes may vary in different individuals [13,14].
Here, we show that many common and rare variants within loci of syndromic disease genes existing in the general population may result in expression of traits and phenotypes closely related to the syndrome of interest. For example, we observe associations of a common intronic variant in FBN1 rs589668 (MAF = 0.25 in European populations) with increases in blood pressure and height and decreased subcutaneous fat distribution. In GTEx [15], this variant is an eQTL strongly associated with decreased expression in whole blood (p = 1.7x10 -37 ), which would be concordant with the known molecular mechanism of FBN1 pathogenesis in MS: pathogenic alleles impairing gene function result in increased height and abnormal fat distribution and increased arterial stiffness [16,17].
Modifiers of penetrance and phenotypic expressivity in Marfan syndrome have been previously proposed [18,19] and a study based on a single Marfan syndrome family also suggested that differences in normal FBN1 expression could contribute to the clinical variability of Marfan syndrome [19]. We observed four additional variants (rs11070641, rs4775760, rs363832 and rs140605) in FBN1 to be associated with high standing/sitting height ratio, a characteristic feature often observed in Marfan syndrome. These variants are reported as benign variants in CLINVAR suggesting that they do not cause primary syndromic disease, but our data suggest they may be modifiers of penetrance for the phenotype of height. Our results suggest that common variants and local haplotype structure around syndromic genes may deserve more attention [20].
Noonan syndrome is caused by mutations in PTPN11 and part of a group of related disorders arising from activating mutations in RAS-MAPK signaling pathway known as RASopathy which display many phenotypes across a variety of organ systems. A wide phenotypic variability and genetic heterogeneity have also been described in individuals with Noonan syndrome in relation to rare variants in PTPN11 [22]. Here, we show that even in the general population, common and rare variants in PTPN11 are independently associated with phenotypes such as hypothyroidism, small birth weight and low percent of body fat observed in some cases of Noonan syndrome [21][22][23]. In GTEx [15], numerous variants in PTPN11, such as rs11066309, rs3741983 and rs11066322 were significantly associated with a decreased expression in atrial appendage, adipose tissue, thyroid and skin and esophagus. Although in consistent with the role of PTPN11 in thyroid function, cancer and autoimmunity [24][25][26], these variants are instead described as eQTLs with TMEM116, ALDH2 and MAPKAPK5-AS1 located up to 500kb upstream of PTPN11, suggesting that the associations observed with rs11066309, rs3741983 and rs11066322 could potentially also arise from associations with these other genes.
Growth retardation, lower BMI and short stature are additional well-known characteristics of Noonan syndrome and display a phenotype-genotype variability of growth patterns in affected individuals [27]. In concordance with this study, we showed that in the general population, common and rare variants in RASA2, SOS2 and MAP2K1 are independently associated with growth characteristics (body mass index, height and growth abnormality) and the association driven by one or more haplotypes in each gene. Among the RAS-MAPK signaling pathway genes tested, we observe a significant association between the related phenotypes (lower body mass index, growth retardation, low percent of body fat) and rs3741983 (PTPN11), rs72681869 (SOS2), rs61755579 (SOS2), and rs112542693 (MAP2K2) reported in CLINVAR as "benign". Although these variants are indeed not sufficient to cause mendelian disease, they may nonetheless contribute to specific phenotypes related to Noonan syndrome when a "pathogenic" variant is present.
When performing genetic testing, allele frequency is often incorporated into an assessment of the pathogenicity of a genetic variant. Common variation in and around JAG1 has previously been associated with such disparate phenotypes as pulse pressure, circulating blood indices, and birthweight, and none of the variants included in our analysis appeared to be directly associated with the component phenotypes of AS. However, the unifying molecular abnormality in AS are defects in vascular formation which lead to each of the component cardiovascular and liver phenotypes of Alagille syndrome [28,29]. The pleiotropic effect detected for common alleles in JAG1 (S7 Table) with multiple different phenotypes, may be linked to the underlying role in vascular formation.
Our study has some limitations. Our analysis was limited to phenotypes with more than 100 cases, and variants with minor allele frequency of at least 0.0001. Therefore, diseases with relatively rare prevalence or variants with extremely rare frequencies were not analyzed. In addition, because our study cohort consist of adults from the general population, specific phenotypes targeting facial and skeletal dysmorphism, such as butterfly vertebrae or broad forehead; specific abnormalities of organs, such as biliary disease were not present. However, to work around the absence of some phenotypes, we used proxy phenotypes or measurements present in the UKBB, such as head circumference as an alternative for broad forehead, education level for ADHD, weight and height at age 10 as proxy for growth abnormality. All things considered, the complexity of matching UKBB phenotypes to HPO terms may simply not capture some phenotypes, despite manual curation. An additional limitation of our study is the fact that, it is not possible to "diagnose" individuals from the data available in the UK Biobank in order to exclude them from analysis.
Key strengths of our study include the ability to systematically test multiple phenotypegenotype association and to highlight phenotypic expressivity of different variants linked to syndromic genes. Our study maps UKBB phenotypes to HPO terms and shows that common and rare variants in genes responsible for Alagille, Marfan, Noonan, and DiGeorge syndromes, are also independently associated with component phenotypes of these syndromes in the general population.
Our findings suggest that within the general population both common and rare variation in syndromic disease genes may be associated with component phenotypes of a syndrome. Further research on the expressivity of alleles in genes in the general population is needed to link our understanding of Mendelian syndromes with complex trait genetics.

Study population and data collection
The study cohort was derived from the UK Biobank (UKBB), a large prospective cohort study with comprehensive health data from over 500,000 volunteer participants in the United Kingdom aged 37-73 years at recruitment in 2006-2010. The cohort has previously been described in detail [30][31][32]. Information on the UK biobank participants was collected at enrollment, and from electronic health record (EHR) information which includes diagnostic codes (ICD10, ICD9) and procedural codes (OPCS) from hospital admission records dating to 1992, and cancer registries. Data collected at the assessment visit included information on a participant's health and lifestyle, hearing and cognitive function, collected through a touchscreen questionnaire and verbal interview. A range of physical measurements was also performed, including blood pressure; arterial stiffness; body composition measures (including impedance); hand-grip strength; ultrasound bone densitometry; spirometry; and an exercise/fitness test with ECG. Samples of blood, urine, and saliva were also collected. Medical phenotypes were aggregated as previously described, incorporating available information including a broad set of medical phenotypes defined using computational matching and manual curation of on hospital in-patient record data (ICD10 and ICD9 codes), self-reported verbal questionnaire data, and cancer and death registry data [33,34].

Phenotypes of target syndromes
We identified phenotypes related to syndromic diseases through the Human Phenotype Ontology (HPO). HPO is an ontology-based system developed using medical literature, and other ontology-based systems such as Orphanet, and OMIM [12]. HPO provides a standardized vocabulary of phenotypic and abnormalities encountered in human diseases. The HPO has link symptoms/phenotypes to diseases or genetic disorders, and the causing genes. As an example, Alagille syndrome (AS) is linked to JAG1, and NOTCH2 genes as well as all the phenotypes or symptoms observed in AS, such as atrial septal defect, hypertelorism, and butterfly vertebra.
HPO terms were directly matched to UKBB phenotypes when phenotypes in both systems had similar terminology. The direct phenotype matching was conducted using a semi-automatic mapping system which combines semantic and lexical similarity between word [35] followed by manual curation. When the HPO terms were not present, we performed an indirect matching by hand to find in the UKBB, the phenotype that best reflects the target HPO terms. For example, abnormality of body structure or body morphology such as abnormal body height, reduced subcutaneous adipose tissues, bone density or broad forehead, were respectively matched to sitting/standing height ratio; body fat percentage; bone mineral density, and head bone area. blood biomarkers measuring liver, and kidney functions such as direct bilirubin, creatinine, Alanine aminotransferase, Alkaline phosphatase, Gamma glutamyl-transferase were used as proxy for liver, or renal function. For psychiatric diseases such as depression and neurodevelopmental disorders such as attention deficit and hyperactivity disorder (ADHD), we used a score of depressive symptoms and self-reported educational level respectively as proxies for these terms.
To increase the number of subjects in some subgroup of phenotypes, we combined subcategories of HPO terms into a group or category. For example, 39 HPO terms representing an abnormality of head, ears, and eyes such as low-set ears, strabismus, macrotia, webbed neck, short neck, abnormality of the eye, microcornea, down-slanted palpebral fissure and other congenital abnormality of ears, were grouped into "Abnormality of head or neck (HP0000152)" and mapped to icd10 targeting congenital malformations of eye, ear, face, and neck and other organs especially facial appearance (ICD10: Q10 to Q18 and Q87). Ten HPO terms for congenital abnormality of cardiovascular system including Ventricular septal defect, Atrial septal defect, Tetralogy of Fallot, Patent ductus arteriosus, Bicuspid aortic valve, Truncus arteriosus, Coarctation of aorta, Tricuspid valve prolapse were combined into abnormality of the cardiovascular system (HP0001626) and mapped to Congenital malformations of the circulatory system (ICD10: Q20 to Q28).

Genotyping data
Genotyping was performed using the Affymetrix UK BiLEVE Axiom array on an initial 50,000 participants; the remaining 450,000 participants were genotyped using the Affymetrix UK Biobank Axiom1 array. The two arrays are extremely similar (with over 95% common content). Quality control and imputation to over 90 million variants, indels and large structural variants was performed [35].

Gene definitions
Using OMIM, and HPO, we identified 26 genes linked to the syndromes of our interest (Table 3). Each gene was defined from 5'UTR to 3'UTR with an extra additional 5kb upstream, and 5kb downstream the gene. To account for variants in regulatory elements of the target gene but located outside of the defined boundary, we additionally include within each target gene, variants located outside of the defined boundary but in eQTL in any tissue with the target gene. For the variant level association, we further extend the boundary of each gene to 50 kb upstream, and 50 kb downstream a gene. We identify a total of 21,712 variants in 26 genes related to Alagille syndrome, Marfan syndrome, Noonan syndrome, and DiGeorge syndrome were selected for our study ( Table 2). The selected variants had a MAF � 0.0001 and an imputation measurement (R2) � 0.6 Statistical analysis SNP level. We performed the association between all 84 identified phenotypes, and all 21,721 variants in all the syndromic diseases genes included in our analysis. For binary traits, logistic regression with adjustment on age, sex, batch, and the top 5 principal components were used. First, regression was used in a situation of unbalanced numbers of cases and controls, especially when the number of cases was very small (less than 200 cases). For continuous traits, we performed linear regression with adjustment on age, sex, batch, and the top 5 principal components. Our analysis was restricted to individuals of European descent, due to the relatively small number of individuals from other ethnic groups in the UKBB. Bonferroni correction based on the number of independent tests was used to correct on multiple testing. Given the high correlation between variants within gene or regions, Bonferroni correction is often stringent when the number of tests considered is number of SNPs time the number of phenotypes. To take in account the correlation between variants, we estimate the number of independent variants in a block of 50 kb with a correlation > 0.8 using the pairwise pruning method implemented in PLINK which estimated 2166 independent variants within our target regions. We apply a threshold of 2.7x10 -07 = 0.05/(2166x84) independent tests. In order to replicate our finding, we looked in the GWAS catalog [36], for sets of variants in genes in association with phenotypes or proxy-phenotypes included in our study. We excluded studies performed in non-Europeans or in the UK Biobank. We then extracted for each remaining phenotype, the association between candidate variants (defined variants significantly associated with the same phenotype in our study or in high LD (R2>0.8) with the associated variant) and the corresponding phenotype Identification of an independent variant-trait association set. To identify a subset of variants independently associated with each phenotype, we performed the stepwise model selection for identification of variants independently associated with a phenotype; implemented in GCTA [37].
Pleiotropy and epistasis assessment. To assess whether variants associated with multiple phenotypes reflect a correlation between phenotypes or are independently associated with each phenotype, we performed a pleiotropy test between each variant and the set of associated phenotypes using pleio, a R package for pleiotropy assessment. The association between each variant and a group of phenotypes were considered significant if the p-value were less than 2.1x10 -04 (p<0.05/233 significantly associated variants). Similarly, to test for interaction between variants within a single gene or across different genes, we performed an epistasis test which consist of testing for the interaction between each pairs of associated variants and the corresponding phenotypes. We used linear regression for continuous phenotype and logistic regression for binary phenotypes with adjustment for age, sex and the first 10 PCs. An interaction term was considered significant if the pvalue were less than 3.7x10 -05 based on Bonferroni correction (0.05/1320 tests).
Gene level. At a gene level, we performed a Sequence Kernel association test (SKAT) using a sequence kernel method as well as a burden test [38,39]. We performed the SKAT test on rare and common variants as well as on rare variants only. To account for the contribution of rare variants and common variants, we use SKAT CommonRare methods in which, rare and common variants are partitioned into two groups to test for the association with the phenotypes; the results of association is then combined using combined multivariate collapsing [38]. A variant was considered rare if the Minor allele frequency was less or equal to 0.05 (MAF �0.05). To account for their possible functional relevance, each variant was weighted in the SKAT test by their CADD score (Combined Annotation Depletion Dependent) [39,40]. Although Gene-based SKAT tests are relatively insensitive, for sensitivity analysis, we also performed SKAT using allelic frequency. Each gene was defined from 5'UTR to 3'UTR with an extra 5kb upstream, and 5kb downstream. To account for variants in regulatory elements of the target genes but located outside of the defined boundary, we additionally include within each target gene, variants in eQTL with the target gene in any tissue but located outside of the defined boundary. We used FDR to correct for multiple testing.
Supporting information S1  Table. Gene level association with different SKAT test P.weight CADD (rare variant test each variants are weight by their CADD score); P.RaCo (SKAT with adaptive sum test of rare and common variants); P.Burden (SKAT burden test with rare and common variants aggregation); SKAT rare (rare variant test only). (XLSX) S7