Valerie Vanhooren

Sylviane Dewaele

Claude Libert

Sebastiaan Engelborghs

Peter Paul De Deyn

Olivier Toussaint

Florence Debacq-Chainiaux

Michel Poulain

Youri Glupczynski

Claudio Franceschi

Sebastiaan Engel- Borghs

Peter Paul De Deyn

Koos Jaspers

Ingrid Van Der Pluijm

Jan Hoeijmakers

Cuiying Chitty Chen

Serum N-Glycan Profile Shift During Human Ageing

Keywords: biomarker, ageing, age-related disease, N-glycan profile

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A C C E P T E D

M A N U S C R I P T

Introduction

The ageing process, like many other biological processes, is regulated by classical signaling pathways and transcription factors. Many of these pathways were first discovered in small, short-lived organisms such as yeast, worms and flies, but a remarkable fraction of these pathways turn out to extend lifespan in mammals as well (Kenyon 2010). Lipoprotein metabolism, genomic instability, oxidative stress, cellular response to damage, inflammatory processes, TOR signaling, signals from the reproductive system, telomeres, insulin/IGF-1

A C C E P T E D M A N U S C R I P T

signaling and sirtuins have all been proposed as regulators of the ageing process (Whalley 2009, Kenyon 2010). For instance, DNA damage and repair plays a very important role in the ageing process. To avoid the deleterious consequences of genetic damage, the DNA of all organisms also encodes an intricate network of specialized DNA repair pathways. The vital importance of basic DNA repair mechanisms for maintaining genomic integrity is underscored by their extreme conservation from unicellular bacteria and yeast to man (Hoeijmakers 2001). Defects in DNA repair systems either cause embryonic death or result in severe pathology later in life. Many patients and mouse mutants with defects in different nucleotide excision repair factors exhibit phenotypes reminiscent of premature ageing.

Defects in the transcription-coupled nucleotide excision repair factors CSA or CSB lead to

Cockayne syndrome (CS) in humans (Mitchell 2003).

After middle age, the body starts to experience a decline in various functions, and this leads to a high risk for ageing-associated diseases, such as cardiovascular disease, cancer, arthritis, cataract, osteoporosis, type 2 diabetes, hypertension and Alzheimer's disease. The incidence of all of these diseases increases rapidly with ageing. These age-related diseases can accelerate the ageing process. Dementia encompasses deteriorations in several cognitive domains [START_REF] Morris | Dementia update 2005[END_REF]. It has been defined as a manifestation of deteriorating brain functions over time as a part of ageing due to death or functional decline of brain cells caused by neurodegeneration or other disease [START_REF] Peng | Is dementia a disease?[END_REF]. Alzheimer's disease (AD), the most common type of dementia, is a progressive, degenerative brain disease accounting for 65-70% of all dementia cases [START_REF] Blennow | Alzheimer's disease[END_REF]. The cause of AD is unknown, but ageing is one causative factor. The incidence of AD increases substantially after the age of 70 years and represents a major health concern and economic burden worldwide. Mutations that slow down ageing also postpone age-related diseases such as dementia (Kenyon 2010).

The molecular mechanisms underlying the ageing process are only beginning to be unraveled, yet there is clear evidence that the rate of ageing differs substantially between members of the [START_REF] Anisimov | Rapamycin extends maximal lifespan in cancer-prone mice[END_REF]. To test new interventions such as drugs or food components, it is necessary to measure the effect of the intervention on the underlying process of ageing. Preferably, this should not be limited to the effect on one of the body's systems or functions, for example, on blood pressure or cholesterol levels, but should reveal whether the intervention slows down the ageing process itself. Currently, the only way to test interventions aimed at extending life in humans is to monitor people until the end of their lives to determine if life span has been affected.

Glycosylation is a form of co-translational and post-translational modification. Glycans serve a variety of structural and functional roles in membrane-bound and secreted proteins. An Nglycan (N-linked oligosaccharide) is a sugar chain covalently linked to an asparagine residue of a polypeptide chain. Analyses of protein sequence databases have revealed that about two thirds of those sequences are likely to be N-glycosylated [START_REF] Varki | Essentials of Glycobiology[END_REF]. Glycans can mediate a wide variety of biological functions by virtue of their mass, shape, charge and other physical properties. Unlike protein sequences, which are primary gene products, glycan chain structures are not encoded directly in the genome but are secondary gene products. Glycans

A C C E P T E D M A N U S C R I P T

are synthesized in many different combinations by a variety of competing and sequentially acting glycosidases and glycosyltransferases without use of a template [START_REF] Varki | Essentials of Glycobiology[END_REF].

Small alterations in the environment can cause dramatic changes in the glycans produced by a particular cell, and so N-glycosylation could mirror the status of the cell. We and others recently showed that the composition of serum N-glycans changes during human ageing [START_REF] Knezevic | Effects of aging, body mass index, plasma lipid profiles, and smoking on human plasma N-glycans[END_REF][START_REF] Liu | N-glycomic changes in hepatocellular carcinoma patients with liver cirrhosis induced by hepatitis B virus[END_REF]. Certain relatively specific alterations in the expression of glycans are also found during oncogenic transformation and progression to malignancy, as well as in other pathological conditions, such as inflammation [START_REF] Kramer | Presynaptic alpha-synuclein aggregates, not Lewy bodies, cause neurodegeneration in dementia with Lewy bodies[END_REF]. The spatial and temporal control of glycan expression implies the involvement of glycans in many physiological and pathological processes [START_REF] Haltiwanger | Role of glycosylation in development[END_REF][START_REF] Ohtsubo | Glycosylation in cellular mechanisms of health and disease[END_REF]. Since N-glycosylation reflects in part the state of the cell, it could be a suitable biomarker not only for healthy ageing but also for age-related disease.

To determine whether changes in the extent and nature of posttranslational modification of proteins could be hallmarks of biological age, we investigated the N-glycosylation status of proteins in human blood. We also examined the serum N-glycosylation fingerprint in patients with age-related disease (dementia) and in patients suffering from Cockayne syndrome, a premature ageing syndrome caused by deficiency in a DNA repair. We defined a particular serum N-glycan profile shift during ageing as a biomarker.

Materials and methods

Serum samples from healthy Belgians

Serum samples from 425 healthy Belgian adults were obtained from the Biobank, Institute Born-Bunge, and from the Transfusion Center of the Red Cross, Belgium, in accordance with Red Cross health standards. All of these people were in relatively good health and had no

A C C E P T E D M A N U S C R I P T
relevant acute or chronic disease affecting the immune system. They were divided in four age groups: 20-49 (n=60), 50-75 (n=160), 76-89 (n=30) and 90-99 years (n=175). Informed consent was obtained from all participants.

Serum samples from healthy Italians

Serum samples were obtained from 84 Italian centenarians [START_REF] Liu | N-glycomic changes in hepatocellular carcinoma patients with liver cirrhosis induced by hepatitis B virus[END_REF]. All were in relatively good health and had no relevant acute or chronic disease affecting the immune system. As controls, we studied three groups of healthy Italian volunteers of different ages: 20-49 (n=35), 50-75 (n=37) and 76-99 years (n=71). Informed consent was obtained from all participants.

Serum samples from dementia patients

Serum samples from 79 patients with autopsy-confirmed dementia were included. All sera samples were obtained from the Biobank, Institute Born-Bunge, Antwerp, Belgium. We recorded the scores of the Mini-Mental State Examination (MMSE) at the time of serum sampling, year of clinical diagnosis, clinical and pathological diagnosis, year of death, date of autopsy (Chen C, 2010). The study was approved by the local ethics committee (CME Middelheim Hospital, Antwerp).

Serum samples from patients with Cockayne Syndrome (CS)

Serum samples were obtained from six CS patients aged 3, 7, 8, 10, 10 and 20 years.

Informed consent was obtained from parents or guardians. CS diagnosis was confirmed by functional Nucleotide excision repair (NER) assays in fibroblasts derived from these patients as previously described (Jaspers et al., 2007) A C C E P T E D M A N U S C R I P T

N-glycan fingerprinting

N-glycan fingerprinting of all the samples was done as detailed before [START_REF] Vanhooren | N-Glycan profiling in the study of human aging[END_REF]. N-glycan profiles were generated from 2-µl serum samples by DSA-FACE technology using a capillary electrophoresis (CE)-based ABI3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). The intra-and inter-assay coefficients of variation (CV)

of the glycan analyses are less than 5%.

Statistical analysis

All quantitative variables are presented as mean ± SD. The unpaired t-test was used for comparison of two groups, and for comparison of more than two groups, one-way ANOVA was performed. All reported P-values are two-tailed, and P-values < 0.05 are considered statistically significant. Statistical analysis was performed with SPSS15.0 software (SPSS, Chicago, IL, USA).

Results

The healthy human ageing biomarker GlycoAgeTest: shift of the serum glycoprofile

Our department developed an ultra-sensitive technique to quantitatively profile and sequence N-linked glycans using standard DNA-sequencing equipment [START_REF] Laroy | Glycome mapping on DNA sequencing equipment[END_REF], and we adapted this technique for the study of human ageing [START_REF] Vanhooren | N-glycan profiles as tools in diagnosis of hepatocellular carcinoma and prediction of healthy human ageing[END_REF]. Using "DNA Sequencer Adapted-Fluorophore Assisted Carbohydrate Electrophoresis" (DSA-FACE), an accurate fingerprint of the glycan composition in serum or plasma can be obtained.

The DSA-FACE procedure requires only 2 µl of serum and involves labeling of N-glycans, profiling, and read-out (Fig. 1). A typical N-glycan profile of human serum and plasma was previously generated and the nine most prominent structures corresponding to glycans in the profile were identified previously [START_REF] Liu | N-glycomic changes in hepatocellular carcinoma patients with liver cirrhosis induced by hepatitis B virus[END_REF]. The sizes of the peaks, which represent Using this approach, we previously performed a thorough study of the N-glycosylation status of serum proteins during ageing in 100 healthy Belgian blood donors (Belgian Red Cross, Ghent, Belgium) [START_REF] Liu | N-glycomic changes in hepatocellular carcinoma patients with liver cirrhosis induced by hepatitis B virus[END_REF]. The volunteers were from five age groups (20, 30, 40, 50 and 60 years), each consisting of 10 males and 10 females. Several glycan structures were found to be closely related to ageing. In particular, peak 1 (an agalactosylated core-α-1,6-fucosylated biantennary; NGA2F) increased gradually with age, while peak 6 (a bigalactosylated, core-α-1,6-fucosylated biantennary; NA2F) decreased (Fig. 2). The increases in peak 1 (agalacto-glycan, NGA2F) during ageing could be due to a greater abundance of beta-galactosidase, which removes two galactoses from Peak 6 (NA2F), or to a lesser abundance of beta-1,4-galactosyltransferase, which adds galactose to NGA2F (Fig 2 ).

For precision and accuracy, a ratio of two variables is better than the use of a single measurement. In the current study, we defined the log ratio of NGA2F (peak1) and NA2F

(peak6) as a potential ageing biomarker, named GlycoAgeTest = log[NGA2F/NA2F]

(log[P1/P6]). Using this formula, we reanalyzed its correlation with age in the Belgian volunteers. As shown in Figure 3, the level of GlycoAgeTest remained steady until 40 years, and then increased. This observation demonstrates that the N-glycan profile is rather stable and consistent up to the age of 40 years. It starts to increase by the age of 50 and even more so after the age of 60 years. The intra-and inter-assay coefficients of variations (CVs) of the glycan analysis are less than 5%.

To see whether the shift in serum glycan profile is further increased at even higher later ages, we evaluated the GlycoAgeTest marker in 425 healthy Belgian volunteers in four age groups, including nonagenarians: 20-49 (n=60), 50-75 (n=160), 76-89 (n=30) and 90-99 (n=175). To see whether the GlycoAgeTest marker is also applicable to other populations, we tested a group of 229 healthy Italian volunteers (female n=160, male n=68) of different ages: 20-49 (n=35), 50-75 (n=37), 76-89 (n=33), 90-99 (n=39) and ≥100 years (n=84). The GlycoAgeTest displayed the same trend of change with age observed in the Belgian volunteers (Fig. 4). This demonstrates that the alterations of glycan structures are closely related to the progress of human ageing at least in these two populations. Noteworthy, the level of GlycoAgeTest in the Italian volunteers was significantly lower (p<0.05) than in their age-matched Belgian.

GlycoAgeTest in patients with dementia and patients with Cockayne syndrome

We analyzed serum N-glycans by DSA-FACE in 79 Belgian patients with autopsy-confirmed dementia (31 of 50-75 years and 48 of 76-89 years). The GlycoAgeTest value was significantly higher (p<0.05) in dementia patients than in age-matched controls (Fig. 5), which suggests that the biological age of dementia patients corresponds to a higher chronological age.

On the other hand, progeroid syndromes can occur early in life and accelerate the ageing process, such as Cockayne syndrome (progeroid syndrome), which is a rare autosomal recessive disorder characterized by progressive multisystem degeneration and segmental premature ageing. The average reported life span of patients with CS is approximately 12.5 years (for extensive reviews on clinical characteristics of Cockayne syndrome: (Nance and Berry, 1992) [START_REF] Licht | Cockayne syndrome group B cellular and biochemical functions[END_REF]. Because many of the traits are reminiscent of normal ageing, CS has been classified as a segmental premature ageing syndrome. To analyze the sensitivity of the GlycoAgeTest for a segmental progeroid syndrome, we determined the serum glycan profile in six CS patients (age of 3 to 17 years). As shown in figure 6, the level of GlycoAgeTest in the CS patients is highly heterogeneous, but even the lowest value is much higher than in the controls. The average score is in the same range as that of people in their nineties. Thus, even though the primary defect in CS is due to defective DNA damage repair, it appears that the profile of N-glycans in serum of these patients corresponds with dramatic accelerated ageing.

Discussion

We This difference could be partly due to common lifestyle between Belgium and Italy, and genetic polymorphisms [START_REF] Melzer | Genetic polymorphisms and human aging: association studies deliver[END_REF][START_REF] Agundez | CYP2D6, NAT2 and CYP2E1 genetic polymorphisms in nonagenarians[END_REF][START_REF] Melzer | Genetic variation and human aging: progress and prospects[END_REF]. It was recently documented that plasma N-glycans is affected by age, body mass index, plasma lipid profiles and smoking [START_REF] Knezevic | Effects of aging, body mass index, plasma lipid profiles, and smoking on human plasma N-glycans[END_REF] Moreover, we found that the GlycoAgeTest value was high in patients with dementia compared to age-matched healthy volunteers, which indicates that the biological age of dementia patients is higher than their chronological age. This was also observed in CS. The fact that a defect in the mechanisms of DNA damage repair, and particularly those dealing with DNA lesions that block transcription can lead to a dramatic alteration in the serum Nglycan profile is very interesting. This implies that indirect links exist between defective functioning of specific DNA repair processes and the synthesis or turn-over of specific carbohydrate chains in the ER and Golgi apparatus. In any case, the findings are consistent with the idea that CS is a (segmental) premature ageing syndrome and that DNA damage repair is fundamental to various cellular and organismal processes, including N-glycosylation.

The one common factor between between dementia and CS is the ageing process, and this common association provides a strong argument that our GlycoAgeTest measures a process related to the ageing process.

We previously showed that N-glycan markers can also clearly distinguish between the normal ageing process and disease because sugar structures are changed in different ways and to different extents. For instance, GlycoHCCTest and GlycoCirrhoTest are serum N-glycan diagnostic biomarkers for liver cancer and cirrhosis, respectively, and they were designed

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Vanhooren et al. Serum N-glycan biomarker 13 using different N-glycans from GlycoAgeTest [START_REF] Callewaert | Noninvasive diagnosis of liver cirrhosis using DNA sequencer-based total serum protein glycomics[END_REF][START_REF] Liu | N-glycomic changes in hepatocellular carcinoma patients with liver cirrhosis induced by hepatitis B virus[END_REF]. As the changes in N-glycans are a consequence of genetic and environmental influences and physiologic responses, the GlycoAgeTest marker has a significant diagnostic potential. It could be used in combination with other existing ageing parameters: due to the intrinsic multicausal and multi-system nature of the ageing process, no single measurement has proven capable on its own of determining biological age [START_REF] Barqawi | Observed effect of age and body mass index on total and complexed PSA: analysis from a national screening program[END_REF][START_REF] Johnson | Recent results: biomarkers of aging[END_REF].

The American National Institute of Ageing (NIA) has proposed the following requirements for a biomarker of ageing [START_REF] Ward | Changes in the serum proteome associated with the development of hepatocellular carcinoma in hepatitis C-related cirrhosis[END_REF]. ( 1

  A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT Vanhooren et al. Serum N-glycan biomarker 9 the relative concentrations of the oligosaccharide structures, can be analyzed statistically to determine the correlation between specific glycans and ageing. This high-throughput technique is robust, reproducible, sensitive and quantitative.

  GlycoAgeTest value gradually increased during ageing and reached the highest level in nonagenarians (90-99 years) (Fig 4).

Figure legends Fig. 1 .

 legends1 Figure legends

Fig. 2 .

 2 Fig. 2. A serum N-glycan fingerprint generated by DSA-FACE of a young and an old Belgian healthy person (upper panel). P1 can be converted to P6 by β-1,4galactosyltransferase and P6 can be converted to P1 by β-galactosidase. The structures of P1

Fig. 3 .

 3 Fig. 3. The serum GlycoAgeTest values of healthy Belgian blood donors aged ≤ 60 years. The

Fig. 4 .Fig. 5 .

 45 Fig. 4. The serum GlycoAgeTest values of a population of healthy Belgian (full lines) and

Fig. 6 .

 6 Fig. 6. The GlycoAgeTest performed on patients with Cockayne Syndrome (n=6). The

  Figure 2

  same species, including humans. In other words, "biological age" of an individual might be different from the chronological age, in which case the biological age should be a better indicator of health and life expectancy. Researchers are trying to identify age-related changes that gauge physical age and predict life span more accurately than by counting years. Due to the rapidly increasing average life expectancy in most developed countries, the prevalence of age-related diseases is increasing. This necessitates development of new effective strategies

for prevention and early diagnosis of such conditions. As ageing in many species appears amenable to interventions, research into the process of ageing and means by which it can be slowed down has gained considerable momentum in recent years

  Belgium and Italy) remained unchanged until the age of 40 years, and then gradually increased with age. Interestingly, a lower serum GlycoAgeTest level was detected in the Italians than in age-matched Belgians. This discrepancy might hint a healthier ageing and longer life in Italians compared to Belgians. It has been reported by the

	EHEMU project (European Health Expectancies Monitoring Unit; Healthy Life Years In The
	European Union: Facts And Figures 2005) that the life expectancy at birth is slightly higher in

profiled serum N-glycans by using DSA-FACE to validate our serum N-glycan ageing biomarker, GlycoAgeTest, in different age groups of healthy human volunteers, patients with dementia, and patients with CS, a genetic disease involving progressive neurodysfunction and many features of premature ageing. We demonstrate that the GlycoAgeTest is better than chronological age at estimating the physiological age of human and that it could be used as a biomarker of healthy ageing in humans.

As we age, the functions and physiological reactions of our body gradually slow down, and this becomes apparent mainly after middle age. We found that the level of our GlycoAgeTest in the two European nations (Italy (78.2 years for men and 84 years for women) than Belgium (76.2 years for men and 81.9 for women) (http://ec.europa.eu/health/archive/ph_information/reporting/docs/hly_en.pdf).

  . It is very interesting to further validate whether the GlycoAgeTest could identify accelerated ageing individual who is associated with lifestyle factors (such as smoking, obesity, lack of exercise and etc). The GlycoAgeTest is being validated further in 14 European nations in the context of the MARK-age project (www.mark-age.eu), a large-scale integrated project supported by the European Commission.

  ) It must predict the rate of ageing, i.e. it should determine how far one has progressed in his or her total life span, and it must predict this more accurately than chronological age. (2) It must monitor a basic process that underlies the ageing process, not the effects of disease. (3) Repeated testing should be harmless and should use, for example, a blood test or an imaging technique. (4) It should be possible to use it in humans and in laboratory animals, such as mice, so that it can be tested in lab animals before being validated in humans. Our GlycoAgeTest fulfills most of these criteria, and it is also extremely robust. Serum samples kept for years in the freezer yielded the same N-glycan profile even after numerous freeze-thaw cycles. Only 2 µl of serum is needed for generating the N-glycan profile, and so the marker can be easily repeated without harming the person.In summary, measuring the shift in serum N-glycan profile could be used to regularly evaluate the overall age-related health status. This not only facilitates the diagnosis of age-related diseases such as dementia, but also the follow-up of different therapies in individuals not suffering from accelerated ageing. As ageing is the biggest risk factor for all age-related diseases, which emerge mainly in middle age or later, detection of an abnormal N-glycan fingerprint and a faster-than-normal rate of ageing in a given would indicate the need for additional diagnostic and preventive measures. With a reliable GlycoAgeTest biomarker to predict biological age and monitor healthy ageing progress, it could become possible to take preventive measures to slow down ageing and age-related diseases.

	M A N U S C R I P T
	Moreover, this marker has already been validated in mice (unpublished data). A C C E P T E D
	The mechanism behind these age-related N-glycosylation alterations requires investigation.
	Preliminary data indicate that the N-glycosylation changes are not limited to a certain serum
	protein but cover all serum glycoproteins. This need further study. Whether the N-
	glycosylation alterations are caused by differences in biosynthesis or clearance or by other
	factors is currently under investigation and should provide us with many interesting answers
	about the mechanism behind our GlycoAgeTest.

ACCEPTED MANUSCRIPTVanhooren et al. Serum N-glycan biomarker

Acknowledgements

We thank Dr. Amin Bredan for editing and Dr. Alexander Buerkle for critical reading the manuscript. This work was supported by a grant from Ghent University, Belgium (BOF No. 01106205), European commission FP6 Linkage (FP6-513866), European Commissionsponsored project Geha (Genetics of Healthy Ageing), (LSHM-CT-2004-503270), Markage (FP7-Health-2008-200880), DNA Repair (LSHG-CT-2005-512113) and LifeSpan (LSHG-CT-2007-036894), National Institute of Health (NIH)/National Institute of Ageing (NIA) (1PO1 AG-17242-02), NIEHS (1UO1 ES011044), and Netherlands Organization for Scientific Research (NWO) through the foundation of the Research Institute Diseases of the Elderly.

Abbreviations

DSA-FACE: DNA sequencer-assisted, fluorophore-assisted carbohydrate NGA2F: agalactosylated, core-α-1,6-fucosylated biantennary NA2F: bigalactosylated, core-α-1,6-fucosylated biantennary CS: Cockayne Syndrome


Aude Saint

Jean-Marc Kaufman

Agnes Ostertag

Martine Cohen-Solal

Anne Boland

Kaatje Toye

Diana Zelenika

Mark Lathrop

Christine De Vernejoul

Maria Martinez

Saint- Pierre Aude
email: pierre@inserm.fr

Kaufman Jean-Marc

Ostertag Agnes

Cohen-Solal Martine

Boland Anne

Toye Kaatje

Zelenika Diana

Lathrop Mark

Vernejoul Marie- Christine

Martinez Maria

Bivariate association analysis in selected samples: Application to a GWAS of two Bone Mineral Density phenotypes in males with high or low BMD

Keywords: Bivariate association, GWAS, BMD, Osteoporosis

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

With the availability of high-density maps of single nucleotide polymorphisms (SNPs), association studies have become popular tools for identifying genes underlying complex human traits and diseases. For most current population-based genome-wide association studies (GWAS) statistical power is often limited due to the complex interplay among factors that influence the etiology of diseases [START_REF] Hirschhorn | Genome-wide association studies for common diseases and complex traits[END_REF] . Increasing sample size and multilocus or multivariate statistical analyses can improve the power for detecting association.

Sample size is often restricted due to genotyping costs and limited sample resources.

Several studies have demonstrated that analyzing samples selected with extreme values can be more powerful than analyzing samples randomly selected from the population 2-4 .

In addition to using selected samples, another approach to increasing association test power is to perform joint analysis of multiple correlated phenotypes. For many common multifactorial traits, several correlated phenotypes are usually recorded for each individual during sample collection, but most often the phenotypes are analyzed separately in a univariate framework. Joint analysis of correlated phenotypes can theoretically provide greater power than that provided by analysis of individual phenotypes [START_REF] Allison | Multiple phenotype modeling in gene-mapping studies of quantitative traits: power advantages[END_REF][START_REF] Amos | Comparison of multivariate tests for genetic linkage[END_REF][START_REF] Almasy | Bivariate quantitative trait linkage analysis: pleiotropy versus co-incident linkages[END_REF][7] . Multivariate analysis can also alleviate the multiple testing problem, caused by testing different traits separately, and thereby improve the ability to detect genetic variants whose effects are too small to be detected in univariate analysis 8 . Several multivariate approaches have been applied to linkage studies of correlated complex phenotypes, as osteoporosis and bone-related phenotypes [START_REF] Wang | Chromosome 2q32 may harbor a QTL affecting BMD variation at different skeletal sites[END_REF][START_REF] Pan | Chromosomal regions 22q13 and 3p25 may harbor quantitative trait loci influencing both age at menarche and bone mineral density[END_REF][START_REF] Wang | Bivariate whole genome linkage analyses for total body lean mass and BMD[END_REF][START_REF] Liu | A bivariate whole genome linkage study identified genomic regions influencing both BMD and bone structure[END_REF] . Similarly, various methods, often based on Generalized Estimating Equations (GEE), have been proposed for performing multivariate association tests on population-or family-based data [START_REF] Lange | A multivariate familybased association test using generalized estimating equations: FBAT-GEE[END_REF][START_REF] Lange | A family-based association test for repeatedly measured quantitative traits adjusting for unknown environmental and/or polygenic effects[END_REF][START_REF] Jung | Bivariate combined linkage and association mapping of quantitative trait loci[END_REF][16][START_REF] Pei | Multivariate association test using haplotype trend regression[END_REF][START_REF] Yang | Bivariate association analysis for quantitative traits using generalized estimation equation[END_REF][START_REF] Zhang | Family-based bivariate association tests for quantitative traits[END_REF][START_REF] Zhang | Univariate/multivariate genomewide association scans using data from families and unrelated samples[END_REF] . Of the two studies that have investigated the power of bivariate association test in population-based data, one applied the restricted bivariate association test that assumes same Quantitative Trait Locus effects on each trait 16,[START_REF] Yang | Bivariate association analysis for quantitative traits using generalized estimation equation[END_REF] . Such constraints in the model may have overestimated or underestimated the relative performance of bivariate over univariate analysis. Finally, GWAS studies using multivariate analysis are rare, especially in samples of subjects selected through their phenotype values, and further investigations using this approach are warranted [START_REF] Abecasis | The power to detect linkage disequilibrium with quantitative traits in selected samples[END_REF] .

To this aim, we evaluated the statistical properties of joint association analysis of two correlated quantitative traits in samples of unrelated subjects through simulation studies, using the Seemingly Unrelated Regression (SUR) bivariate model that allows for different QTL effects on traits. The evaluation was conducted under different situations according to the sample selection design, genetic effects and residual correlation between the traits. We demonstrate the efficacy of SUR-based bivariate test by applying it to simultaneous GWAS analysis of two correlated bone phenotypes, Bone Mineral Density (BMD) at the Lumbar Spine and at the Femoral Neck, which are major risk factors of osteoporosis.

METHODS

SUR-based bivariate model

The Seemingly Unrelated Regression (SUR) model [START_REF] Zellner | An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias[END_REF] is a generalization of a classical linear regression model that consists of several regression equations with potentially different sets of explanatory variables. It thus allows for a differential effect of explanatory variables on phenotypes as well as the possibility that some variables might be associated with only one trait. Let's N be the total number of unrelated subjects (i=1,..,N), each having observations on two phenotypes y ji (j=1,2). Consider a system of 2 equations, where the jth equation is of the form: y j = X j x β j + e j ; y j is a Nx1 vector of the phenotypic values, X j is a (K j +1)xN matrix of explanatory variables with K j representing the number of explanatory variables in the model for phenotype j excluding the intercept; 0 1 ( , ,..., )

j j j j t j K β β β β =
is the (K j +1)x1 vector of coefficients and e j is a Nx1 vector of the residuals errors. The system of SUR can be written as:

1 1 1 1 2 2 2 2 0 0 y e X X y e β β ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ = + ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝ ⎠ (1)
The SUR model allow for cross-equation correlation of the residual terms. The covariance matrix of all the residuals is assumed to be normally distributed with mean 0 and covariance matrix ( )

t N E ee I = Ω = Σ ⊗ where N
I is a NxN unit matrix and Σ a 2x2 matrix with the following form:

2 1 12 2 1 2 2 rE rE σ σσ σ σ σ ⎛ ⎞ × Σ = ⎜ ⎟ × ⎝ ⎠ (2) 2 1
σ and 2 2 σ are the residual variances of Y 1 and Y 2 respectively and rE is the residual correlation between Y 1 and Y 2 .

The SUR model is estimated using the generalized least squares method where the covariance matrix Ω is first estimated using ordinary least squares regression in system (1). Linear restrictions on coefficients can be tested by an F test. The F statistic for systems of equations is:

F = [ (Cβ') t (C cov (β') C t ) -1 (Cβ')] /2
where, C is the matrix of restrictions on coefficients. Under the null hypothesis, the F statistic has a central Fisher distribution with 2 and 2xN-K degrees of freedom where K is the total number of estimated coefficients (K=K 1 +K 2 +2). The goodness of fit of the whole system can be measured by the McElroy's r-square (R 2 ). R 2 is the proportion of co-variance due to X taking into account the residual matrix covariance Ω [START_REF] Mcelroy | Goodness of Fit for Seemingly Unrelated Regressions[END_REF] .

Here, we applied the SUR model to test association to two continuous phenotypes in unrelated subjects genotyped at one SNP marker, and X j is the Nx1 vector of genotypes at the SNP. Under an additive model, the genotype for each individual i, noted g i , is coded as a function of the number of minor alleles, that is, 0, 1 or 2. We computed the SUR model free of constraints on the regression coefficients, that is, β1 and β2 were freely estimated. Under the null hypothesis of no association to either one or both phenotypes, the F statistic has a central Fisher distribution with 2 and 2 x (N-2) degrees of freedom. Separate association analyses of Y 1 and Y 2 can be conducted using traditional univariate linear regression model: y j = g x β j + e j , where y j , g, and β j are as described above but now e j is assumed to follow a normal distribution N (0, 2 j σ ). The null hypothesis of no association (β j = 0) can be tested against the alternative (β j ≠ 0) with a Student statistic (t-test) with N-2 degrees of freedom.

Simulation study:

We considered genetic models of complex traits and specifically tried to generate correlated data mimicking as much as possible our real Bone Mineral Density (BMD) GWAS data (see below). Since a strong (~0.5) and positive phenotypic co-variation exists for BMD values at the Lumbar Spine (LS) and at the Femoral Neck (FN) [START_REF] Livshits | Genetics of bone mineral density: evidence for a major pleiotropic effect from an intercontinental study[END_REF] , we generated data for two positively correlated quantitative phenotypes. Further, in real datasets, as causal loci usually contribute a small proportion to the total phenotypic correlation, residual correlation approximates phenotypic correlation between traits. It is also more realistic to assume that the investigator has a priori knowledge on the magnitude and sign of the co-variation of the studied phenotypes than on the magnitude and sign of the QTL effect on each phenotype. Therefore, in all our scenarios, the sign of the residual correlation (rE) was positive, but the sign of the induced QTL correlation (rG) was either positive or negative. Also, our BMD GWAS study used a sampling design, with extreme truncate selection of unrelated males, aiming to improve power.

Therefore, we also generated samples of subjects drawn from the extremes of the phenotype(s) population distribution.

The main scenarios and parameter settings are shown in Table 1. The different settings allowed us to generate data for a QTL having same or different effect on the two positively correlated phenotypes, and the two sources of co-variation (QTL and residual) have same or opposite sign. Briefly, we assumed a bi-allelic QTL having additive effects (a j ) on Y j (j=1,2), with minor and major allele frequency q and p, respectively. The QTL contribution to Y j is the trait-specific QTL heritability, h 2 j . Here, we focussed our power investigation to QTLs explaining a relatively small part of the trait variance, i.e., from 0.5% to 3% which, for complex traits, seemed to us more realistic. The genotypic means (m jk ) of Yj are equal to 2q x a j , (q-p) x a j and -2p x a j when k, the number of minor alleles, is equal to 0, 1 and 2 respectively and with a j = √[h 2 j ./2pq]. We varied the sign of a j : both were of same or opposite sign and the QTL correlation (rG) was, thus, equal to +1 or -1, respectively. We first generated samples of subjects unselected for their traits values (denoted as Su). Second, we generated subjects selected from the 2.5% (i.e., trait value ≤-2) and 30% (i.e., trait value >0.5) left and right tail of the population distribution of Y 1 (denoted as S1). Third, we included Y 2 in the selection design, that is, we selected subjects from the 2.5% and 30% left and right tail of the population distribution of Y 1 and Y 2 (denoted as S2). These truncate selection criteria (trait value ≤-2 or >0.5) are the values we have used in our real BMD GWAS. Under S1 and S2, we generated samples with equal number of subjects drawn from the left (N/2) and the right (N/2) side of the phenotypes distributions.

Traits values of N (300, 1000) unrelated subjects were generated as follows. For a given combination of parameter values (rE, h 2 1 , h 2 2 , rG), we first draw QTL alleles from a binomial distribution with parameter q, and built genotypes under Hardy-Weinberg equilibrium. Then, conditionally on the generated genotype, g k (k=0,1,2), we jointly drew the values of Y 1 and Y 2 via a bivariate normal distribution with mean (m 1k , m 2k ) t and variance matrix Ω, given in equation (2). Third, under sampling S1 or S2, we applied the corresponding truncate selection, that is individuals not fulfilling the selection criteria were withdrawn from the sample. Steps 1 to 3 were repeated until reaching the required left and right truncated sample sizes of (N/2) subjects.

Each replicate was analyzed with SUR-based bivariate and with two separate univariate analyses using the systemfit package of R software (http://www.r-project.org/) using the genotypes at the QTL, that is, the SNP is the causal variant. The mean and standard deviation of each association statistic (F test and t 1 , t 2 tests) were derived from K replicates. Power and type I error rates of each association test were calculated as the proportion of replicates with a test statistic exceeding a given theoretical threshold (Rα) value, at nominal significance levels, α=5%, 1%, 0.1% and 10 -5 . Type 1 errors were estimated in the settings were h 2 1 =h 2 2 =0 with K=20 000 replicates. Power rates were derived with K=1 000 replicates. To compare the performance of bivariate and that of univariate association analysis, we computed the proportion of replicates where t 1 and t 2 were both lower than Rα. One minus this proportion estimated the probability to detect association to either one of the two phenotypes. To adjust for the two univariate association tests, we applied the Bonferroni correction, that is, we used the theoretical thresholds Rα/2.

RESULTS

SIMULATION STUDY

Tables 2 and 3 present the mean (and sd) association statistic of the SUR-based bivariate (F test) and of the traditional univariate tests (t test), respectively when N=1 000 for 66 scenarios under the alternative hypothesis and when q=0.4. For a given QTL heritability value, the results did not vary, as expected, with q.

Bivariate association statistics: In randomly selected samples, the results in Table 2 show several well-established power figures. First, mean F statistics of bivariate association analysis increase with the size of the trait-specific QTL heritability (h 2 1 and/or h 2 2 ) irrespective of rG and rE. Second, the power is highest in presence (rG≠0) than in absence (rG=0) of pleitropic effects: the highest power is achieved when rG=-1, that is, when the correlation induced by the QTL effect and the residual correlation are opposite in sign. Third, the results also confirm that the power of bivariate association test varies with the size of the residual correlation: when rG=0 or rG=-1, the power increases with rE; conversely, when rG=+1 it decreases with rE. These general trends are observed irrespective of the sampling selection designs. Applying extreme truncate selection increases the power of bivariate association analysis, but the optimal selection design depends on the true genetic model. When rG=0 or rG=-1, extreme selection on one trait (S1) is more efficient than extreme selection on both traits (S2). Conversely, when rG=+1, S2 is more efficient than S1. Overall, under Su or S1, the highest mean F statistics are obtained when rG=-1, irrespective of rE. Under S2, the highest power is achieved when rG=+1 or when rG=-1, depending on the size of rE. Interestingly, when the traits are moderately (rE=0.20) correlated, mean F statistics have greater values when rG=+1 than when rG=-1.

Univariate association statistics

Table 3 shows again several well-established power figures. In randomly selected samples, the power of univariate analysis increases with the QTL heritability (h 2 1 /h 2 2 ) and varies little with the size of the residual correlation, rE. For phenotype Y 1 , under a given QTL heritability (h 2 1 ) value, the mean statistic values of all models are similar in the randomly selected samples. Applying extreme truncate selection increases the power of univariate association analysis of Y 1 . Under S1, the power remains similar whichever rG. Under S2, the power is the highest and the lowest for the pleiotropic models rG=+1 and rG=-1, respectively. When rG=-1 or rG=0, the power of univariate association analysis is greater under S1 than under S2. The reverse trend is obtained when rG=+1. For phenotype Y 2 , the power of univariate analysis depends on rG and rE. Further, applying extreme selection does not always lead to a gain in power. Indeed, when rG=-1 the power of univariate analysis is the greatest in the unselected samples (Su). When rG=0 the mean t statistic values in the selected samples are biased and inflated. The magnitude of the bias is greater under S2 than under S1. Under S1, the bias increases with rE.

Overall, applying selection criteria on one or both traits is an optimal sampling design when rG=+1: the power of each separate univariate analysis is improved over that in randomly selected samples. When rG=-1, applying extreme truncate selection leads to both a substantial gain and decrease in power for Y 1 and Y 2 , respectively. For the situations in which the QTL does not exert pleiotropic effects (rG=0), the highest power of univariate analysis of Y 1 is obtained in the selected samples. However, the mean t statistic values for Y 2 , the trait no associated to the QTL, are also increased. Type I error rates of separate univariate analyses may thus be inflated, especially in selected samples and when the residual correlation is high. 1.B). However, the false-positive rates of univariate association analyses for detecting association to either or both the two traits are, as expected, inflated: the estimated rates are roughly two times higher than the theoretical rates. Applying a Bonferroni correction (denoted as U_b) leads to slightly conservative significance levels, especially when the residual correlation between the traits is strong.

Power comparisons :

The power to detect association to either or both of the two traits using SUR-based bivariate analysis was compared to the power of separate univariate analysis of Y 1 and Y 2 adjusted for multiple testing by the Bonferroni correction (denoted as U_b). Figure 1.A shows the power curves (at significance of 10 -5 ) against the QTL heritability (h 2 1 , h 2 2 ) when N=1 000, for moderately (rE=0.2) or strongly (rE=0.6) correlated traits. Power curves under S1 and S2 are shown in Figure 1.B, when h 2 1 =h 2 2 =0.005, N=1 000 and rE=0.2 or 0.6.

In randomly selected samples (Figure 1.A), the relative advantage of SUR-based bivariate over univariate association analysis is more obvious when rG=-1 and/or the traits are strongly correlated (rE=0.6) but also when rG=+1 and the traits are moderately correlated (rE=0.2). Under S1 (Figure 1.B), SUR-based bivariate is slightly less powerful than univariate analysis when rG=+1 and rE=0.6 or when rG=0 and rE=0.2.

For strongly correlated traits, the power rates are equal to 94.5% (SUR) vs. 29.3%

(U_b) when rG=-1; 44.0% (SUR) vs 32.3% (U_b) when rG=0; 36.8% (SUR) vs 39.9% (U_b) when rG=+1. For moderately correlated traits, the power rates are equal to 64.6%

(SUR) vs 31.7% (U_b) when rG=-1; 32.9% (SUR) vs 34.9% (U_b) when rG=0; 43.7% (SUR) vs 32.6% (U_b) when rG=+1. Under S2 (Figure 1.B), SUR-based bivariate shows same or slightly lower power than univariate analysis except when rG=-1 or when rG=0 and rE=0.6 where it outperforms univariate test. As already noted above selecting on Y 1 (S1) is the most efficient sampling design when rG=-1 or when rG=0 and the traits are strongly correlated (rE=0.6). Selecting on both traits (S2) is the most efficient design when rG=+1. Overall, when rE=0.6, the power of SUR is the greatest (94.5%) when rG=-1 and under S1, while the power of univariate analysis is the greatest (56.8%) when rG=+1 and under S2. When rE=0.2, the power of SUR and univariate analysis are both the greatest (72.5% and 72.9%) when rG=+1 and under S2. As shown in Supplementary Table 2, all these trends are confirmed under various parameter settings.

ANALYSES OF EMPIRICAL BMD-GENOME-WIDE ASSOCIATION DATA

BMD GWAS data: Subjects were recruited from the Network in Europe on Male

Osteoporosis study [START_REF] Pelat | Complex segregation analysis accounting for GxE of bone mineral density in European pedigrees selected through a male proband with low BMD[END_REF][START_REF] Kaufman | Genome-wide linkage screen of bone mineral density (BMD) in European pedigrees ascertained through a male relative with low BMD values: evidence for quantitative trait loci on 17q21-23, 11q12-13, 13q12-14, and 22q11[END_REF] . Subjects selected from this cohort were unrelated males > 18

and < 68 years of age. In addition, the subjects were selected by bone densitometry (measured at the Lumbar Spine and at the Femoral Neck) criteria, having either low BMD (LS-Z-scores ≤-2, n=175) or high BMD (both LS-and FN-Z-scores >0.50, n=155). Further details of the study sample are provided in Supplementary Table 3.

Genotyping was carried out at the Centre National de Génotypage (CNG, Evry, France) using the Illumina 370K platform. SNPs and DNA data were subjected to standard quality control analyses with PLINK [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF] (details are provided in Supplementary Methods).

Association analysis: Our primary analysis was the joint association analysis of LS- Our study used a design, with extreme truncate selection of unrelated males, aiming to improve power. The approach of studying samples drawn from the extremes of the population distribution of BMD has been used in several linkage studies of BMD variation [START_REF] Kaufman | Genome-wide linkage screen of bone mineral density (BMD) in European pedigrees ascertained through a male relative with low BMD values: evidence for quantitative trait loci on 17q21-23, 11q12-13, 13q12-14, and 22q11[END_REF][START_REF] Sims | Genetic analyses in a sample of individuals with high or low BMD shows association with multiple Wnt pathway genes[END_REF] , but rarely in association studies [START_REF] Kung | Association of JAG1 with bone mineral density and osteoporotic fractures: a genome-wide association study and followup replication studies[END_REF] , and to our knowledge, never in samples drawn from the population of males. Due to our relatively small GWA sample size, no SNP showed evidence of association to either one or both BMD phenotypes at genome-wide significance threshold of 1.7x10 -7 (0.05/ 298 783 SNPs). However, we used an extreme truncate selection design that, as shown by our simulation studies, has increased power over unselected samples. Our SUR-based bivariate association analyses identified strong association (P<8.4x10 -6 ) with 3 genomic regions (6q22.1, 15q14 and 22q11). These SNPs have not yet been reported to be associated with bone density in previous GWAS [START_REF] Richards | Collaborative meta-analysis: associations of 150 candidate genes with osteoporosis and osteoporotic fracture[END_REF][START_REF] Rivadeneira | Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies[END_REF][START_REF] Styrkarsdottir | Multiple genetic loci for bone mineral density and fractures[END_REF] . Two of them, on 15q14-15 and 22q11, are located in genes that are known to be expressed in skeletal muscle [START_REF] Doege | Characterization of human glucose transporter (GLUT) 11 (encoded by SLC2A11), a novel sugar-transport facilitator specifically expressed in heart and skeletal muscle[END_REF][33] : GLUT 11 encoded by SLC2A11 on 22q11 and RYR3, on 15q14-15. Because muscle contraction has a major impact on bone density, this might represent an indirect role of these genes on bone density.These genetic variants, whether they are site-specific or possibly shared (pleiotropic), may warrant further follow-up genetic studies on BMD and other bone-related phenotypes.

DISCUSSION

We have evaluated the performance of bivariate association analysis based on the Seemingly Unrelated Regression (SUR) model, which allows different genetic models for different traits. To our knowledge, this is the first study to specifically derive the power and the relative performance of bivariate association analysis in selected samples of unrelated subjects. Our main results coincide with well-known power figures [START_REF] Almasy | Bivariate quantitative trait linkage analysis: pleiotropy versus co-incident linkages[END_REF][7][8] and confirmed that bivariate association analysis outperforms univariate analysis when the QTL exerts pleiotropic effects and the relative increase in power is greatest when correlation of the QTL is opposite in sign to the residual correlation. The most powerful sampling selection design varied with the genetic model, specifically with the size and the direction of the induced-QTL correlation. Applying truncate selection on one trait was found the most efficient sampling design when the genetic and the residual correlations are opposite in signs. The same most efficient design was found when the QTL does not exert pleiotropic effects: the power of the SUR-based bivariate association test was found as good as or better than that of univariate association test, depending on the size of the residual correlation. Finally, when the QTL exerts pleiotropic effects and both sources (QTL and residual) of co-variation are of same sign, applying selection criteria on both traits was found the optimal sampling selection design. Under this sampling design, the performance of SUR-based bivariate test relatively to univariate analysis decreases with the size of the residual correlation.

So far, two studies have investigated the power of bivariate association in unselected population-based data, and they both applied bivariate association test based on Generalized Estimating Equations 16,[START_REF] Yang | Bivariate association analysis for quantitative traits using generalized estimation equation[END_REF] . The former applied a general GEE-based model that allows, as the SUR model, for different QTL effects on the two traits. The second study used a GEE-based bivariate model that assumed same QTL effects on the phenotypes. Our results are congruent with those reported by the first study. The restricted bivariate test estimates, as the univariate test, a single parameter (i.e., the SNP regression coefficients on each trait are all set equal). Under the restricted bivariate model, the gain in power of bivariate analysis is enhanced and reduced when the QTL has similar effect and when it affects one trait only, respectively. Clearly, rarely, knowledge of this magnitude about a complex trait is known a priori. Thus, we do not recommend using restricted bivariate models even in unselected data.

Our bivariate genome-wide association analysis of Lumbar Spine and Femoral Neck BMD values, conducted in a sample of unrelated males with low BMD (LS Zscores ≤-

2) and high BMD (LS and FN Zscores >0.5), consistently demonstrated the advantage of the SUR-based bivariate test over separate univariate analysis. All of the top hits in univariate analysis also showed strong evidence of association in bivariate analysis.

Conversely, additional SNP associations were detected with the bivariate method that did not reach nominal significance in single-trait analyses: this was achieved without adjusting significance of univariate analyses for multiple testing.

In conclusion, our results showed that SUR-based models are useful to detect association for correlated phenotypes. However, our results also showed that similar power levels can be achieved whether the QTL exerts or not pleiotropic effects. Thus, disentangling pure pleiotropic from residual covariation remains a challenge even in bivariate association analysis. Bertocchini F, Ovitt CE, Conti A et al: Requirement for the ryanodine receptor type 3 for efficient contraction in neonatal skeletal muscles. Embo J 1997; 16: 6956-6963.

Titles and legends to figures

Figure 1: Power rates at α=10 1 Minor allele; 2 Minor allele frequency; 3 r-square of the whole system taking into account the residual (co)variance matrix; 4 rank of the identified SNP; 5 r-square from linear regression; [START_REF] Amos | Comparison of multivariate tests for genetic linkage[END_REF]

  When the QTL/SNP has no effect on Y 1 and Y 2 , the values of the mean and standard deviation of both bivariate and univariate association tests are close to the theoretical values, regardless of the residual correlation, minor allele frequency of the studied SNP and of the selection sampling design (Supplementary Table1.A).

Type I error rates:

Indeed, SUR-based bivariate and each separate univariate association tests have correct type I error rates (Supplementary Table

Table 2 :

 2 Mean (and sd) of the SUR-based bivariate association statistic (F test) in

	-5 of SUR-based bivariate analysis and univariate analysis
	adjusted for multiple testing by Bonferroni correction (U_b), in samples of N=1 000
	subjects and under various parameters settings: QTL heritability (h 2 1 /h 2	2 ), sign of the
	induced genetic correlation (rG), residual correlation (rE).
	(A) Power estimates against QTL heritability for moderately (rE=0.2) or strongly
	(rE=0.6) correlated traits, in randomly selected samples (Su)
	(B) Power estimates under extreme selection (S1 or S2) for moderately (rE=0.2) or
	strongly (rE=0.6) correlated traits and QTL heritability (h 2	1 =h 2	2 =0.005)
	Figure 2: Overlap in significance of results from bivariate and univariate (Best_U)
	association analysis.
	(A) Top 100 hits in SUR-based bivariate association test: -log 10 P-values of univariate
	analysis against -log 10 P-values of SUR-based bivariate analysis
	(B) Top 100 hits in univariate association test: -log 10 P-values of SUR-based bivariate
	analysis against -log 10 P-values of univariate analysis

Table 3 :

 3 Mean (and sd) of the traditional univariate association statistic (t test) in

	samples of N=1 000 subjects for various parameter settings: QTL heritability (h 2	1 /h 2	2 ),
	sign of the induced genetic correlation (rG), residual correlation (rE), and sampling
	selection design.

from Univ. Of FN Top from SUR analysis Univ. LS Univ. FN SUR Top from Univ. of LS

Table 4 :

 4 Association results: Top 10 most associated SNPs from SUR-based bivariate or from separate univariate analysis of LS and FN BMD

1 :

  Unadjusted univariate P values.

					A -rE=0.2								A -rE=0.6
		1.0							1.0
		0.9							0.9
		0.8							0.8
		0.7							0.7
	Power	0.4 0.5 0.6							0.4 0.5 0.6
		0.3							0.3
		0.2							0.2
		0.1							0.1
		0.0							0.0
		0.5%	1.0%	0.5%	0.5%	1.0%	0.5%	0.5%	1.0%	0.5%	1.0%	0.5%	0.5%	1.0%	0.5%	0.5%	1.0%
		0.0%	0.0%	0.5%	1.0%	1.0%	0.5%	1.0%	1.0%	0.0%	0.0%	0.5%	1.0%	1.0%	0.5%	1.0%	1.0%
		rG=0		rG=1			rG=-1		rG=0		rG=1			rG=-1
					SUR U_b

Su: unselected sample; S1: sample selected on Y 1 distribution; S2: sample selected on Y 1 and Y

distributions.

a j = √[h

j ./2pq], where q is the Minor Allele Frequency

ACKNOWLEDGMENTS

Part of this work was supported by the Network in Europe on Male Osteoporosis (European Commission grant QL6-CT-2002-00491), by the Flemish Fund for Scientific Research (FWO Vlaanderen grants G.0331.02 and G.0662.07); by the Société Française de Rhumatologie (SFR), and by the French National Agency of Research (ANR).

Conflict of interest: none

QTL Model

Heritability (h 2 j ) [START_REF] Hirschhorn | Genome-wide association studies for common diseases and complex traits[END_REF] effect size (aj) rG q h 2 j rE N Sampling Design