Eleonor I Fransson
email: eleonorfransson<e.fransson@ucl.ac.uk

G David Batty

Adam Tabák

Eric J Brunner

Meena Kumari

Martin J Shipley

Archana Singh-Manoux

Mika Kivimäki

Adam G Tab

Mika Kivim

Association between change in body composition and change in inflammatory markers: an 11-year follow-up in the Whitehall II Study

Keywords: Weight, Body Mass Index, waist circumference, inflammation, CRP, IL-6 Adult, Biological Markers, blood, Body Composition, physiology, Body Height, Body Mass Index, Body Weight, C-Reactive Protein, metabolism, Cohort Studies, Female, Follow-Up Studies, Humans, Inflammation, blood, epidemiology, metabolism, Interleukin-5, blood, Male, Middle Aged, Obesity, blood, complications, metabolism, Waist Circumference, Weight Gain, physiology

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

 2 7Low-grade systemic inflammation is associated with atherosclerosis and other chronic conditions such as diabetes ( , ). C-reactive 1 2 protein (CRP, a plasma protein synthesized by the liver) and interleukin-6 (IL-6, a cytokine that governs inflammatory cascades and modulates CRP) are key biochemical markers of inflammation ( , ). It is well established that overweight and obese persons have, on 3 4 average, higher concentrations of CRP and other inflammatory markers than their leaner counterparts. ( , ). It has also been shown that 5 6 intentional weight loss leads to reduced inflammatory marker levels in the short term ( , ). However, the long-term associations of 7 8 change in weight or central obesity with changes in inflammatory markers in larger populations are poorly understood ( , ). 9 10

Accordingly, we assessed how changes in weight, BMI, and waist circumference over a period of 11 years were associated with changes in CRP and IL-6 in a large cohort of middle-aged British adults ( ). We also examined whether weight-related changes in inflammatory 11 markers differed between normal weight, overweight and obese individuals.

MATERIALS AND METHODS

Study population

Whitehall II, a prospective cohort study, was established in 1985. Participants were all office staff, aged 35 55 years at enrolment and working in 20 London-based Civil Service departments. In total, 10,308 people (6895 men and 3413 women), 73 of those invited, were % included at baseline (1985 1988) ( ). Initial cohort members were subsequently invited to participate in seven data collection phases.

-11

In the present analysis, we utilised data from phases three (1991 1993) and seven (2002 2004) the only occasions when CRP and ---IL-6 were measured. After excluding participants who reported having had a cold or flu in the 14 days preceding blood collection, or had a CRP over 10 mg/L at either phase three or seven, a total of 3522 participants (2496 men and 1026 women) had values on inflammatory markers and anthropometric measures from both phases three and seven; this being our analytical sample (Supplemental table S1). The mean follow-up time between phase three and phase seven was 11.3 years (range 9.5 12.9 years).

-

The Whitehall II study has approval from the University College London Medical School Committee on the Ethics of Human

Research. Informed consent was gained from all participants.

Anthropometry

Data on weight, height and waist circumference were collected by trained study nurses at clinical screening sessions both at phase three and seven. Weight was measured to the nearest 0.1 kg. Height was measured with participants standing erect with their head in the Frankfort plane, and was recorded to the nearest 0.1 cm. BMI was calculated as weight (kg)/height (m) . Waist circumference was 2 measured halfway between the costal margin and the iliac crest; measurements were recorded to the nearest 0.1 cm.

Inflammatory markers

Fasting serum was collected between 0800 and 1300 hours and was stored at 70 C. Samples from phases three and seven were ånalysed at the same time. CRP was measured using a high-sensitivity immunonephelometric assay in a BN ProSpec nephelometer (Dade Behring, Milton Keynes, UK). IL-6 was measured using a high-sensitivity enzyme-linked immunosorbent assay (R&D Systems, Oxford, UK). Values lower than the detection limit (0.154 mg/L for CRP and 0.08 pg/mL for IL-6) were assigned a value equal to half the detection limit.

Covariates

Covariates included in the present study, measured at baseline were: age, ethnicity (White, South Asian, Black, other), socioeconomic position based on civil service employment grade (low, intermediate, high), smoking (never smoker, ex-smoker, current smoker), and leisure-time physical activity (none or mild physical activity vs. moderate activity at least 3 times or vigorous activity at least 1 2 times per week). Coronary heart disease (CHD) was defined by using the Multinational Monitoring of trends and determinants in cardiovascular disease (MONICA) criteria ( ). Diabetes was defined by a fasting glucose of 7 0 mmol/L or more; a 2-h postload glucose of 11 1 12

• •
mmol/L or more; self-report of doctor diagnosis; or use of diabetes medication ( ). 13

Statistical analysis

The distributions of CRP and IL-6 at baseline and follow-up were skewed, so the geometric means and standard deviations of logged values are presented for these variables. The distributions for the 6( ) changes in CRP and IL-6 between phases three and se ven were 13 not skewed, so the arithmetic means and standard deviations are presented for these variables. -

RESULTS

In the analytical sample, mean age at baseline (phase three) was 49.4 (SD 6.0) years; men constituted 71 of the sample; and most = % study participants were white ( ). Mean weight increased in this sample by 3.0 (SD 5.8) kg over 11 years of follow-up. The Table 1 = increase in BMI and waist circumference was 1.4 (SD 2.0) kg/m and 7.7 (SD 7.4) cm for all participants. On average, CRP increased by were noted in the associations between change in weight, BMI or waist circumference, and IL-6.

In analyses stratified by baseline BMI categories, the relationship between change in the anthropometric measures and change in CRP was stronger in those overweight and obese at baseline, as compared with those who were lean or normal weight at baseline ( ). Table 2 This observation was unchanged irrespective of whether we used change in weight (P 0.002), BMI between change in waist circumference and baseline BMI category was statistically significant (P 0.002).

=

To test the robustness of our findings, we repeated the main multivariable linear regression analyses of weight change and change in inflammatory markers after excluding participants who had any CHD up to and including phase three (n 266), type 2 diabetes at phase = three (n 78), as well as those who reported taking lipid-lowering drugs any time during follow-up (n 420)(note that lipid lowering

= =
medication may affect levels of inflammation). In these analyses our conclusions were essentially unchanged (results not shown but available upon request).

DISCUSSION

We found that CRP and IL-6 levels increased over an 11-year follow-up period in a middle-aged population. The size of this increase in inflammatory markers was associated with the change in weight, BMI or waist circumference. These associations were more pronounced among those who were overweight or obese at baseline relative to normal weight persons.

The strengths of this study are use of multiple indicators of both adiposity and inflammation, and a large sample size. With little difference between the characteristics of participants at baseline and those included in this analytical sample, there was no strong evidence of selection bias (supplementary table S1). However, our study uses observational data, which can never provide complete inferences about causality.

In a UK community-based cohort study of 1222 participants followed over 9 years, an average increase in CRP of 0.09 mg/L per kg increase in weight was observed, which is consistent with our findings ( ). Our findings also accord with a French study of 1099 9 middle-aged participants which found significant correlations between change in BMI, waist-hip ratio and waist girth, and change in CRP, haptoglobin, orosomucoid, and between change in waist-hip ratio and waist girth and change in white blood cells, over a period of 5 years ( ). 10

To our knowledge, this is the first study to report that those who are obese may have a greater absolute increase in CRP per unit increase in subsequent weight, BMI and waist circumference, than those who are normal weight. This observed effect is plausible. Obese persons are known to have enlarged adipocytes ( , ), and adipose tissue in obese persons shows higher expression of several 8 15 tissue macrophage type, from the anti-inflammatory alternatively-activated type to the pro inflammatory classically-activated type ( ). 18

Furthermore, a higher incidence of pro-inflammatory co morbidities, such as osteoarthritis, has been observed among obese persons ( ). 19 Additional research is needed to examine whether these and other mechanisms may be responsible for the interaction effect observed in our study.

TABLE 1

Baseline characteristics and changes between baseline at follow-up ( ) in obesity measures and inflammatory markers of the study sample with full measurements on inflammatory markers and Δ anthropometry at both baseline (phase three, 1991 1993) and follow-up (phase seven, 2002 2004), the Whitehall II study.

--

All

  Multivariable linear regression models were used to evaluate the relationships of change in weight, BMI and waist circumference (exposure variables) with changes in CRP and IL-6 (outcome variables). Throughout these analyses we simultaneously took into account the effect of a range of potential confounding factors. To examine whether the association was dependent on baseline obesity status, we stratified analyses by the following baseline categories of BMI: underweight/normal weight (BMI<25.0); overweight (25.0 BMI<30.0);≤ 3 7and obese (BMI 30.0) ( ). Statistical interaction was evaluated by including product terms between baseline BMI category and the ≥ 14 exposure variable of interest in the respective models (i.e., weight change, BMI change or change in waist during the follow-up). These analyses were performed separately for the two outcome measures.All analyses were performed using SAS version 9.1 (SAS Institute Inc., Cary, NC, USA, 2002 2003).

= 2 =

 2 0.53 (SD 1.81) mg/L and IL-6 by 0.45 (SD 1.79) pg/mL.= =The crude changes in mean CRP and IL-6 during the follow-up by quartiles of change in weight, BMI, and waist circumference, are shown in Supplemental FigureS1. These data suggest a linear association between increase in inflammatory markers and increase in weight, BMI and waist circumference. Linear regression showed a one kg increase in body weight over the 11 years of follow-up to be associated with a mean increase in CRP of 0.08 (95 confidence interval (CI) 0.07 0.09) mg/L, while the corresponding increase in IL-6 % was 0.04 (95 CI 0.03 0.05) pg/mL ( ). The relationship between change in weight and BMI in relation to change in CRP was % -

  Baseline BMI category x change in weight interaction ' = (P 0.003) or waist circumference (P <0.0001) as the exposure Baseline BMI category x change in BMI interaction ' ' = Baseline BMI category x change in waist interaction ' ' variable in the models. For example, the average change in CRP was 0.31 (95 CI 0.22 0.40) mg/L per unit change in BMI in the obese % group, while the corresponding figure in the overweight and lean/normal weight groups were 0.22 (95 CI 0.17 0.27) mg/L and 0.18 (95 % -CI 0.14 0.22) mg/L, respectively. A similar pattern, albeit weaker, was observed in relation to change in IL-6, but only the interaction % -

production in obese adipose tissue ( , ). Mouse models have demonstrated that diet-induced obesity leads to a shift in the adipose

Ackowledgements:

We thank all participating men and women in the Whitehall II Study; all participating Civil Service departments and their welfare, personnel, and establishment officers; the Occupational Health and Safety Agency; and the Council of Civil Service Unions. The Whitehall II Study team comprises research scientists, statisticians, study coordinators, nurses, data managers, administrative assistants and data entry staff, who make the study possible.

Whitehall II study is supported by Medical Research Council; British Heart Foundation; Wellcome Trust; Health and Safety Executive;

Foundation. MiK is supported by the Academy of Finland.

Footnotes:

None of the authors declared conflicts of interest.

Disclosure statement:

This is an un-copyedited author manuscript copyrighted by The Endocrine Society. This may not be duplicated or reproduced, other that for personal use or within the rule of Fair Use of Copyrighted Materials (section 107, Title 17, U.S. Code) without permission of the copyright


Maria Melchior
email: maria.melchior@inserm.fr

Marie Choquet
email: marie.choquet@inserm.fr

Yann Le Strat
email: yann.lestrat@inserm.fr

Christine Hassler
email: christine.hassler@inserm.fr

Philip Gorwood
email: p.gorwood@ch-sainte-anne.fr

Parental alcohol dependence, socioeconomic disadvantage and alcohol and cannabis dependence among young adults in the community

Keywords: alcohol dependence, cannabis, family history, socioeconomic position, young adults, epidemiology

HAL is

Introduction

In industrialized countries, alcohol and illicit drug-related problems affect 5% of the population (2;19;27) and cost 1-3% of gross domestic product (GDP) (13;17;20;21;23), constituting a public health and societal challenge. Individuals with a family history of substance dependence are at high risk (1;28;35); however, for reasons that are unknown, intergenerational transmission is heterogeneous, that is not all individuals with family history of substance dependence have substance-related problems themselves [START_REF] Vanvoorst | Are relations between parental history of alcohol problems and changes in drinking moderated by positive expectancies? Alcoholism[END_REF].

Using data from a community-based study, we test the hypothesis that intergenerational transmission of substance dependence varies with exposure to socioeconomic disadvantage, previously shown to predict alcohol and drug use disorders (15;16;34). Parental substance-related problems are captured by history of alcohol dependence, the most frequent form of substance dependence (2;27) and a common correlate of illicit drug dependence [START_REF] Kessler | Lifetime prevalence and age-of-onset distributions of DSM-IV disorders in the National Comorbidity Survey Replication[END_REF]. Analyses control for factors associated both with family socioeconomic and substance dependence characteristics and with youths' risk of alcohol and cannabis dependence, including immigrant status, history of sexual abuse [START_REF] Fergusson | The developmental antecedents of illicit drug use: evidence from a 25-year longitudinal study[END_REF], and mental health difficulties such as conduct disorder and depression (18;34).

Material and Methods

The SAGE (Susceptibility Addiction Gene Environment) study examines factors associated with psychiatric disorders and alcohol and drug use among young adults in France (n= 3,056 youths in postgraduate training in North-Eastern France, March-April 2007, mean age: 20, sd: 1.4, 60.1% male, 79% response rate) [START_REF] Strat | First positive reactions to cannabis constitute a priority risk factor for cannabis dependence[END_REF]. The study received approval from France's national body supervising ethical data collection (CNIL).

Alcohol dependence was assessed using the well-validated 10-item WHO AUDIT questionnaire (7;10) and defined as a score of >=13 (males) or >= 12 (females). Cannabis dependence was assessed using a 10-item questionnaire derived from the Diagnostic Interview for Genetic Studies (DGIS) [START_REF] Nurnberger | Diagnostic interview for genetic studies. Rationale, unique features, and training. NIMH Genetics Initiative[END_REF] and defined as >=3 symptoms among the 7 DSM IV criteria. This questionnaire is concurrent with a semi-structured clinical interview [START_REF] Strat | First positive reactions to cannabis constitute a priority risk factor for cannabis dependence[END_REF].

Parental (mother's and father's) alcohol dependence was ascertained using the 13-item abridged version of the SMAST (Short Michigan Alcoholism Screening Test), which has good reliability and validity [START_REF] Crews | Using adapted short MASTs for assessing parental alcoholism: reliability and validity[END_REF] and defined as >=3 symptoms [START_REF] Selzer | A self-administered Short Michigan Alcoholism Screening Test (SMAST)[END_REF]. Participants with one or two parents with alcohol dependence had positive parental history (85.4%: fathers, 12.6%: both parents).

Following previous research, socioeconomic disadvantage was ascertained as a composite of parental employment status and educational level (4;46) and defined as parental unemployment or <secondary education.

Lifetime history of sexual abuse was ascertained using one item on attempted or completed unwanted sexual intercourse. Adolescent conduct disorder was measured using a self-report questionnaire derived from the DIGS (37;41) and defined as >=3 of the 15 DSM-IV criteria before age 15 [START_REF]Diagnostic and statistical manual of mental disorders[END_REF]. Current depression was assessed using the 10-item Adolescent Depression Rating Scale (ADRS) [START_REF] Revah-Levy | The Adolescent Depression Rating Scale (ADRS): a validation study[END_REF] and defined as >=3 positive symptoms. Immigrant status was defined as at least one parent born not in France (58.3% in this group had two parents born abroad). Immigrant participants primarily came from North Africa (Algeria, Morocco, Tunisia: 38.6%) or Europe (37.0%), but due to insufficient statistical power we could not distinguish these subgroups in the analyses.

To test the hypothesis that socioeconomic disadvantage exacerbates the intergenerational transmission of alcohol and cannabis dependence, we created a combined measure of youths' familial risk (low risk: no parental alcohol dependence n=2,476; intermediate risk: parental alcohol dependence but no socioeconomic disadvantage n=339; high risk: parental alcohol dependence and socioeconomic disadvantage n=69). First, analyses were controlled for sex. Next, we additionally controlled for all covariates. Data were analyzed using logistic regression models in the SAS statistical software package (SAS V9, Carey, North Caroline).

Results

The prevalence rates of alcohol and cannabis dependence among study youths were 5.8% and 7.3% (correlation coefficient: 0.26, p-value <0.0001). Youths whose parents had a history of alcohol dependence were more likely to have alcohol or cannabis dependence themselves (fullyadjusted ORs: alcohol dependence: 1.66, 95% 1.20-2.48, cannabis dependence: 1.48, 95% CI 1.02-2.15). Socioeconomic disadvantage was not independently associated with either study outcome (fully-adjusted ORs: alcohol dependence: 0.68, 95% CI 0.38-1.23, cannabis dependence: 0.92, 95% CI 0.57-1.47). Table 1 shows youths' characteristics in relation to alcohol and cannabis dependence. As depicted in Figure 1, we observed a gradient-like association between familial risk and youths' alcohol or cannabis dependence. In sex-adjusted regression analyses (Table 2), a family history of alcohol dependence was associated with alcohol or cannabis dependence, especially in the presence of socioeconomic disadvantage (ORs: 2.64, 95% CI 1.22-5.70 and 2.27, 95% CI 1.10-4.70). In fully-adjusted analyses (Table 2), the ORs were decreased, however the likelihood of cannabis dependence in relation to high familial risk remained elevated and statistically significant. Formal tests of additive statistical interactions [START_REF] Rothman | Modern epidemiology[END_REF], conducted using a method proposed by Andersson et al. [START_REF] Andersson | Calculating measures of biological interaction[END_REF], showed that the Relative Excess Risk due to Interaction (RERI) between parental alcohol dependence and socioeconomic disadvantage was much different from 0 but did not reach statistical significance, probably due to a small number of cases (RERI: for alcohol dependence: 1.03, 95% CI -0.99-3.06, for cannabis dependence: 0.48, 95% CI -1.33-2.31).

Discussion

In a community-based sample, young adults with parental history of alcohol dependence were disproportionately likely to have alcohol or cannabis dependence, particularly if they experienced socioeconomic disadvantage. This social disparity in the intergenerational transmission of substance dependence partly reflected risk factors of substance dependence including history of sexual abuse and mental health difficulties. Addictive behaviors probably result from a combination of heritable and environmental risks and youths who cumulate both risks may be especially vulnerable.

The study's limitations are: 1) a sample of postgraduate students; 2) a cross-sectional design;

3) participant reports of parental alcohol dependence. In France, approximately 50% of youths achieve postgraduate training and our sample may not include individuals who experience severe socioeconomic hardship (22;38). Nevertheless, SAGE study participants represent diverse socioeconomic backgrounds and have levels of alcohol and cannabis disorders similar to national samples [START_REF] Legleye | Usages de drogues des étudiants, chômeurs et actifs de 18-25 ans[END_REF]; thus our results should apply to other settings. Youths' substance-related problems are unlikely to predict parental alcohol dependence or socioeconomic disadvantage; however participants' reports may be influenced by substance use. Reassuringly, informant reports of substance dependence are highly specific (>90%) [START_REF] Milne | The validity of the family history screen for assessing family history of mental disorders[END_REF] and our results are concurrent with those of prospective studies [START_REF] Agrawal | Are there genetic influences on addiction: evidence from family, adoption and twin studies[END_REF]. Still, additional studies using multiple assessments of family history of alcohol and drug dependence would be useful.

Our study's main strength is a large sample of community-based young adults. The period of transition between adolescence and adulthood is key in the emergence of long-term patterns of substance abuse [START_REF] Chen | The natural history of drug use from adolescence to the mid-thirties in a general population sample[END_REF], yet young adults are difficult to include in epidemiological studies and data on this demographic group are few. Respectively 5.8% and 7.3% of study participants had alcohol and cannabis dependence, adding to evidence that a non-negligible proportion of youths does not desist from problematic substance use upon entering adulthood (33;39;40). Youths' rates of substance use have increased in recent years (8;12) and better understanding of lifelong risk trajectories is needed.

Our study included immigrant youths, who are rarely studied in France [START_REF] Lert | Functional limitations and overweight among migrants in the Histoire de Vie study (Insee, 2003)[END_REF]. Compared to non-immigrants, this group had lower rates of alcohol dependence but comparable rates of cannabis dependence, implying specificity in substance use in relation to immigration status [START_REF] Jayakody | Illicit and traditional drug use among ethnic minority adolescents in East London[END_REF]. Immigrants are disproportionately exposed to socioeconomic disadvantage, which justifies close monitoring of their health in a way that accounts for community of origin, conditions of migration, and acculturation in the host country [START_REF] Venema | Health of migrants and migrant health policy, The Netherlands as an example[END_REF].

In our study, the cumulative effect of parental alcohol dependence and socioeconomic disadvantage on youths' substance dependence was partly explained by risk factors such as history of sexual abuse, adolescent conduct disorder, and depression. The mechanisms linking familial context to these risk factors may include poor parenting (6;9;24;29) as well as high levels of family stressors such as financial difficulties, divorce or family conflict (9;42). Additionally, the intergenerational transmission of substance dependence may also reflect genetic influences, the expression of which may be enhanced in detrimental environmental conditions [START_REF] Agrawal | Are there genetic influences on addiction: evidence from family, adoption and twin studies[END_REF]. It is also possible that parental alcohol dependence leads to family socioeconomic disadvantage; however, data needed to test test this hypothesis were not available in our study. In the future, the moderation of genetic risk of substance dependence by socioeconomic disadvantage should be tested.

Conclusion

Parents appear to have a lasting influence on their offspring's substance use well into adulthood. Improvements in families' socioeconomic conditions could reduce the transmission of substance dependence to the next generation.

1

  Parental and socioeconomic risk: low=no parental alcohol dependence, intermediate=parental alcohol dependence, no socioeconomic disadvantage; high: parental alcohol dependence and socioeconomic disadvantage.

Figure 1

 1 Figure 1 Prevalence of alcohol and cannabis dependence among young adults according to familial risk level: the SAGE study (%).A. Alcohol dependence

Table 1 . Demographic, family and mental health characteristics of young adults in relation to alcohol and cannabis dependence: the SAGE study (n=2,884).

 1

	Prevalence (%)	Alcohol dependence	Cannabis dependence
		OR (95% CI)	OR (95% CI)

Table 2 . Familial risk level in relation to alcohol and cannabis dependence among young adults: the SAGE study (logistic regression analysis OR, 95% CI) Alcohol dependence Cannabis dependence Sex-adjusted models

 2

	(n=2884)
	Familial risk : Low	1.0	1.0
	Intermediate	2.09 (1.39-3.15)	1.77 (1.20-2.60)
	High	2.64 (1.22-5.70)	2.27 (1.10-4.70)
	Sex: Female	1.0	1.0
	Male	3.46 (2.31-5.18)	2.09 (1.57-2.89)
	Fully-adjusted models (n=2851)
	Familial risk 1 : Low	1.0	1.0
	Intermediate	1.61 (1.04-2.48)	1.52 (1.02-2.28)
	High	2.06 (0.88-4.82)	1.43 (1.62-3.32)
	Sex: Female	1.0
	Male	3.29 (2.10-5.12)
	Immigrant status: Immigrant	1.0
	Non-immigrant	2.19 (1.36-3.53)
	History of sexual abuse: Absent	1.0
	Present	2.42 (1.10-5.31)
	Adolescent conduct disorder: Absent	1.0
	Present	3.80 (2.65-5.46)
	Depression: Absent	1.0
	Present	1.83 (1.26-2.64)

Family risk: low=no parental alcohol dependence, intermediate=parental alcohol dependence, no socioeconomic disadvantage; high: parental alcohol dependence and socioeconomic disadvantage.

Acknowledgments:

The authors thank the schools and youths who participated in the SAGE study.

They also wish to express their gratitude to France Lert for her comments on a previous version of the manuscript and Hermann Nabi for methodological help. The SAGE study was funded by the Institut de Recherche sur les Boissons (IREB) and the Mission Interministérielle de Lutte contre la Drogue et la Toxicomanie (MILDT) in France. Maria Melchior is the recipient of a 'Young Researcher' award from France's Agence Nationale de la Recherche (ANR).

Conflict of interest: none.


Archana Singh-Manoux
email: <archana.singh-manoux@inserm.fr

Commentary: is it time to redefine cognitive epidemiology?

Keywords: Achievement, Age Factors, Data Collection, methods, Educational Status, Epidemiologic Methods, Humans, Intelligence, Research

In this issue of the IJE, Deary and Johnson discuss the association between education and intelligence and argue that their use in

when in fact the confounder or the common cause in the analyses ought to be intelligence. Epidemiologists generally agree that the association between an exposure and a health outcome contains two elements; the true causal effect of the exposure on health and the effect due to a shared common cause, the confounder ; . Analyses are adjusted for the confounder (common cause) in order to remove [2 3 ] its spurious effect on the association between an exposure and the health outcome in question. If intelligence is the real confounder of these associations then, in the absence of a direct measure of intelligence, is it correct to adjust for education? The answer is yes. Hernan and colleagues use directed acyclic graphs to show that when the hypothesised causal confounder (in this case intelligence) is not measured, adjusting for the measured variable (education) is correct as it acts as a surrogate confounder and adjustment for it removes most of the effect of the common cause (intelligence) from the association between the exposure and the outcome .

[3 ]

Two-way relationship between education and intelligence

It is possible that there is a two-way association between education and intelligence; although I am not sure that the two criteria, temporality and heritability, discussed by Deary and Johnson are sufficient to enable us to determine the causal direction of the association . Temporality alone is not a sufficient criterion, and the heritability estimate of most human traits at around 50 does not help us [START_REF] Rothman | Causation and causal inference in epidemiology[END_REF] [4 ] %

resolve issues of causality. The author s suggest that the temporal cascade between intelligence and education will be clearer when

' "
repeated measures of each are available in longitudinal studies. However, at present statistical methods that allow causality to be " established in the presence of reciprocal effects do not exist.

Natural experiments tell us something about causality. A large body of research supports the existence of the Flynn effect, the trend for mean IQ scores to increase between generations, leading to much discussion of the explanations (improved nutrition, better care of children, better education, increasing familiarity with tests, heterosis, etc.) for the malleability of intelligence scores . The increase in [START_REF] Dickens | Heritability estimates versus large environmental effects: the IQ paradox resolved[END_REF] test scores can be dramatic as evident in , which shows three repeat assessments on a test of reasoning (AH4-I), also a test of Figure 1 [6 ] fluid intelligence with a large loading on the general intelligence factor, in men and women of the Whitehall II study, a cohort of British civil servants . Gender differences in childcare, nutrition, familiarity with multi-choice questions or heterosis appear unlikely. However, [START_REF] Marmot | Cohort Profile: the Whitehall II study[END_REF] we know that education levels in these women have improved; 22.6 of the oldest women (birth years 1930 1934) completed high school % education compared to 66.5 of the youngest women (birth years 1950 52) in this cohort. Does this prove that education causes % intelligence? No, but how important is it for epidemiology to establish this causal direction? Epidemiology is a practical discipline concerned with identifying modifiable determinants to improve health at the population level. While intelligence lies beyond the reach of policy makers, education is a much easier target. Many middle-and low-income countries set improving education as a major policy goal with a view to improving both the economic and health prospects of the population. Policies directed at education span the range from minimum school leaving age to increasing of the population with a university degree.

%

Defining cognitive epidemiology

The final point made by Deary and Johnson is that the direction of the association between education and intelligence is important for the newly emerging field of cognitive epidemiology , defined previously as the use of cognitive ability test scores as risk factors for " " "

human health and disease outcomes, including mortality . Given the association between intelligence and education, extensively

" [8 ]
discussed by Deary and Johnson, this definition of cognitive epidemiology puts it squarely in the domain of social epidemiology, a discipline concerned with the social distribution of determinants of health. Location in this broader church, rather than the micro-discipline of cognitive epidemiology, will avoid a narrow focus on intelligence that ignores its associations with markers of social position such as education, income, and occupation.

My concern is also with the cognitive in this definition of cognitive epidemiology. Cognitive function and intelligence are not " "

interchangeable concepts. For one, cognitive function is not a unitary concept like intelligence; it is composed of multiple aspects involved in the input, storage, processing and output of information Not all cognitive functions decline with age at the same rate. Some, like [START_REF] Lezak | Neuropsychological assessment[END_REF] vocabulary, show little decline. Others, such as measures of executive function, are more likely to be affected by cardiovascular risk factors. Decline in memory, besides its importance for Alzheimer s disease, might have a stronger association with mortality ' [START_REF] Sabia | Does cognition predict mortality in midlife? Results from the Whitehall II cohort study[END_REF] Furthermore, most epidemiologists study cognitive function as an outcome rather than an exposure, making the definition proposed by Deary and colleagues rather restrictive.

Continuing increases in life expectancy imply fundamental changes to the population structure. Impaired cognitive status is one of the biggest challenges of the future due to its impact on both the individual and society. We now know that decline in multiple domains of cognitive function starts in midlife, evident in those under 50 in the Whitehall II study ( ). We also know that the inter-individual Figure 1 differences steadily increase with age, even though the determinants of this increasing heterogeneity are not entirely clear. For these [11 ] reasons, more than the impact of cognition on health, cognitive epidemiology should principally be the study of the pathophysiology of cognitive decline, involving etiologic models to identify the biological, behavioural, psychological, and social risk factors. Although age itself is not modifiable, given its heterogeneous impact on cognitive decline it is time to deconstruct this effect in order to meet the challenge of ageing populations. Epidemiology s contribution to medicine is the study of the causes, distribution, and control of disease in ' populations, cognitive epidemiology need not be any different.

Ackowledgements:

I would like to thank Aline Dugravot, Jane Ferrie, Herman Nabi, S verine Sabia and Silvia Stringhini for comments on a previous draft of this é paper. I acknowledge support from the European Science Foundation via a European Young Investigator Award and the National Institute

" "
on Aging, NIH (R01AG013196; R01AG034454).

References:

Figure 1

 1 Figure 1 Decline (and 95 CI) in reasoning (Alice-Heim 4-I) estimated from 3 repeat measures over 10 years as a function of year of birth in 5247 % men and 2233 women in the Whitehall II study.


Stéphane Guindon
email: guindon@stat.auckland.ac.nz

Frédéric Delsuc

Jean-François Dufayard

Olivier Gascuel

Stephane Guindon

Estimating maximum likelihood phylogenies with PhyML

Keywords: DNA and protein sequences, molecular evolution, sequence comparisons, phylogenetics, statistics, maximum likelihood, Markov models, algorithms, software, PhyML

Our understanding of the origins, the functions and/or the structures of biological sequences strongly depends on our ability to decipher the mechanisms of molecular evolution. These complex processes can be described through the comparison of homologous sequences in a phylogenetic framework. Moreover, phylogenetic inference provides sound statistical tools to exhibit the main features of molecular evolution from the analysis of actual sequences. This chapter focuses on phylogenetic tree estimation under the maximum likelihood (ML) principle. Phylogenies inferred under this probabilistic criterion are usually reliable and important biological hypotheses can be tested through the comparison of different models. Estimating ML phylogenies is computationally demanding though and careful examination of the results is warranted. This chapter focuses on PhyML, a software that implements recent ML phylogenetic methods and algorithms. We illustrate the strengths and pitfalls of this program through the analysis of a real data set.

1 Introduction.

In statistics, models are mathematical objects designed to approximate the processes that generated the data at hand. These models can be more or less complex, depending on the type of data and the available knowledge on the process that generated them. Each model has parameters which values need to be estimated from the data. Least-squares, maximum a posteriori estimation (MAP) or maximum likelihood (ML) are the main statistical frameworks suitable for this task. The least-squares criterion has been widely used in phylogenetics in order to build trees from matrices of pairwise distances between sequences. The last two criteria, ML and MAP, both rely on the probability that the data were generated according to the selected model. This conditional probability is the likelihood of the model. The next section gives a short description of the type of models that are used for phylogenetic inference.

Mathematical description of sequence evolution.

In molecular evolution, we commonly assume that homologous sequences evolve along a bifurcating tree or phylogeny. The topology of this tree describes the different clades or groups of taxa. The tree topology is generally considered to be the most important parameter of the whole phylogenetic model. The second parameter is the set of branch lengths on a given topology. The length of each branch on this topology represents an amount of evolution which corresponds to an expected number of nucleotide, codon or amino-acid substitutions. The last component of the model is a mathematical description of the process that generates substitutions during the course of evolution. Several assumptions are made about this process. We first consider that the sites of the alignment, i.e., its columns, evolve independently (see Note 1) and under the same phylogeny. We also assume that the substitution process is the same at the different sites of the alignment. This process is modelled with Markov chains. The two main components of a Markov model of substitution are (1) a symmetrical matrix that describes the relative speed at which the different substitution events occur (e.g., transition and transversion rates are generally distinct), and (2) the frequencies of the different states to be considered (i.e., nucleotides, codons or amino-acids). The symmetrical matrix is often referred to as the exchangeability matrix. If combined with the vector of state frequencies, the resulting matrix corresponds to the generator of the Markov process and the vector of state frequencies is the stationary distribution.

A difficult optimisation problem.

ML tree estimation consists of finding the phylogenetic model, i.e., the tree topology, branch lengths and parameters of the Markov model of substitution, that maximises the likelihood. Calculating the likelihood of a given phylogenetic model can be done efficiently using Felsenstein's pruning algorithm [START_REF] Felsenstein | Evolutionary trees from DNA sequences: a maximum likelihood approach[END_REF]. Unfortunately, finding the ML model is a difficult problem. The difficulty mostly comes from the very nature of the model itself. Indeed, while branch lengths and parameters of the substitution model are continuous variables, the tree topology is a discrete parameter. Hence, it is not surprising that, for most data sets, the likelihood function defines a rugged landscape with multiple peaks (see (2) however). Searching for the ML phylogenetic model in such conditions therefore rely on sophisticated optimisation methods, which combine both discrete and continuous optimisation procedures.

These methods are heuristics. As opposed to exact algorithms, heuristics do not guarantee to find the best (i.e., ML) solution. Despite this, simulation studies [START_REF] Huelsenbeck | Success of phylogenetic methods in the four-taxon case[END_REF][START_REF] Swofford | Phylogenetic inference[END_REF] suggest that the heuristics designed for phylogenetic estimation perform well overall. The most efficient methods are those that provide the best trade-off between their ability to maximise the likelihood function and the time spent to achieve this. Several methods/programs use deterministic approaches: given an initial non-optimal solution, the heuristic always follows the same path of intermediate solutions to reach the estimated ML one. As a consequence, given a particular input (i.e., data and model settings), these methods always produce the same output (i.e., the estimated ML phylogenetic model). Other approaches implement non-deterministic heuristics. These methods do not always follow the same path of intermediate solutions to reach the estimated ML tree. Hence, the same input potentially leads to distinct outputs (see Note 2).

Searching through the space of tree topologies.

Another important distinction between tree building methods lies in the type of operations (or moves) that are used to explore the space of tree topologies. The three principal moves are: nearest neighbour interchange (NNI), subtree pruning and regrafting (SPR) and tree bisection and reconnection (TBR) (see [START_REF] Swofford | Phylogenetic inference[END_REF]). Each of these moves permits an exhaustive exploration of the space of tree topologies. However, the neighbourhood of trees defined by applying TBR moves to any given topology is much wider than the neighbourhood defined by SPR moves, which is itself larger than the neighbourhood defined by NNI moves.

Therefore, TBRs are more efficient than SPRs or NNIs in jumping across very distinct tree topologies in just one step. The same also holds for SPR vs. NNI moves. These differences result in various abilities to escape local maxima of the likelihood function. Consider a suboptimal peak of the likelihood surface and the tree topology found at this peak. Every potential NNI move will only allow to reach similar topologies. Therefore, such moves will sometimes fail to reach a higher peak of the likelihood surface. Hence, SPR operations are more efficient than NNIs in finding the highest peaks of the likelihood surface and TBR moves are more efficient than SPR ones. However, the large neighbourhood of trees defined by TBR compared to SPR or NNI also generally implies greater run times.

Basically, many non-optimal solutions are evaluated when using TBR moves as compared to SPRs or NNIs. Hence, the best methods with respect to likelihood maximisation also tend to be the slowest ones.

1.4 PhyML v3.0: new features.

PhyML (5) is a software that estimates ML phylogenies from alignments of nucleotide or amino acid sequences. It provides a wide range of options that were designed to facilitate standard phylogenetic analyses. This chapter focuses on PhyML v3.0. This version of the program provides two important advances compared to the previous releases. Indeed, PhyML v3.0 now proposes three different options to search across the space of phylogenetic tree topologies. It also implements a new method that evaluates branch supports. These new options are presented below.

Escaping local maxima.

PhyML originally relies on a deterministic heuristic based on NNI moves. Traditional greedy approaches first evaluate the gain of likelihood brought by every possible NNI applied to the current topology. Only the best move is then used to improve the phylogeny at each step of the tree building process. Hence, most of the calculations evaluate moves that will not be used subsequently. To avoid such waste of potentially useful information, PhyML applies several 'good' NNI moves simultaneously. Moreover, these local and simultaneous changes of the tree topology are accompanied by adjustments of every branch length in the phylogeny (see Note 3). Hence, most of the calculations that are performed during one step of the algorithm are actually used to improve the tree. This is what makes PhyML faster than popular algorithms such a fastDNAml [START_REF] Olsen | fastDNAml: a tool for construction of phylogenetic trees of DNA sequences using maximum likelihood[END_REF], and simulation results [START_REF] Guindon | A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood[END_REF] demonstrated the accuracy of this approach. PhyML also outperforms the standard NNIbased greedy searching algorithms in terms of maximisation of the likelihood function.

However, the analysis of real data and the comparison with SPR-based algorithms shows that PhyML occasionally gets trapped in local maxima. This shortcoming is more and more obvious as the number of sequences to analyse gets large (i.e., > 50-100). Hence, while the phylogenetic models estimated with PhyML are generally good in terms of likelihoods, models with greater likelihoods can often be found. This is the reason why the release 3.0 of PhyML proposes new SPR-based tree searching options. The methods implemented in PhyML v3.0 are inspired by Hordijk and Gascuel [START_REF] Hordijk | Improving the efficiency of SPR moves in phylogenetic tree search methods based on maximum likelihood[END_REF] work on the topic. Their approach essentially relies on using a fast distance-based method to filter out SPR moves that do not increase the likelihood of the current phylogenetic model. Hence, the likelihood function is only evaluated for the most promising moves. This strategy proves to be efficient in finding trees with high likelihoods. While the computational burden involved with SPRs is heavier than with simultaneous NNIs, this new approach is clearly less prone to be stuck in local maxima of the likelihood function (see Note 4). PhyML 3.0 also proposes an intermediate option that includes both NNIs and SPRs.

Simultaneous NNIs are first applied following the original algorithm, until no additional improvement is found. A single round of SPR moves are then tested: each subtree is pruned, regrafted, filtered and only the most promising moves are actually evaluated. If one or more SPR moves increase the likelihood, the best one is applied to the current tree. A new round of simultaneous NNIs then starts off after this step. Simultaneous NNIs and SPRs therefore alternate until a maximum of the likelihood function is reached. This approach is generally faster than the SPR-only one but slower than the NNI-based heuristic. Also, while this strategy performs better than the NNI-only one in optimising the likelihood function, the SPR-only search often outperforms this mixed approach.

Fast tests for branch support.

PhyML v3.0 also provides users with a fast approximate likelihood ratio test (aLRT) for branches [START_REF] Anisimova | Approximate likelihood-ratio test for branches: a fast, accurate, and powerful alternative[END_REF], which proves to be a good alternative to the (time-consuming) bootstrap analysis. The aLRT is closely related to the conventional LRT, with the null hypothesis corresponding to the assumption that the tested branch has length 0. Standard LRT uses the test statistics 2(L 1 -L 0 ), where L 1 is the log-likelihood of the current tree, and L 0 the log-likelihood of the same tree, but with the branch of interest being collapsed. The aLRT approximates this test statistics in a slightly conservative but practical way as 2(L 1 -L 2 ), where L 2 corresponds to the second best NNI configuration around the branch of interest.

Such test is fast because the log-likelihood value L 2 is computed by optimising only over the branch of interest and the four adjacent branches, while other parameters are fixed at their optimal values corresponding to the best ML tree. Three branch supports computed from this aLRT statistics are available in PhyML v3.0: (1) the parametric branch support, computed from the χ 2 distribution (as usual with the LRT); (2) a non-parametric branch support based on a Shimodaira-Hasegawa-like procedure (9); (3) a combination of these two supports, that is, the minimum value of both. The default is to use SH-like branch supports.

The rational behind the aLRT clearly differs from non-parametric bootstrap, as detailed in [START_REF] Anisimova | Approximate likelihood-ratio test for branches: a fast, accurate, and powerful alternative[END_REF]. Basically, while aLRT values are derived from testing hypotheses, the bootstrap proportion is a repeatability measure; when the bootstrap proportion of a given clade is high, we are quite confident that this clade would be inferred again if another original data sample was available and analysed by the same tree-building method (which does not mean that the clade exists in the true tree). Also, computing aLRT values is much faster than getting bootstrap supports, as PhyML is run just once, while bootstrap requires launching PhyML 100 to 1,000 times. In fact, computing aLRT branch supports has a negligible computational cost in comparison with tree building. Note however that SH-like branch supports are non-parametric, just as are the bootstrap proportions. In fact, they often provide similar results as the bootstrap, pointing out the same poorly supported branches of the phylogeny, as illustrated in Fig. 1.

The aLRT assesses that the branch being studied provides a significant gain in likelihood, in comparison with the null hypothesis that involves collapsing that branch but leaving the rest of the tree topology identical. Thus, the aLRT does not account for other possible topologies that would be highly likely but quite different from the current topology. This implies that the aLRT performs well when the data contains a clear phylogenetic signal, but not as well in the opposite case, where it tends to give a (too) local view on the branch of interest and be liberal. Note also that parametric χ 2 branch supports are based on the assumption that the evolutionary model used to infer the trees is the correct one. In that respect, the aLRT parametric interpretation is close to Bayesian posteriors.

As the later, the resulting test is sometimes excessively liberal due to violations of the parametric assumptions.

Let us now focus on the practical aspects that go with ML phylogenetic model estimation using PhyML v3.0. The next section presents the inputs and outputs of the program, the different options and how to use them.

Program usage.

PhyML v3.0 has two different user interfaces. The default is to use the PHYLIP-like text interface (Fig. 2) by simply typing 'phyml' in a command-line window or by clicking on the PhyML icon (see Note 5). After entering the name of the input sequence file, the user goes through a list of sub-menus that allow her/him to set up the analysis. There are currently four distinct sub-menus:

1. Input Data: specify whether the input file contains amino acid or nucleotide sequences. What is the sequence format (see Section 2.1) and how many data sets should be analysed.

Substitution Model: selection of the Markov model of substitution (see Section 2.2).

3. Tree Searching: selection of the tree topology searching algorithm (see Section 2.4).

4. Branch Support: selection of the method that is used to measure branch support (see Section 2.3).

'+' and '-' keys are used to move forward and backward in the sub-menu list. Once the model parameters have been defined, typing 'Y' (or 'y') launches the calculations.

The meaning of some options may not be obvious to users that are not familiar with phylogenetics. In such situation, we strongly recommend to use the default options. As long as the format of the input sequence file is correctly specified (sub-menu Input data), the safest option for non-expert users is to use the default settings.

The alternative to the PHYLIP-like interface is the command line. Users that do not need to modify the default parameters can launch the program with the 'phyml -i your input sequence file name' command. The list of all command line arguments and how to use them is given in the 'Help' section which is displayed after entering the 'phyml help' command. Command lines are specially handy for launching PhyML v3.0 in batch mode. Note however that some options are only available through the PHYLIPlike interface. Hence, options that are not listed in the 'Help' section may be accessible through the interactive text interface.

Inputs / outputs.

PhyML reads data from standard text files, without the need for any particular file name extension. Alignments of DNA or protein sequences must be in PHYLIP sequential or interleaved format (see Fig. 3a). The first line of the input file contains the number of species and the number of characters, in free format, separated by blanks. One slight difference with PHYLIP format concerns sequence name lengths. While PHYLIP format limits this length to ten characters, PhyML can read up to hundred character long sequence names. Blanks and the symbols "(),:" are not allowed within sequence names because the NEWICK tree format makes special use of these symbols (Fig. 3b).

Another slight difference with PHYLIP format is that actual sequences must be separated from their names by at least one blank character. These sequences must not be longer than 10 6 amino acid or nucleotide characters and a given data set can have up to 4×10 3 of them. However, the size of the largest data set PhyML v3.0 can process depends on the amount of physical memory available. To avoid overflows, PhyML v3.0 pauses when the estimated amount of memory that needs to be allocated exceeds 250Mb. The user can then decide whether she/he wants to continue or cancel the analysis (see Note 6).

An input sequence file may also display more than a single data set. Each of these data sets must be in PHYLIP format and two successive alignments must be separated by an empty line. Processing multiple data sets requires to toggle the 'M' option in the Input Data sub-menu or use the '-n' command line option and enter the number of data sets to analyse. The multiple data set option can be used to process re-sampled data that were generated using a non-parametric procedure such as cross-validation or jackknife (a bootstrap option is already included in PhyML). This option is also useful in multiple gene studies, even if fitting the same substitution model to all data sets may not be suitable.

Gaps correspond to the '-' symbol. They are systematically treated as unknown characters "on the grounds that we don't know what would be there if something were there" (J. Felsenstein, PHYLIP main documentation). The likelihood at these sites is summed over all the possible states (i.e., nucleotides or amino acids) that could actually be observed at these particular positions. Note however that columns of the alignment that display only gaps or unknown characters are simply discarded because they do not carry any phylogenetic information (they are equally well explained by any model). PhyML v3.0 also handles ambiguous characters such as R for A or G (purines) and Y for C or T (pyrimidines). Tables 1 and2 give the list of valid characters/symbols and the corresponding nucleotides or amino acids. PhyML v3.0 can read one or several phylogenetic trees from an input file. This option is accessible through the Tree Searching sub menu or the '-u' argument from the command line. Input trees are generally used as initial ML estimates to be subsequently adjusted by the tree searching algorithm (see Section 2.4). This option is also helpful when one wants to evaluate the likelihood on a particular set of (possibly competing) phylogenetic trees. Trees should be in standard NEWICK format (Fig. 3b). They can be either rooted or unrooted and multifurcations are allowed. Taxa names must, of course, match the corresponding sequence names.

Single or multiple sequence data sets may be used in combination with single or multiple input trees. When the number of data sets is one (n D = 1) and there is only one input tree (n T = 1), then this tree is simply used as input for the single data set analysis. When n D = 1 and n T > 1, each input tree is used successively for the analysis of the single alignment. If n D > 1 and n T = 1, the same input tree is used for the analysis of each data set. The last combination is n D > 1 and n T > 1. In this situation, the i-th tree in the input tree file is used to analyse the i-th data set. Hence, n D and n T must be equal here.

Table 3 presents the list of files resulting from a PhyML v3.0 analysis. Basically, each output file name can be divided into three parts. The first part is the sequence file name, the second part corresponds to the extension ' phyml ' and the third part is related to the file content. When launched with the default options, PhyML v3.0 only creates two files: the tree file and the model parameter file. The estimated ML tree is in standard NEWICK format (Fig. 3b). The model parameters file, or statistics file, displays the ML estimates of the substitution model parameters, the likelihood of the ML phylogenetic model, and other important information concerning the settings of the analysis (e.g., type of data, name of the substitution model, starting tree (see Section 2.4), etc.). Two additional output files are created if bootstrap supports were evaluated. These files simply contain the ML trees and the substitution model parameters estimated from each bootstrap replicate. Such information can be used to estimate sampling errors around each parameter of the phylogenetic model. The best ML tree file is only created when the estimation of the phylogeny resulted from multiple random starting trees. It contains the tree with the highest likelihood that was found among all the trees estimated during the analysis (see Section 2.4).

Substitution models.

PhyML implements a wide range of substitution models: JC69 (10), K80 (11), F81 (1),

↔ C = A ↔ T = C ↔ G = G ↔ T .
This model corresponds to HKY85 (default) or K80 if the nucleotide frequencies are all set to 0.25. '010020' and '012345' correspond to TN93 and GTR models respectively. The digit string therefore defines groups of relative substitution rates. The initial rate within each group is set to 1.0, which corresponds to F81 (JC69 if the base frequencies are equal).

Users also have the opportunity to define their own initial rate values (this option is only available through the PHYLIP-like interface). These rates are then optimised afterwards (option 'O') or fixed to their initial values. The CUSTOM option can be used to implement all substitution models that are special cases of GTR. PhyML v3.0 also implements Ziheng Yang's discrete gamma model [START_REF] Yang | Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods[END_REF] to describe the variability of substitution rates across nucleotide or amino acids positions. Users can specify the number of substitution rate categories (Substitution Model sub-menu option 'C' or '-c' argument from the command line) and choose to estimate the gamma shape parameter from the data or fix its value a priori (option 'A' or '-a'). The program also handles invariable sites (option 'V' or '-v'). Here again, the value of this parameter can be estimated in the ML framework or fixed a priori by the user (see Note 8)

Branch support.

PhyML v3.0 proposes two main options to assess the support of the data for non-terminal branches in the phylogeny. The most popular approach relies on non-parametric bootstrap.

This option is available through the PHYLIP-like interface (Branch Support sub-menu) or the '-b' argument from the command line. Users only have to specify the number of replicates that will be generated to work out the bootstrap values. The ML output tree (see Section 2.1) will then display both branch lengths and bootstrap values. It is very important to keep in mind that bootstrap values are displayed on the ML tree estimated from the original data set. There is no consensus tree reconstruction involved here. Note however that PhyML also outputs every tree estimated from the re-sampled data sets.

These trees can then be used to build a consensus tree, using the program CONSENSE from the PHYLIP package for instance. However, there is no guarantee that this consensus tree is also the ML tree (but both generally have similar topologies). PhyML v3.0 also implements the approximate likelihood-ratio tests (aLRT) described in Section 1.6. Just like bootstrap support, aLRT options are available from the Branch Support sub-menu in the PHYLIP-like interface. The 'A' key is used to choose among four distinct options as aLRT supports can be assessed through different tests. The default is to use SH-like branch supports. It is also possible to test aLRT values against a χ 2 distribution or to use a combination of these last two options. Expert users also have the opportunity to retrieve the aLRT values themselves in order to examine the differences of likelihood among competing topologies. These four options are also available using the '-b' argument from the command line as explained in the 'Help' section.

Tree searching algorithms.

PhyML v3.0 implements three heuristics to explore the space of tree topologies (see Section 2). All methods take as input a starting tree that will be improved subsequently. The default is to let PhyML v3.0 build this starting tree using BioNJ [START_REF] Gascuel | BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data[END_REF]. However, users can also give their own input tree(s) in NEWICK format (see Section 2.1). Under the default settings, the starting tree is improved using simultaneous NNIs. PhyML v3.0 also proposes NNI+SPR and full-SPR heuristics that both search more thoroughly the space of tree topologies. Simultaneous NNIs, NNI+SPR and full-SPR options are available from the Tree Searching sub-menu or the '-s' argument from the command line.

Full-SPR searches can also be combined with multiple random starting trees (Tree Searching sub-menu, options 'S' and 'R'). Multiple random starting solutions are useful to evaluate the difficulty of the optimisation problem given the data at hand. Indeed, if the estimated solution (the ML phylogenetic model in our case) strongly depends on the value that is used to initiate the optimisation process (the starting tree), then the function to optimise is probably not smooth and local optima are likely to be commonplace. Hence, multiple random starting trees provide a diagnosis tool to evaluate the smoothness of the likelihood surface and the robustness of the inferred tree. The default is to use five random starting trees. Five ML phylogenetic models are then estimated and the best phylogenetic tree (i.e., the tree with the highest likelihood) is printed in the ' phyml best tree.txt' file (Table 3). When all these trees are identical, this tree is likely to be the ML solution.

Recommendations on program usage.

From the user perspective, the choice of the tree searching algorithm among those provided by PhyML v3.0 is probably the toughest one. The fastest option relies on local and simultaneous modifications of the phylogeny using NNI moves. More thorough explorations of the space of topologies are also available through the NNI+SPR and full-SPR options.

As these two classes of tree topology moves involve different amounts of computation, it is important to determine which option is the most suitable for the type of data set or analysis one wants to perform. Below is a list of recommendations for typical phylogenetic analyses.

1. Single data set, unlimited computing time. The best option here is probably to use a full-SPR search. If the focus is on estimating the relationships between species (see Section 3), it is a good idea to use more than one starting tree to decrease the chance of getting stuck in a local maximum of the likelihood function. Note however that NNI+SPR or NNI are also appropriate if the analysis does not mainly focus on estimating the evolutionary relationships between species (e.g. a tree is needed to estimate the parameters of codon-based models later on). Branch supports can be estimated using bootstrap and aLRT.

2. Single data set, restricted computing time. The three tree searching options can be used depending on the computing time available and the size of the data set. For small data sets (i.e., < 50 sequences), NNI will generally perform well provided that the phylogenetic signal is strong. It is relevant to estimate a first tree using NNI moves and examine the reconstructed phylogeny in order to have a rough idea of the strength of the phylogenetic signal (the presence of small internal branch lengths is generally considered as a sign of a weak phylogenetic signal, specially when sequences are short). For larger data sets (> 50 sequences), full-SPR search is recommended if there are good evidence of a lack of phylogenetic signal. Bootstrap analysis will generally involve large computational burdens. aLRT branch supports therefore provide an interesting alternative here.

3. Multiple data sets, unlimited computing time. Comparative genomic analyses sometimes rely on building phylogenies from the analysis of a large number of gene families. Here again, the NNI option is the most relevant if the focus is not on recovering the most accurate picture of the evolutionary relationships between species.

More time-consuming heuristics (NNI+SPR and full-SPR) should be used when the topology of the tree is an important parameter of the analysis (e.g., identification of horizontally transferred genes using phylogenetic tree comparisons). Internal branch support is generally not a crucial parameter of the multiple data set analyses. Using aLRT statistics is therefore the best choice.

4. Multiple data sets, limited computing time. The large amount of data to be processed in a limited time generally requires the use of the fastest tree searching and branch support estimation methods Hence, NNI and aLRT are generally the most appropriate here.

Another important point is the choice of the substitution model. While default options generally provide acceptable results, it is often warranted to perform a pre-analysis in order to identify the best-fit substitution model. This pre-analysis can be done using popular software such as Modeltest [START_REF] Posada | Modeltest: testing the model of DNA substitution[END_REF] or ProtTest [START_REF] Abascal | Prottest: selection of best-fit models of protein evolution[END_REF] for instance. These programs generally recommend the use of a discrete gamma distribution to model the substitution process as variability of rates among sites is a common feature of molecular evolution. The choice of the number of rate classes to use for this distribution is also an important one. While the default is set to four categories in PhyML v3.0, it is recommended to use larger number of classes if possible in order to best approximate the patterns of rate variation across sites [START_REF] Galtier | Markov-modulated Markov chains and the covarion process of molecular evolution[END_REF]. Note however that run times are directly proportional to the number of classes of the discrete gamma distribution. Here again, a pre-analysis with the simplest model should help the user to determine the number of rate classes that represents the best trade-off between computing time and fit of the model to the data.

Example.

As an illustration of using PhyML v3.0, we focus on reconstructing placental mammal phylogeny from the comparison of complete mitochondrial genomes. In such situation, taxonomic purposes are clearly more important than deciphering the evolutionary processes. This particular class of problems usually involves huge computational burdens.

The next sections present the data and goes through the different steps of the analysis.

The data.

Reconstructing the evolutionary history of major placental mammal lineages has been a long standing phylogenetic challenge. For long restricted to the study of morphological characters, placental mammal phylogenetics has widely benefited from the advent of molecular techniques giving access to a large number of informative characters from both mitochondrial and nuclear genomes. After a period of relative confusion mainly due to restricted taxon and gene sampling, phylogenetic analyses of placental mammal relationships based on complete mitochondrial genomes [START_REF] Lin | Four new mitochondrial genomes, and the stability of evolutionary trees of mammals[END_REF][START_REF] Reyes | Congruent mammalian trees from mitochondrial and nuclear genes using bayesian methods[END_REF] and concatenated nuclear genes [START_REF] Murphy | Resolution of the early placental mammal radiation using bayesian phylogenetics[END_REF][START_REF] Delsuc | Molecular phylogeny of living xenarthrans and the impact of character and taxon sampling on the placental tree rooting[END_REF][START_REF] Amrine-Madsen | A new phylogenetic marker, apolipoprotein B, provides compelling evidence for eutherian relationships[END_REF] have converged towards congruent solutions in striking contrast to morpho-anatomical data. According to these new phylogenies, four major groups of placental mammals have been recognised: 1) Afrotheria also known as the African clade (Elephants, Sirenians, Hyraxes, Aardvark, Elephant Shrews, Golden Moles and Tenrecs), 2) Xenarthra (Armadillos, Anteaters and Sloths) a clade of South American endemics, and two distinct groups that comprise most of today's mammalian diversity with 3) Euarchontoglires (Lagomorphs, Rodents, Tree Shrews, Flying Lemurs and Primates) and 4) Laurasiatheria (Insectivores, Bats, Pangolins, Carnivorans, Perissodactyls, Artiodactyls and Cetaceans).

As stated above, phylogenetic reconstruction of deep placental mammal relationships based on mitochondrial genomes had first been hampered by the peculiar evolutionary properties of this molecule. Indeed, the mitochondrial genome evolving about four times more rapidly than the nuclear genome, the recovery of the earliest divergences (37) proved to be difficult due to substitutional saturation. In particular, the first analyses of complete mitochondrial genomes created a controversy surrounding the origin of rodents that were initially found to be paraphyletic [START_REF] D'erchia | The guinea-pig is not a rodent[END_REF]. Indeed, murid rodents (mice and rats) emerged first among placentals in most mitogenomic trees [START_REF] Reyes | Complete mitochondrial DNA sequence of the fat dormouse, Glis glis: further evidence of rodent paraphyly[END_REF][START_REF] Reyes | Long-branch attraction phenomenon and the impact of among-site rate variation on rodent phylogeny[END_REF]. This unexpected finding was actually the result of a long-branch attraction artifact due to the fast evolving murids [START_REF] Philippe | Rodent monophyly: pitfalls of molecular phylogenies[END_REF][START_REF] Sullivan | Are guinea pigs rodents? the importance of adequate models in molecular phylogenetics[END_REF][START_REF] Lin | Four new mitochondrial genomes, and the stability of evolutionary trees of mammals[END_REF]. Thus, for illustrative purposes, we assembled a complete mitochondrial genome data set from 38 species representing all placental orders (Table 4). The concatenation of the 12 H-stranded mitochondrial protein-coding genes (NADH Dehydrogenase subunits 1, 2, 3, 4, 4L and 5; Cytochrome c Oxydase subunits I, II and III, ATP Synthase F0 subunits 6 and 8; and Cytochrome b) led to an alignment of 3,507 amino acid sites after removing ambiguously aligned positions.

NNI heuristic search.

A ML phylogenetic model was first estimated using simultaneous NNI moves to improve a BioNJ starting tree (default option) under the mtMAM substitution model. The shape parameter (α) of a gamma distribution with four bins as well as the proportion of invariants (p-inv) were both estimated from the data. All together, the substitution model is noted mtMAM+Γ4+I. The estimated phylogeny is presented in Fig. 4. This analysis, including the calculation of aLRT-based branch supports, took about 25 minutes to complete on an AMD Opteron 250 2.4 Ghz processor running Linux. The estimated model parameters (α = 0.69; p-inv = 0.34) indicates strong among site rate heterogeneity within the concatenation. The log likelihood of the whole phylogenetic model is -74040.38. Interestingly, three out of the four recognised major placental clades are recovered as monophyletic: Afrotheria, Xenarthra and Laurasiatheria. The fourth group, Euarchontoglires, appears paraphyletic at the base of the tree with murid rodents (Mouse and Rat) emerging first and successively followed by Guinea Pig, Squirrel, Lagomorphs (Hare and Pika) and Primates (Lemur and Human). In fact, this topology differs from what is expected to be the most accurate placental phylogeny because of the position of the murid rodent clade. This peculiar topology is reminiscent of early mitogenomic studies.

It is likely to be the result of a long-branch attraction artifact that is caused by the high evolutionary rate of murid mitochondrial genomes relative to other rodents. Worth noting, however, is the fact that the BioNJ starting topology also presents this apparent rooting artifact (not shown). This observation suggests a potential influence of the starting tree on the outcome of the NNI heuristic search. To check whether the likely artefactual topology found in Fig. 4 is a consequence of using the BioNJ topology as a starting tree for the NNI heuristic search, we conducted the same analysis using a maximum parsimony (MP) tree as the starting topology. Using this strategy results in a slightly better model with respect to the likelihood (-74040.25). The corresponding ML topology only differs from the one in Fig. 4 in the position of Bats (Fruit-eating Bat and Flying-Fox) that now emerge second within Laurasiatheria (not shown). However, as the MP topology also suffers from the rooting artifact, it is still possible that the NNI heuristic search from this starting tree may have reached a local maximum of the likelihood function.

SPR heuristic search.

The maximum likelihood phylogeny reconstructed by applying full-SPR moves to the BioNJ starting tree under the mtMAM+Γ4+I model is presented in Fig. 5. Building the tree took about 1 hour 25 minutes to complete on the same computer, which is almost exactly one hour more than the search using NNI moves. As previously, aLRT branch supports were also calculated and the time needed to work out these values was approximately 5% of the total computing time. The estimated model parameters (α = 0.68; p-inv = 0.34) are almost identical to the ones obtained using the NNI search strategy. This result confirms that ML estimation of model parameters is relatively insensitive to topological differences induced by the different heuristics. However, SPR moves converged to a model that has a greater log likelihood (-74021.74) than the NNI-based one (-74040.38).

The SPR topology presents major differences with the NNI one since the respective monophyly of each of the four major placental clades is now recovered (Fig. 5). Indeed, the ML topology is no longer rooted on the Mouse/Rat ancestral branch with Euarchontoglires being monophyletic and Afrotheria now appearing as the earliest placental offshoot. This topology is almost fully compatible with the new placental phylogeny inferred from large concatenated data sets of mainly nuclear genes [START_REF] Murphy | Resolution of the early placental mammal radiation using bayesian phylogenetics[END_REF][START_REF] Delsuc | Molecular phylogeny of living xenarthrans and the impact of character and taxon sampling on the placental tree rooting[END_REF][START_REF] Amrine-Madsen | A new phylogenetic marker, apolipoprotein B, provides compelling evidence for eutherian relationships[END_REF]. This result illustrates that, in this particular case, NNI-based heuristic is trapped in a local maximum because of insufficient tree space exploration conditioned by the BioNJ or MP starting topologies. Note however that the NNI-based and SPR-based ML phylogenetic models are not statistically different according to a SH test [START_REF] Shimodaira | Multiple comparisons of log-likelihoods with applications to phylogenetic inference[END_REF]. Hence, while simultaneous NNIs get trapped in a local maximum, the estimated ML solution is not different from the one found by SPRs from a statistical point of view. In other words, NNI tree searching is probably trapped in this local optimum because of the lack of phylogenetic signal for certain parts of the tree.

3.4 SPR heuristic search starting from random trees.

SPR searches can also be initiated with random starting trees (see section 2.4). The full-SPR-based tree search procedure was thus repeated 10 times, each analysis starting from a different random phylogeny. In eight of these 10 replicates, the ML topology of Fig. 5 compatible with the new placental phylogeny was recovered (log-likelihood: -74021.74).

However, in the two remaining replicates (20%), the SPR heuristic search converged to the alternative suboptimal topology (log-likelihood: -74040.25) that was previously found when using a MP starting tree combined with the NNI heuristic search (see section 3.2).

These results suggest that the likelihood surface for this placental mitogenomic data set is dominated by at least two peaks relatively close in likelihood but distant in the space of tree topologies, one probably corresponding to the ML topology and the other to the suboptimal topology caused by the rooting artifact.

Assessing statistical support for internal edges.

Statistical supports for branches displayed by each of the two competing topologies were measured using 100 non-parametric bootstrap replicates (BP; [START_REF] Felsenstein | Confidence limits on phylogenies: an approach using the bootstrap[END_REF]) and the approximate likelihood ratio test for branches (aLRT; (8)). BP and aLRT values (more precisely, pvalues obtained from SH tests) are reported on Fig. 4 and 5 for nodes that show at least 50% bootstrap support. For both topologies, the statistical support is almost maximal for the respective monophyly of Afrotheria, Xenarthra and Laurasiatheria. The monophyly of Euarchontoglires recovered in the SPR topology (Fig. 5) is not statistically supported, as is also its paraphyly induced by the rooting artifact in the NNI topology (Fig. 4).

In fact, the differences between the two topologies only involve nodes with weak support from the data.

Bootstrap and aLRT values largely agree on most branches of both trees. Indeed, branches with bootstrap supports close to 100 also have aLRT values close to 1.0 in most cases. However, a few branches are well supported according to aLRT values but have small bootstrap proportions (see the branch at the root of the Xenarthra, Afrotheria and Laurasiatheria clades, with BP=58 and aLRT=1.0). Such differences between aLRT and bootstrap values are probably the consequence of a lack of phylogenetic signal to resolve specific parts of the phylogeny. More work still needs to be done to fully understand such result though.

4 Notes.

1. The hypothesis of site independence is relaxed in certain models. For instance, codon-based models impose a constraint on groups of columns that belong to the same codon site. Felsenstein and Churchill (44) also proposed a model where the rates of substitutions at adjacent sites are correlated. Also, several models describe the evolution of pairs of interacting nucleotides among ribosomal RNA molecules [START_REF] Schniger | A stochastic model for the evolution of autocorrelated DNA sequences[END_REF][START_REF] Muse | Evolutionary analyses of DNA sequences subject to constraints on secondary structure[END_REF][START_REF] Tillier | High apparent rate of simultaneous compensatory basepair substitutions in ribosomal rna[END_REF].

2. It is not very clear which of non-deterministic or deterministic approach is the best for phylogenetic inference. Note however that deterministic methods are generally faster than stochastic optimisation approaches and sufficient for numerous optimisation problems [START_REF] Aarts | Local search in combinatorial optimization[END_REF]. On the other hand, stochastic methods have the ability to find several near-optimal trees, which gives an idea of the inferred tree variability.

3. Under the default settings, PhyML modifies the tree topology using NNI moves and simultaneously optimises branch lengths. However, when the tree topology estimate is stable (i.e., no improvement of the likelihood can be found by modifying the current tree topology), the optimisation concentrates on branch lengths and parameters of the Markov model in order to save computing time. NNI moves with optimisation of the central and the four adjacent branch lengths are also systematically tested during the very last optimisation step. This last step frequently finds a modification of the tree topology that was not detected by the other approximate (but fast) tree topology search methods (i.e., NNI, NNI+SPR or full-SPR).

4. The SPR-search strategy implemented in PhyML 3.0 actually relies on filtering SPR moves using the parsimony criterion instead of a distance-based approach. Indeed, parsimony and likelihood are closely related from a statistical perspective and the analysis of real and simulated data (Dufayard, Guindon, Gascuel, unpublished) have demonstrated the benefits of using a parsimony-based filter.

5. PhyML binary file must be located in a directory listed in the PATH environment variable if the program is launched from a command-line window. The program can also be launched by typing './phyml' provided that PhyML binary file is in the current directory. Launching PhyML by clicking on the corresponding icon is not recommended. In case PhyML can not find the sequence data file when launched by clicking on the icon, we suggest using a command-line window.

6. Questions regarding the amount of memory required can be eluded using the '-DBATCH' flag when compiling the program. This option is available through a simple modification of the Makefile. It is highly recommended to use this option when launching PhyML in batch mode or when comparing run times of different programs.

7. The estimation of bases or amino-acid frequencies relies on a iterative method that takes into account gaps and ambiguous characters. The frequencies of the nonambiguous characters (bases or amino-acids) at step n are functions of the counts of the non-ambiguous characters plus the counts of the ambiguous characters weighted by the probabilities of the non-ambiguous characters estimated at step n -1. These probabilities correspond to the frequencies estimated at step n -1. The same approach is also used in PAML [START_REF] Yang | PAML : a program package for phylogenetic analysis by maximum likelihood[END_REF] and PHYLIP (12) programs.

8. PhyML uses an original method that simultaneously estimates the gamma shape parameter and the proportion of invariants (Guindon and Gascuel, unpublished). This method relies on the observation that the two parameters show a strong quasi-linear positive relationship. Basically, the gamma shape parameter is first estimated using a standard one-dimensional optimisation method. The proportion of invariants is then deduced from the linear relationship with the gamma shape parameter. Hence, the estimation of the proportion of invariants does not rely on time-consuming optimisation methods. For each sub-menu, the defaults options can be altered by entering the relevant keys. The data analysis is launched once the 'y' (or 'Y') key is entered.

F84 ( 12

 12 ), HKY85 (13), TN93 (14) GTR[START_REF] Lanave | A new method for calculating evolutionary substitution rates[END_REF][START_REF] Tavaré | Some probabilistic and statistical problems on the analysis of DNA sequences[END_REF] and CUSTOM for nucleotides ; WAG[START_REF] Whelan | A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach[END_REF],Dayhoff (18), JTT[START_REF] Jones | The rapid generation of mutation data matrices from protein sequences[END_REF], Blosum62[START_REF] Henikoff | Amino acid substitution matrices from protein blocks[END_REF], mtREV[START_REF] Adachi | MOLPHY version 2.3. programs for molecular phylogenetics based on maximum likelihood[END_REF], rtREV[START_REF] Dimmic | rtREV: an amino acid substitution matrix for inference of retrovirus and reverse transcriptase phylogeny[END_REF], cpREV[START_REF] Adachi | Plastid genome phylogeny and a model of amino acid substitution for proteins encoded by chloroplast DNA[END_REF], DCMut[START_REF] Kosiol | Different versions of the Dayhoff rate matrix[END_REF], VT (25) and mtMAM[START_REF] Cao | Conflict among individual mitochondrial proteins in resolving the phylogeny of eutherian orders[END_REF] for amino acids. Nucleotide equilibrium frequencies are estimated by counting the occurrence of A, C, G and T s in the data (seeNote 7).These frequencies can also be adjusted in order to maximise the likelihood of the phylogenetic model (Substitution Model sub-menu, option 'F' ; command line argument : '-f e') or deduced from the actual sequences. Amino acid equilibrium frequencies are either deduced from the actual sequences (Substitution Model sub-menu, option 'F' ; command line argument : '-f e') or given by the substitution models themselves (default option).The CUSTOM option provides the most flexible way to specify the nucleotide substitution model. The model is defined by a string made of six digits. The default string is '000000', which means that the six relative rates of nucleotide changes:A ↔ C, A ↔ G, A ↔ T , C ↔ G, C ↔ Tand G ↔ T , are equal. The string '010010' indicates that the rates A ↔ G and C ↔ T are equal and distinct from A

Figure 1 .

 1 Figure 1. Comparison of bootstrap and SH-like branch supports, using C8 alpha chain precursor data set from TREEBASE. This data set displays a strong phylogenetic signal overall as most internal branches (13/19) have support close to their maximum values with both approaches. SH-like supports are non-significant for 5 branches (noted with '△') that also have low bootstrap proportions, while one branch with a relatively low bootstrap proportion (72 in italic) is strongly supported according to the aLRT test (SH-like value : 0.97).

Figure 2 .

 2 Figure 2. Text-based interface to PhyML. The PHYLIP-like text interface to the program is organised in sub-menus. The '+' and '-' key are used to cycle through them.

Figure 3 .

 3 Figure 3. Typical DNA sequence alignments (a) and input trees (b). Sequence names do not contain any blank character and at least one blank separates each name from the corresponding sequence. Trees are in standard NEWICK format.

Figure 4 .

 4 Figure 4. ML phylogeny reconstructed using simultaneous NNI moves. The log-likelihood of the corresponding phylogenetic model is -74040.38. Numbers in the tree correspond to non-parametric bootstrap supports (100 replicates) and p-values of the approximate likelihood ratios (SH-test). Values are reported only for nodes with BP > 50 and stars indicate nodes that received 100% support. S.-e. Elephant Shrew is for Shorteared Elephant Shrew.

Figure 5 .Table 1 .LTable 2 .

 512 Figure 5. ML phylogeny reconstructed using SPR moves. The log-likelihood of the corresponding phylogenetic model is -74021.74. Numbers in the tree correspond to non-parametric bootstrap supports (100 replicates) and p-values of the approximate likelihood ratios (SH-test). Values are reported only for nodes with BP > 50 and stars indicate nodes that received 100% support. S.-e. Elephant Shrew is for Short-eared Elephant Shrew.

Figure 2

 2 Figure 2.

Figure 3 .Figure 5 .

 35 Figure 3.

Acknowledgements

This work was supported by the 'MITOSYS' grant from ANR. The chapter itself is the contribution 2007-XXX of the Institut des Sciences de l'Evolution (UMR5554-CNRS)