Background: The FOXC2 gene on 16q24 is mutated in lymphoedema distichiasis (LD), in which varicose veins (VV) are a common feature. We hypothesised that this gene might be implicated in the development of VV in the normal population, therefore, after performing a classical twin study, we tested for linkage and association in white women. We also tested for linkage with haemorrhoids (H), as a separate venous anomaly at the same locus.
Methods: A total of 2060 complete female twin pairs aged 18–80 years from the St Thomas’ Adult UK Twin registry replied to questions on VV and H as part of a broader postal survey of 6600 twins (62% response rate). Dizygotic female twin pairs were tested for linkage and association to the candidate marker D16S520 (1903 individuals genotyped), which is located about 80 kb from FOXC2.
Results: Casewise concordance rates were significantly higher for monozygotic than dizygotic twins for both phenotypes (VV 67% v 45%; p = 2.2×10−6; H 68% v 59%; p = 0.01; H including during pregnancy 73% v 64%; p = 2.1×10−4), corresponding to additive genetic heritabilities in liability of 86% (95% confidence interval (CI) 73% to 99%) for VV and 56–61% for H (95% CI 43% to 73%). The presence of VV and H were significantly correlated. We found significant evidence of linkage to the marker for VV (MLSASP = 1.37, p = 0.01; GLMASP/DSP Z = 3.17 p = 0.002), but no association. Both linkage and association tests were negative for H. The combined phenotype of having VV and H did not show any evidence of linkage or association.
Conclusion: These results demonstrate VV and H to be heritable, related conditions, and the data strongly suggest FOXC2 to be implicated in the development of VV in the general population.
- ASP, affected sibling pairs
- DSP, discordant sibling pairs
- DZ, dizygotic
- H, haemorrhoids
- LD, lymphoedema distichiasis
- MZ, monozygotic
- PIC, polymorphism information content
- VV, varicose veins
- lymphoedema distichiasis
- varicose veins
Statistics from Altmetric.com
- ASP, affected sibling pairs
- DSP, discordant sibling pairs
- DZ, dizygotic
- H, haemorrhoids
- LD, lymphoedema distichiasis
- MZ, monozygotic
- PIC, polymorphism information content
- VV, varicose veins
Varicose veins (VV) are probably the most common vascular abnormality. While there is known to be an environmental influence on the development of VV, there have been reports of possible familial clustering and estimates of heritability based on unusual families. Hauge and Gundersen1 suggested multifactorial inheritance in a very large pedigree, while Matousek and Prerovsky2 inferred a multifactorial heritability in liability of about 50%. No genes have so far been directly implicated as candidates for the development of VV, and no large scale studies have been performed to obtain an accurate estimate of heritability in the general population.
Lymphoedema distichiasis (LD) is a form of lymphoedema with associated features, including distichiasis (extra eye lashes from the meibomian glands), ptosis, congenital heart disease, and cleft palate.3 The gene was localised to markers D16S520 and D16S3037 in 1999,4 where venous insufficiency was noted in the families used in the analysis, a finding also noted by Rosbotham et al.5 The gene itself was identified as FOXC2 by Fang et al,6 and several mutation studies have shown that the vast majority of mutations are small insertions or deletions, almost certainly exerting an effect via haploinsufficiency.7–9 FOXC2 is a transcription factor that is involved in development, and it is assumed that the end organ effects seen in LD patients are present from birth, even though all the clinical effects are not. The lymphoedema, for example, is almost always of pre-pubertal or pubertal onset. Although studies had reported on families with LD by 2001,8,9 the full extent of venous abnormality in the disorder was not realised until the report of Brice et al in 2002,3 which demonstrated that varicose veins were present in half the individuals examined, and that this was the third most frequent feature after distichiasis and lymphoedema. This finding suggested FOXC2 as a possible candidate gene for VV in the wider population, and to test this hypothesis we examined the St Thomas’ twin cohort for linkage and association to the marker D16S520.We also performed the same analyses on haemorrhoids (H) as a potentially associated venous abnormality, and tested whether the two conditions were related in the general population. The phenotypic definition of H did not include those individuals who only manifested during pregnancy, unless stated as such below.
MATERIALS AND METHODS
The study subjects comprised female only twin pairs recruited from the St. Thomas’ UK adult twin registry.10 This is a volunteer registry with twins, ranging from 18 to 80 years of age at first interview, drawn from national media campaigns. The twins completed a range of standardised self and nurse administered questionnaires, and tests providing details of their lifestyles and their medical, drug, obstetric, and gynaecological histories. Their zygosity was determined by a standard questionnaire, and where ambiguous, status was confirmed by genotyping.
The phenotypes were obtained from self administered questionnaire. For varicose veins (VV) the question asked was “Do you suffer from varicose veins?”, and for haemorrhoids the question asked was “Have you ever been affected by piles or haemorrhoids?” A stringent definition of haemorrhoids (H) was used, by which subjects with haemorrhoids were assigned as affected, while those who had suffered from haemorrhoids only during pregnancy were coded as missing. The broader definition of haemorrhoids used included women who reported having had haemorrhoids during pregnancy.
Based on the observed phenotypic association between VV and H, a new combined phenotype (VH) was defined. Individuals were defined as being affected if they reported having VV and H. Individuals with H during pregnancy only were coded as missing.
Genotyping was performed using DNA extracted from venous blood samples. The microsatellite marker D16S520 was genotyped using standard ABI Prism (Applied Biosystems) fluorescence based genotyping methodologies.11,12 Specifically, the marker locus was amplified in 10 μl single plex PCRs in 384 well microtitre plates. Amplification products were pooled and precipitated, and then combined with loading buffer, formamide, and an internal size standard (GeneScan-500; Applied Biosystems); products were then separated by size and were detected using ABI Prism 377 automated sequencers (Applied Biosystems). Twin zygosity and family relationships were rigorously investigated and discrepant pairs discarded from further analyses. The estimated genotyping error rate was <1%.
Casewise concordance was calculated as the probability that a twin is affected given that the other twin is affected. On the assumption that the effects of environmental factors shared by twins are independent of zygosity,13 sometimes known as the “equal environments assumption”, an excess in the casewise concordance estimate in monozygotic (MZ) compared with dizygotic (DZ) twins suggests a genetic aetiological component to a phenotype.
The relative size of the contribution of genetic and environmental factors can be modelled using maximum likelihood variance components or regression based methods. To estimate the proportions of variance in liability attributable to genetic and environmental effects, we used logistic regression to model VV and H, using an extension of the Defries and Fulker method for quantitative traits14 applied to binary traits.15
By assuming an underlying multifactorial normal distribution of small additive risk factors, the phenotypic variance in liability to VV was decomposed into additive polygenic (A) and dominant (D) genetic components, and environmental components shared by both twins (C) and unique to each twin (E). The statistical significance of the variance components A, C, D, and E were assessed by testing whether the regression coefficients, re-parameterised to estimate the attributable proportion of variance, could be dropped from the full model (ACE or ADE) without deterioration in model fit.
The proportion of phenotypic variance attributable to genetic variation (for the best fitting model) is termed “heritability”. For a binary trait, the estimated heritability applies to a postulated underlying latent liability trait with a normal distribution, in which the threshold for affected status is specified by the disease prevalence.15,16 It should be noted that the validity of a liability model is contingent upon the reality of the existence of many small additive risk factors for the trait in question.17 Heritability estimates were adjusted for age, as MZ twins were on average slightly older than DZ twins for the VV sample, and both traits are age related.
The casewise concordance and heritability analyses were implemented in the statistical software package Stata (Stata/SE, Stata, College Station, TX, USA), and the routines used are available upon request (from TA).
Polymorphism information content
To determine the degree of polymorphism of the candidate marker D16S520, the polymorphism information content (PIC) was calculated based on the observed allele frequencies. Generally, marker informativeness increases with the number of alleles at the locus. The PIC of a marker is given by:
Linkage to the marker was tested using twins concordant for the trait (that is, affected sibling pairs; ASP), and sibling pairs discordant for the trait (discordant sibling pairs; DSP). ASP were used to test for linkage using the maximum lod score (MLSASP) method implemented in Genehunter 2.0.18,19 This is a model free, single locus test of linkage for qualitative phenotypes by examining the sharing of marker alleles identical by descent (IBD) among affected sibling pairs. Under the null hypothesis, the estimated proportion of sibling pairs with 0, 1, or 2 marker alleles IBD is 0.25, 0.5, and 0.25. If the candidate marker is linked to a locus predisposing to the phenotype, the observed IBD distribution should have a significant excess of sibling pairs sharing two alleles IBD compared with the null hypothesis of independent Mendelian segregation.
In addition, we contrasted the IBD allele distribution for ASP and DSP using optimal Haseman and Elston methods implemented in Stata using generalised linear modelling20 to generate a Z score (GLMASP/DSP) and p value. Specifically, ASP pairs (n = 102) were coded 1 and DSP coded 0 (n = 266) and this binary trait was regressed upon sibling pair IBD probability using a binomial distribution, robust Huber/White estimate of variance and an identity link function, which assumes an additive genetic model.
MZ twins were not included in the analysis because, if genetically identical, by definition they share two alleles IBD at the marker locus and are therefore uninformative for linkage.
Marker D16S520 has 14 alleles and shows high allelic heterogeneity. A case–control study was performed to test allelic frequencies for significant association with the VV and H phenotypes using χ2 tests. The analyses were conducted using CLUMP.21
Power to detect linkage using variance components linkage test
The power to detect linkage under the assumption of a liability model was carried out using the Genetic Power Calculator software available online.22 We looked at two scenarios. Firstly, we assumed no recombination between the FOXC2 gene and candidate marker D16S520 (that is, θ = 0). Secondly, we set the recombination rate as θ = 0.05, using a Morgan map function and the observed ∼5 cM/Mb recombination rate for chromosome 16q24. The statistical significance level (type I error α) was set at 0.05.
The power to detect linkage is illustrated in table 1. The power is very similar for the two recombination rates assumed, with θ = 0.05 reflecting an approximate physical distance of 80kb for region 16q24. The table indicates the study is well powered to detect linkage for QTLs of moderate to large effect size ranging from 5 to 10% or greater.
A total of almost 2000/3000 female twin pairs replied to a postal survey asking about VV and H (60–66% response rate). Of these, 900 complete DZ pairs were also genotyped for candidate marker D16S520. A description of the twin pairs is shown in table 2. As expected, those unaffected by VV and H were on average younger (46–48 years) than the affected groups (49–55 years). No obvious bias in response was observed for VV, but those responding to questions on haemorrhoids were on average almost 2 years younger than those who did not reply.
Both H, and H including during pregnancy, were significantly related to VV with odds ratios (OR) of 1.42 (95% idence interval (CI) 1.081 to 1.874; p = 0.012) and 1.37 (1.065 to 1.764; p = 0.014) respectively. A new combined phenotype, VH, was therefore defined and tested for linkage and association to D16S520.
The summary of the casewise concordance rates for VV and H are shown in table 3. For both phenotypes, the casewise concordance rates were significantly higher in MZ than DZ twins. This implies a substantial genetic contribution to the aetiology of VV (66% v 42%; p = 1.2×10) and H (68% v 60%; p = 0.02), and VV and H including during pregnancy (74% v 64%; p = 1.2×10).
The results for modelling the heritability of variation in liability are shown in table 4. Best fitting models for the variation in liability to VV and H included parameters for additive genetic (A) and unique environmental factors (E). The estimated additive polygenic heritability value was 90% for VV (95% confidence interval 77 to 100%), 64% for H (52 to –76%) and 59% for H including during pregnancy (49 to 69%).
The summary of linkage results are shown in table 5. There was strong evidence for linkage to the marker for VV, with an excess sharing of two alleles IBD at marker D16S520 for affected sibling pairs (see above). The estimated distribution for sibling pair allele sharing IBD at IBD0: IBD1: IBD2 for ASP (z0:z1:z2 = 0.12:0.5:0.38) was significantly different (GLMASP/DSP Z = −3.17, p = 0.002) from the DSP distribution (z0:z1:z2 = 0.28:0.49:0.23). Testing for linkage using ASP compared with the null hypothesis distribution (0.25: 0.5: 0.25) gave a lod score of MLSASP = 1.37 (p = 0.01). The GLMASP/DSP test was more powerful in this instance, because the study had sufficient numbers to be able to compare ASP and DSP directly (rather than ASP versus the null distribution), in turn reflecting the high prevalence of VV in aging populations.
There was no evidence for linkage to the marker for H or the combined phenotype VH. The degree of polymorphism for marker D16S520 is relatively high, with a PIC value of 0.83.
To test whether the marker was in linkage disequilibrium with a locus that predisposes individuals to VV, H or VH, a case–control association study was performed using a χ2 test. The marker is multi-allelic, with 14 allele groups (table 6), and approximately 87% falling into five allele groups (allele 2, 3, 4, 5, and 6). Table 7 shows that there was no evidence for association with any of the phenotypes.
These data show that VV in the normal population are linked to a functional variant (or variants) in the vicinity of D16S520 on chromosome 16q24. A likely candidate gene in close proximity (80 kb) to D16S520, which may account for the linkage, is FOXC2. It is possible that other genes in the vicinity could be the candidate, including FOXF1, but only FOXC2 has any known association with VV. Further SNP mapping is required to characterise the linkage disequilibrium in the region and to further refine the location of possible functional polymorphism(s).
The estimated heritability in liability for VV is high and in part may explain the high power of this study to detect linkage.23 Similarly, due of the relatively low trait heritability of haemorrhoids and the correspondingly low power to detect linkage (table 1), the possibility of haemorrhoids also being linked to the region cannot be ruled out, although these data show no evidence for this.
FOXC2 is known to be expressed in the paraxial mesoderm and somites of the early vertebrate embryo,24 and is expressed in the heart and blood vessels at later stages.25,26 The mutation of this gene in LD has also clearly demonstrated its importance in the development of the lymphatic system.6–9 All the mutations in the FOXC2 gene that produce LD are consistent with haploinsufficiency,7–9, although a recent publication has suggested that there might be subtle phenotypic differences between patients with mutations that occur within the Forkhead domain and those that occur after this region.27 The syndromic nature of LD, together with a considerable difference in frequency of clinical manifestations (96% distichiasis, 83% lymphoedema, 4% cleft palate,3) suggests that some pathways are much more readily affected by reduction of FOXC2 levels than others. Venous development now appears to be a pathway in which FOXC2 has a developmental role that is only just being recognised, with half of LD patients having VV.3 Venous reflux may well be due to valve abnormalities, and this might suggest a specific role for FOXC2 in development of both venous and lymphatic valves, as there is known to be lymphatic reflux in LD (Brice et al3; Peter Mortimer, personal communication), and the mouse strain heterozygous for a null FOXC2 mutation shows incompetent interlymphangion valves.28 Which particular genes are controlled by FOXC2 in these pathways remains to be elucidated.
The fact that there is linkage to a candidate marker for the FOXC2 gene suggests there is a functional variant(s) within, or in the vicinity of, the FOXC2 gene, which predisposes individuals to varicose veins. The failure to observe association between D16S520 and varicose veins is most probably explained by a lack of linkage disequlibrium between the marker and FOXC2 due to the physical distance (80 kb) involved.29 Alternatively, observed linkage without association might also arise in allelic heterogeneity, as different functional variants could be involved in generating different phenotypic associations in subgroups of people, disrupting any potential linkage disequlibrium. However, in the first instance, a lack of association between the two loci is the more likely explanation.
Possible limitations of this study include whether (a) self diagnosis of VV and H might have produced some misclassification of phenotyping, as no definitions were given; and (b) whether the data from twins is generalisable to the population at large. Although self diagnosis will always tend to be less reliable than clinical diagnosis, varicose veins are easier to self diagnose than many conditions, and it is a condition with which most people are likely to be familiar. The main analytical consequence for inaccurate phenotypic categories would be the loss in power to detect linkage. Secondly, the twins in this study have been shown to be very similar to the population in terms of lifestyle variables such as prevalence of smoking and alcohol consumption, and mean levels and variance for bone mineral density and blood pressure.30 We have no reason to believe that twins might differ from the population in terms of varicose veins.
In conclusion, these results demonstrate that VV and H are heritable, related conditions and the data strongly suggest that FOXC2 may be implicated in the development of VV in the general population. The latter finding is an encouraging result for complex trait research, as it illustrates that at least in some situations there is a continuity between health and disease, and also that candidate genes for syndromes can be of direct relevance to the population as a whole.
This study was funded in part by Wellcome Trust and the CDRF. We would also like to thank Gemini Genomics/Sequenom Inc for providing genotypic data and L Cherkas for help with questionnaires. This research was conducted within the network of the London IDEAS Genetic Knowledge Park. S Jeffery represents the Lymphoedema Consortium, comprising P Mortimer, A Child, G Brice, S Mansour, R Mellor, C Sholto-Douglas Vernon, P Sharpe, M Sarfarazi, A Smith, and K Burnand.
The first two authors contributed equally to this work.