Dr Elise Buisson
email: elise.buisson@univ-avignon.fr

L Achour

K Guediri

S Bissati

H Abdelkrim

A Eddoud

Changes in weed species composition in irrigated

Keywords: agroecosystems, camel rangeland, cereal center-pivot irrigation plots, date palm plantations, greenhouses, market gardening

published or not. The documents may come

Introduction

Global changes influence weed distribution, species dominance, and other aspects of weed communities in agroecosystems (e.g. climate change: [START_REF] Zicka | Weed Biology and Climate Change[END_REF]. Changes in agricultural practices associated with intensified crop production are by far the leading root cause of weed dynamics: this follows directly from the ecological implications of increased seeding and planting, increased fertilizer and pesticide use, and simplified crop rotation [START_REF] Kleijn | On the relationship between farmland biodiversity and land-use intensity in Europe[END_REF]. On the one hand, the difficulty in controlling weeds is one of the major reasons for which farmers are forced to relinquish cultivated plots and to plow new ones; on the other hand, some weed control practices have proven detrimental to many of the plant and animal communities that are associated with agricultural environments [START_REF] Flohre A | Does soil biota benefit from organic farming in complex vs. simple landscapes? Agriculture, Ecosystems[END_REF] while weeds can play an important role in supporting diversity for higher trophic levels [START_REF] Marshall Ejp | The role of weeds in supporting biological diversity within crop fields[END_REF]. Whereas these changes occur on weed communities, a thorough knowledge of the composition of weeds, along with their spatial and temporal distribution, is a prerequisite for good weed management in agroecosystems.

In Algeria, several studies have increased our understanding of the composition and distribution of weed flora in various northern regions of the country, with a Mediterranean climate, such as: on the Mediterranean central coast [START_REF] Chevassut | Contribution à l'étude des groupements de mauvaises herbes de la région d'el Harrach[END_REF][START_REF] Boulfekhar N | Etude phytosociologique de certains groupements de mauvaises herbes dans la plaine de la Mitidja, Algérie septentrionale[END_REF][START_REF] Fenni M | Contribution à l'étude des groupements messicoles des Hautes Plaines sétifiennes[END_REF][START_REF] Abdelkrim H | Contribution à la connaissance des groupements de mauvaises herbes des cultures du secteur algérois : Approche syntaxonomique et phénologique[END_REF], and on the Mediterranean north-eastern Algeria [START_REF] Fenni M | Contribution à l'étude des groupements messicoles des Hautes Plaines sétifiennes[END_REF] and north-western Algeria [START_REF] Kazi Tani Ch | Aspects floristiques de la flore des champs du domaine phytogéographique oranais (Nord-Ouest algérien) et persistance d'espèces rares et endémiques[END_REF] coast. However, many agricultural regions of the country remain under-characterized. Indeed, there currently exist no studies on the weed flora of the Saharan Algeria.

Algeria has a flora comprising about 3139 species, of which 653 species (accounting for 12.6% of the total number of species) are endemic and distributed among nearly 150 families.

Among the endemic flora, more than 25% of the species are associated with crops (Quézel & Santa, 1962-1963). Botanical and plant ecology studies in the Sahara regions are scarce and mainly oriented towards the native flora of natural areas, such as that of camel rangelands [START_REF] Ozenda P | Flore du Sahara, 3rd edition[END_REF][START_REF] Chehma | Productivité fourragère et capacité de charge des parcours camelins du Sahara septentrional Algérien[END_REF] and conservation and narrow endemic species [START_REF] Brito | Unravelling biodiversity, evolution and threats to conservation in the Sahara-Sahel[END_REF][START_REF] Migliore | Surviving in mountain climate refugia: new insights from the genetic structure and diversity of the relict shrub Myrtus nivellei (Myrtaceae) in the Sahara desert[END_REF].

Our present work aims to show whether the composition of the weed flora associated with irrigated crops in southeastern Saharan Algeria has changed over a period of twenty years from 1990 to 2010. We studied a set of cereal fields, greenhouses and palm plantations.

We hypothesized that the weed flora of cereal fields and palm plantations would change little over the study period because cereal cultivation practices are unchanged since the 1980's and interventions are relatively minor in palm plantations after their establishment. We expected more changes in greenhouses where practices potentially evolve more rapidly.

Materials and Methods

Study area

The study was conducted in Hassi Benabdellah, located in northeastern part of Saharan Algeria, about 800 km from Algiers (32°01'N; 5°28'E; Fig. 1). The northern Sahara, which covers an area of 1 million km², receives precipitation in the winter and yields 100±50 mm of rainfall per year [START_REF] Valles | Soils and brine geochemistry and mineralogy of hyperarid desert playa, Ouargla Basin, Algerian Sahara[END_REF]. It is essentially a transition zone between the northern African Mediterranean steppes and the central Sahara desert, and it was selected mainly for its agricultural potential. Indeed, many important crops, such as dates, greenhouse gardens, and cereals grown using center-pivot irrigation (sometimes called central pivot irrigation or circle irrigation), are all farmed in this region. A relatively recent import to the area, intensive crop cultivation, first appeared in the region in the late 1960s and early 1970s, mainly occupying what had once been open camel rangeland. Although irrigation is applied all year long in the date palm plantations (surface irrigation; once a week from November to March and 2 to 3 times a week the rest of the year), it is only applied during the growing season on cereal plots (center pivot irrigation; 866 mm/m²) and market gardening greenhouses (surface irrigation; twice a week from November to February and 3 times a week from March to May). None of the plots selected for this study were treated with herbicides.

Sampling method

The present study covers the twenty-year period beginning in 1990, and handles the full twenty years in two stages. The first stage relies on a synthesis of all the data from previous studies mentioning exhaustive weed flora sampling in the ITCMI (Institut Technique des Cultures Maraichères et Industrielles) and ITDAS (Institut Technique de Développement de l'Agriculture Saharienne) research stations between 1990 and 2002. For that period, one presence-absence plant list was compiled per year surveyed per crop type (Table 1; [START_REF] Itcmi | Rapport d'activité de recherche de la station expérimentale de Hassi Ben AbdAllah-Ouargla pour la période allant de 1990 à 2002[END_REF]. The second stage of the analysis, covering the period from October 2003 to June 2010, is based on weed surveys conducted by the authors according to the following protocol. For both stages, to match the seasonal planting and harvest schedule, surveys were carried out starting in October of each year (x), and repeated in February and June of the following year (x + 1), to give a total of three surveys per year. The data from each set of three surveys were then taken together to fully characterize each October x to June x+1 crop year. We used the Base de données des plantes d'Afrique for plant nomenclature [START_REF] Cjb | Base de données des plantes d'Afrique -CJB Conservatoire et Jardin Botaniques de la[END_REF]. The botanical attributions for weed species are listed in Appendix 1 in Supplementary Material.

Sampling in cereal plots

Three 30-ha cereal center-pivot irrigation plots, hereafter referred to as cereal plots, cultivated with either hard (Triticum turgidum L. subsp. durum (Desf.) Husn.) or soft (Triticum aestivum L. subsp. aestivum) wheat, were selected. These plots started to be cultivated in the mid 1980's. In each cereal plot, three 100m-long and 2m-wide transects were randomly placed and sampled by establishing a presence-absence plant list three times per year between wheat seeding and harvesting. The three samples taken on the same transect year were pooled for analysis (n = 63; 6 sites × 3 transects × 7 years). The pre-2003 data, collected over 8 years, were considered together with the post-2003 data in the multivariate analyses (with n = 63 + 8 = 71).

Sampling in market gardening greenhouses

Six 400-m² greenhouses, cultivated since the early 1980's, were selected. In each greenhouse, three 50m-long and 1m-wide transects were randomly placed and sampled by establishing a presence-absence plant list three times per year, although the exact time of year at which sampling occurred varied according to crop type (and therefore crop cycle). The three samples taken within the transect year were pooled for analysis (n = 126; 6 sites × 3 transects × 7 years). The pre-2003 data, which came from for 13 years (Table 1), were considered together with the post-2003 data in the multivariate analyses (with n = 126 +13 = 139).

Sampling in date palm plantations

For the palm cultivation portion of the study, three one-hectare date palm (Phoenix dactylifera L.) plantations without understory crops were selected. These plantations had been established in 1972. In each plantation, three 100m-long and 10m-wide parallel transects were randomly placed and sampled by establishing a presence-absence plant list three times per year. The three samples taken on the same transect from October in the year x to June in the year x+1 were pooled for analysis. This means that a given crop year was statistically characterized by a total of 3 sites × 3 transects × 7 years, n = 63). The pre-2003 samples consist of data collected over 9 years (Table 1), and are considered in the multivariate analyses (with n=72, following from 9 + 63 = 72) in which the pre-2003 and post-2003 datasets are handled together.

Plant species study

In order to analyze the characteristics of the species found on the various plots, all species were classified according to (1) their family; (2) their chorotype and (3) their life-form according to Raunkiaer's life form using the following floras: Quézel & Santa (1962-1963), [START_REF] Ozenda P | Flore du Sahara, 3rd edition[END_REF], [START_REF] Tison | Flora Gallica : Flore de France[END_REF] and Danin (2015) (see Appendix 1 in Supplementary Material). The chorotypes were simplified in 7 categories: 1. Alien species from Americas, Australia and South Africa, 2. A complex chorotype containing species with biregional, triregional, and multiregional origins, 3. Euro-Siberian (ES), 4. Irano-Turanian stricto sensu (IT), 5. Mediterranean stricto sensu (M), 6. Saharo-Arabian stricto sensu (SA), and 7. Sudanian (S).

Data analysis

We analyzed differences in percentages of weed species / transect for each crop type separately, depending on the explanatory variable "chorotypes" (or "life forms)' using generalized linear models with binomial error distribution. When the explanatory variable was significant, we conducted pairwise comparisons with Tukey's pairwise comparison tests using function 'glht' in the package "multcomp" in R (R Core Team 2018). For all tests, Pvalues lower than 0.05 were considered as statistically significant.

We analyzed differences in number of weed species / transect along time for each crop types separately. Two Mixed-effects Generalized Linear Models were compared for each crop type: one with "years" as an explanatory variable and one without. For these analyses, only the data from 2003 to 2010 were used. If the first model had a lower AIC (Akaike Information Criterion), time was considered as having a significant effect and we conducted pairwise comparisons with Tukey's pairwise comparison tests using function 'glht' in the package "multcomp" in R (R Core Team 2018). The same procedure was carried out to assess the differences in percentages of weed species / transect along time for particular chorotypes or life forms, considering crop types separately.

The relationships between explanatory variables (crop types or years) and weed assemblages were explored using Correspondence Analyses (CA). To describe the vegetation corresponding to each crop type, a first CA was run on all vegetation surveys (282 samples × CAs were run 1) on the palm plantation samples (72 points × 123 species), 2) on the cereal plot samples (71 points × 110 species) and 3) on the greenhouse samples (139 points × 32 species). These analyses were followed by a Multiple Response Permutation Procedure analysis (MRPP) to determine the statistical significance between the crop types or the years.

On the figures, for clarity, only the species the most correlated to axes 1 and 2 are shown. All of the analyses were conducted with R 3.4.4 (R Core Team 2018), univariate analyses with its package "lme4" and "multcomp" and multivariate analyses with its packages "ade4" and "vegan".

Results

Weed composition depending on crop types

The number of weed species varied from 119 species in the palm plantations (9000 m² explored every year) to 104 in cereal plots (1800 m² explored every year) and 31 species in greenhouses (900 m² explored every year). In total, during the two sampling periods (1990-2002 and 2003-2010), we identified 162 weed species from 40 botanical families, three of which made up for almost half of the species in the study (45.1% of the species): Poaceae (20.4%), Asteraceae (14.2%) and Brassicaceae (10.5%) of the species. Poaceae (19.2%, 35.5% and 22.7% of the weed species in cereal plots, greenhouses and palm plantation respectively), Asteraceae (18.7%, 16.1%, 14.3%) were the dominant families in all crop types. The third most represented family was Brassicaceae in cereal plots (14.4% of the species). In greenhouses and palm plantations Amaranthaceae was the third most represented family (12.9% and 7.6% respectively) (see Appendix 1 in Supplementary Material).

The MRPP analysis, as run on all the samples, showed a significant difference in weed species composition between the crop types (observed delta: 0.37; expected delta: 0.60; p-value<0.001). The first two first axes of the CA explained 31.0% of the total inertia. As shown in Fig. 4, axis 1 (16.4%) delineated cereal plot samples correlated with Anisantha madritensis, Avena barbata, Daucus carota, Emex spinosa, Fagonia glutinosa, Launaea capitata, Raphanus raphanistrum, Sinapis arvensis, and Trigonella polycerata, from greenhouse samples correlated with Amaranthus retroflexus, Convolvulus arvensis, Cynodon dactylon, Cyperus rotundus, Digitaria nodosa and Senecio vulgaris (Fig. 4; the botanical attributions for weed species are listed in Appendix 1 in Supplementary Material). Axis 2 (14.6%) separated these samples from the palm tree plantation samples correlated with Cressa cretica, Tamarix gallica, Juncus maritimus, Euphorbia chamaecyse, Sonchus maritimus, and Phragmites australis (Fig. 2).

Weed chorotypes and life-forms depending on crop types

The study area, located in the northeastern part of Saharan Algeria, was characterized by a set of species from various origins. Overall, 47% of the species belong to the complex chorotype, 25% are Saharo-Arabian species and 18% are Mediterranean species. The other categories (Alien, Euro-Siberian, Irano-Turanian, Saharo-Arabian and Sudanian) accounted for 10% of the species. The most common weed species on the three crop types were: Chenopodium In all three crop types, the most common chorotype category was the complex chorotype (Fig. 3). The second most common chorotype category in the cereal plots was Mediterranean (e.g. Lolium multiflorum, Avena barbata, Melilotus indicus, Sinapis arvensis), while it was Aliens (e.g. Erigeron canadensis, Cyperus rotundus) in the greenhouses and Saharo-Arabian (e.g. Moltkiopsis ciliate, Ifloga spicata, Launaea capitate, Tetraena alba) in the palm plantations (Fig. 3).

Overall, 'therophytes' was the most represented life-form (Fig. 4). Palm plantations were the crop type with the greatest diversity of life-forms, including perennials and a significant amount of hemicrytophytes (Fig. 4; Appendix 1 in Supplementary Material). Few species were geophytes and perennials, mainly found in the palm plantation. More than 60% of the species exhibiting a Mediterranean range or belonging to our complex chorotype category were therophytes. This life-form represented less than 50% of the species with a Saharo-Arabian range, while other species with this origin mainly belonged to chamaephytes or hemicryptophytes (only one geophyte: Euphorbia guyoniana).

Weed species richness changes with time

The number of weed species in cereal plots tended to decrease from 1990 to 2003 and then significantly increased from 2003 to 2010 (df=6, F=22.01, p-value<0.001). In greenhouses, the number of weed species tended to be constant around 18.1±0.6 weed species / transect from 1990 to 2010, with a few punctual increases (e.g. in 2004-05: 22.4±1.2; df=6, F=5.21, p-value<0.001). The number of weed species in palm plantations tended to increase from 1990 to 2003 and then stabilized around 52.3±3.0 weed species / transect (df=6, F=1.42, p-value=0.224) (Fig. 5).

Weed chorotypes and life-forms changes with time

Significant variations of percentages of therophytes were observed in the greenhouses, with no specific trends (Table 2). Variations of percentages of therophytes in the cereals plots and of therophytes and hemicryptophytes in the palm plantations were not significant.

Significant variations of percentages of species in the complex chorotype category were observed in the cereal plots and in greenhouses, with no specific trends (Table 2; Fig. 6).

Percentages of alien species significantly increased in the greenhouses, while percentages of Saharo-Arabian species decreased in the palm plantations (Table 2; Fig. 6).

Weed composition changes with time

The CA, as run on cereal plot data, showed a clear change in species composition with time and this was confirmed by the MRPP analysis showing a significant difference between the dates (observed delta: 0.17; expected delta: 0.37; p-value<0.001); the first two first axes of the CA were found to explain 37.6% of the total inertia. As shown in Fig. 7 The CA, as run on the greenhouse data, only showed that some years (1990-91, 2004-05, 2006-07) were different from others (Fig. 8) and this was confirmed by the MRPP analysis showing a significant difference between the dates (observed delta: 0.22; expected delta: 0.30; p-value<0.001); the first two first axes of the CA were found to explain 32.2% of the total inertia. Axis 1 (17.1%) delineated the 2004-05 samples, correlated with Dactyloctenium aegyptium and Solanum nigrum, from the samples from other years. Axis 2 (15.1%) distinguished the oldest samples (1990-1991) correlated with Cutandia dichotoma and Spergularia marina (Fig. 8).

The CA, as run on the palm tree plantation data, also showed a clear change in species composition with time and this was confirmed by the MRPP analysis showing a significant difference between the dates (observed delta: 0.37; expected delta: 0.49; p-value<0.001); the first two first axes of the CA (axis 1: 16.6% and axis 2: 11.3%) explained 27.9% of the total inertia. As shown in Fig. 9, Agathophora alopecuroides and Stipagrostis plumosa were mainly found in the early years, followed by Anethum graveolens and Echium pycnanthum at the end of the 1990's and by Glaucium corniculatum and Ogastemma pusillum the start of the 21 st century. Around the late 2000's, Cressa cretica, Cyperus rotundus, Phragmites australis Tamarix sp. became more common (Fig. 9).

Discussion

Overall weed composition

The three most common families encountered in this study, namely Poaceae, Brassicaceae and Asteraceae were such a dominating presence because they are important families of the Algerian flora, especially in the Sahara [START_REF] Ozenda P | Flore du Sahara, 3rd edition[END_REF], and they are adapted to disturbed environments [START_REF] Holm Lg | The World's Worst Weeds: Distribution and Biology[END_REF][START_REF]Mauvaises herbes des céréales de la Chaouia (Maroc) II. Aspect écologique[END_REF]. Indeed, Holm et al. (1977) report that Poaceae, Cyperaceae, Asteraceae, Polygonoceae, Amarantaceae, Brassicaceae and Fabaceae constitute the major crop weed families. Other studies in semi-arid Mediterranean and tropical environments have also shown that Amarantaceae, Asteraceae, Brassicaceae, Fabaceae, Poaceae, Polygonoceae, and Solanaceae are very common in the vegetation and seed bank [START_REF] Dorado | Weed seedbank response to crop rotation and tillage in semiarid agroecosystems[END_REF][START_REF] Salama Fm | Weed flora of common crops in desert reclaimed arable lands of southern Egypt[END_REF]. Our results are also comparable to those reported by studies in Mediterranean Algeria [START_REF] Boulfekhar N | Etude phytosociologique de certains groupements de mauvaises herbes dans la plaine de la Mitidja, Algérie septentrionale[END_REF], Morocco [START_REF] Bouhache | Aspects floristico-agronomiques des mauvaises herbes de la région du Loukkos (Maroc)[END_REF] and Egypt (Salama et al., 2016). The most common weed species recorded in this study (Chenopodium murale, Lolium multiflorum, Melilotus indicus, Lysimachia arvensis,

Sonchus oleraceus, Setaria verticillata, Polypogon monspeliensis and Erigeron canadensis)

are recognized as major weeds elsewhere in the world (e.g. USDA, NRCS 2015), with their status as invasive depending on the region considered. Chenopodium murale, Melilotus indicus, Sonchus oleraceus and Polypogon monspeliensis were also found to be dominant weed species in wheat, millet or alfa-alfa fields studied by [START_REF] Salama Fm | Weed flora of common crops in desert reclaimed arable lands of southern Egypt[END_REF].

The proportion of weeds exhibiting Mediterranean or Saharan ranges was high, accounting for 43% of the species, but 47% of the weeds belong to our complex chorotype category, with biregional, triregional, and multiregional origins. According to [START_REF] Ozenda P | Flore du Sahara, 3rd edition[END_REF], the flora of northern Sahara is characterized by a preponderance of the Saharan species with only some incursions of Mediterranean species. This apparent disparity between our observations and Ozenda's prediction can probably be explained by habitat differences.

Ozenda described the flora of all Saharan habitats, including e.g. camel rangelands which experience recurrent disturbances due mainly to drought and grazing and which are thus characterized by an adapted Saharan flora. On the other hand our study was conducted in agroecosystems which are highly disturbed environments. Cultivated fields represent a particular kind of artificial ecosystem that is subject to tillage, pesticides, fertilizers, irrigation and other interventions, thus experiencing an exceptional frequency and intensity of disturbances intended to maintain favorable conditions for production [START_REF] Smith Cw | Crop Production: Evolution, History, and Technology[END_REF].

Artificial watering, in particular, provide conditions favorable to Mediterranean species [START_REF] Steinmaus | Growth analysis and canopy architecture of velvetleaf grown under light conditions representative of irrigated Mediterranean-type agroecosystems[END_REF]. These disturbances lead to a banalization of the weed flora [START_REF] Kazi Tani Ch | Aspects floristiques de la flore des champs du domaine phytogéographique oranais (Nord-Ouest algérien) et persistance d'espèces rares et endémiques[END_REF] and a depletion of endemic species in this type of environment [START_REF] Jauzein P | Biodiversité des champs cultivés : l'enrichissement floristique[END_REF]. With the Mediterranean and complex chorotype species being mainly therophytes (over 60%), these provenances create a dominant therophytic component in the Saharan weed assemblage.

Differences in weed composition between crop types

The weed flora is mainly related to soil and climatic conditions, but also varies in composition and number of weed species depending on crop types and the agricultural practices associated with these crops [START_REF] Zimdahl Rl | Fundamentals of Weed Science[END_REF][START_REF] Salama Fm | Weed flora of common crops in desert reclaimed arable lands of southern Egypt[END_REF]. These differences are explained by the differences in crop management sequence, particularly where it concerns tillage. In greenhouses, tillage occurs continuously throughout the growing season, as part of incorporating mineral fertilizers, weeding, and hoeing (to aerate the soil). In addition, greenhouse crops (Solanaceae, Cucurbitaceae, Fabaceae) offer more possibilities for rotation compared to the two other crop types, which are limited to a single botanical family (Poaceae in cereal plots and Arecaceae in palm plantations). Greenhouse crops therefore have a greater diversity of agricultural practices, and this can explain their relatively low number of weed species [START_REF] Pyšek | Effects of abiotic factors on species richness and cover in Central European weed communities[END_REF][START_REF] Armengot | Landscape and land-use effects on weed flora in Mediterranean cereal fields[END_REF]. It can also explain the presence of Amaranthus retroflexus, Chenopodium murale and Senecio vulgaris, which have a very rapid life cycle. Ball (1992) also reports that in irrigated market gardening crops, tillage and crop rotation are the main factors explaining the floristic composition of weeds and the seed bank, which in turn explain the number of weed species of the vegetation in the plot.

In the cereal plots, tillage only occurs once a year at the time of crop establishment (tillage + sowing) in the fall (October). This provides time for Saharan species from the surroundings camel rangelands to colonize the cereal plots: e.g. Atractylis delicatula, Carthamus eriocephalus, Ifloga spicata, and several species of Launaea. The cereal plots are characterized by a second group of weed species, of Mediterranean origin, probably sowed with wheat seeds [START_REF]New weed species recently introduced into Morocco[END_REF]): e.g. Avena barbata, Daucus carota, Emex spinosa, Raphanus raphanistrum, Sinapis arvensis, and Trigonella polycerata.

In palm plantations without understory crops, tillage is limited to a single hydroagricultural development at the time of tree planting in which ditches are created to ensure irrigation of the date palms [START_REF] Zaid A | Date palm cultivation[END_REF]. Once the plantation is established, the soil is not tilled anymore. Palm plantations are thus characterized by perennial species: chamaephytes (Cynodon dactylon), geophytes (Cyperus rotundus, Cressa cretica, Colchicum gramineum, Oxalis pes-caprae), hemicryptophytes (Phragmites australis, Erodium glaucophyllum, Juncus maritimus, Sonchus maritimus), and phanerophytes (Tamarix gallica). They are also characterized by weed species that prefer moist soil and can tolerate salt, such as Cressa cretica, Euphorbia chamaecyse, and Sonchus maritimus. Salt is a result of the high evaporation rate characterizing arid and semiarid climates, an effect that causes salt to accumulate in the top layer of the soil (Al-Amoud, 2010). It can also come from irrigation with marginal groundwater resources and from a rise in the water table [START_REF] Bouarfa | Salinity patterns in irrigation systems, a threat to be demystified, a constraint to be managed: Field evidence from Algeria and Tunisia[END_REF].

Changes in weed composition over time

In cereal plots, old samples are characterized by species found on the camel rangelands of the region, such as Ammosperma cinereum, Atractylis delicatula, Brocchia cinerea, Diplotaxis acris, Launaea mucronata, Malcolmia africana, Malva aegyptia and Paronychia arabica.

With time, a mix of weeds (e.g. Sonchus asper) and segetal species (e.g. Papaver rhoeas, Galium tricornutum, and Scandix pecten-veneris) colonized these plots. These species are probably sown with the wheat seeds coming from the northern part of the country as has been shown to occur in Morocco [START_REF]New weed species recently introduced into Morocco[END_REF]. While some species native to the camel rangelands progressively disappear from the cereal plots (Atractylis delicatula, Carduus getulus, Cornulaca monacantha, Paronychia arabica, and Sisymbrium reboudianum) others manage to maintain themselves (e.g. Henophyton deserti and Tetraena alba), although species of both groups belong to various life-forms (chamaephytes, hemicryptophytes and therophytes). On the other hand, a study run from 1985 to 2005 in southern Spain showed that weed diversity did not change over time [START_REF] Hernandez Plaza | Tillage system did not affect weed diversity in a 23-year experiment in Mediterranean dryland[END_REF]. This might be explained by the fact that the fields had been cultivated with wheat for 10 years prior to the start of the study. This is unlike our own study, initiated just a few years after the start of cultivation in the mid 1980's. In order to better understand the drivers of change on the cereal plots, a thorough survey among farmers would need to investigate wheat seed provenances and changes in practices that might have gone unnoticed. A study of the species sown along with wheat seeds would also be useful, but might not to be feasible as one would need samples from seeds sown from 1990 to 2010. Finally, a seed bank study in cereal plots of various ages would provide information on native rangeland species depletion.

In the early years, indigenous species were already scarce in greenhouses which is probably due to the fact that these plots had been cultivated for 10 years before the start of the study. This crop type also induce the most disturbances of the three studied crop types. This and the frequent crop rotations (Solanaceae, Cucurbitaceae, Fabaceae) most probably explain the lack of patterns of weed communities through time. Indeed, the weed flora has not changed much in greenhouses over the last twenty years and is characterized by anthropophilous or widespread species (e.g., Amarantus retroflexus, Chenopodium murale, Erigeron canadensis, and Sonchus oleraceus) probably sown with crop seeds. Indeed, crop seeds are often bought from the northern part of the country and are not pure. Studying the species sown along with crop seeds would be useful to better understand weed species composition. A fertilization survey among farmers would also be helpful as they may have various kinds of practices, such as the use of chemical fertilizers, local organic fertilizers (ovine, caprine, mixed manure) or organic fertilizers from other parts of the country (ovine, caprine, bovine, mixed manure).

The vegetation of palm plantations formerly consisted of native species from local camel rangelands, such as Agathophora alopecuroides and Stipagrostis plumosa [START_REF] Ozenda P | Flore du Sahara, 3rd edition[END_REF].

Although the plantations were set in 1972, their establishment consisted of a single hydroagricultural development that disturbed the soil locally [START_REF] Zaid A | Date palm cultivation[END_REF]: the disturbance was limited to a few square meters around each planted tree which were spaced 11m apart. While palm tree plantations used to be well-drained, there is now a lack of maintenance of drainage systems, which combined with the water table rise, makes soils wetter. Also, the water from the water table is increasingly salty, and salt gradually accumulates in the soil over time [START_REF] Bouarfa | Salinity patterns in irrigation systems, a threat to be demystified, a constraint to be managed: Field evidence from Algeria and Tunisia[END_REF]. Therefore, with time, weed flora characterizing humid and salty soils and consisting of such species as Cressa cretica, Cyperus rotundus and Phragmites australis, has colonized the area.

Conclusions

The weed flora after 20 years was very different among the three studied crop types, mainly because their management sequences are different [START_REF] Salama Fm | Weed flora of common crops in desert reclaimed arable lands of southern Egypt[END_REF]. Weed monitoring over 20 years shows that weed composition is changing over time, particularly in cereal plots and in palm plantations. Whereas changes are likely driven by lack of drainage maintenance and salt accumulation in the case of palm plantations, more studies are needed to understand weed species composition in the cereal plots and greenhouses. This first study on crop weed composition in the Algerian Sahara is only a first step on the path to improved management, and an analysis of weed provenance (e.g., crop seed types and use of organic fertilization) and weed traits should be completed next.

Table 1. Year of surveys of the 100 palm plantations, 9 cereal plots and 50 greenhouses between 1990 and 2003. One presence-absence plant list was available per year, compiled from three annual surveys: to match the seasonal planting and harvest schedule, surveys were carried out starting in October of each year (x), and repeated in February and June of the following year (x + 1), to give a total of three surveys per year. The data from each set of three surveys were then taken together to fully characterize each October x to June x+1 crop year [START_REF] Itcmi | Rapport d'activité de recherche de la station expérimentale de Hassi Ben AbdAllah-Ouargla pour la période allant de 1990 à 2002[END_REF]

  murale, Lolium multiflorum, Melilotus indicus, Lysimachia arvensis, Sonchus oleraceus, Setaria verticillata, Polypogon monspeliensis and Erigeron canadensis, all of which are therophytes. Five of them belong to the complex chorotype (Chenopodium murale, Lysimachia arvensis, Polypogon monspeliensis, Setaria verticillata, Sonchus oleraceus), two are Mediterranean (Lolium multiflorum, Melilotus indicus), and one is an alien species from America (Erigeron canadensis).

  , axis 1 (23.9%) delineated old samples (1990-2000) correlated with Ammosperma cinereum, Atractylis delicatula, Brocchia cinerea, Diplotaxis acris, Launaea mucronata, Malcolmia africana and Paronychia arabica from other samples (Fig. 7). Axis 2 (13.7%) separated the samples gathered between 2003 and 2005, correlated with Callipeltis cucullaris, Lysimachia arvensis and Rostraria cristata from the more recent samples (2008 to 2010) correlated with Calendula arvensis, Frankenia pulverulenta, Galium tricornutum, Henophyton deserti, Papaver rhoeas, Scandix pecten-veneris, Sonchus asper, Stipagrostis obtusa and Tetraena alba (Fig. 7).

FigureFig. 1 .

 1 Figure legends are the same for the coloured online version and for the B&W version of the figures.

Fig. 2 .

 2 Fig. 2. Correspondence Analysis run on all data (282 samples × 162 species). Only the most correlated species are shown. Samples are grouped according to the crop types.

Fig. 3 .

 3 Fig. 3. Mean percentage of species / transect in each chorotype, considering crop types separately. For this analysis, only the data from 2003 to 2010 were used: cereal plots (sample unit =200m²), palm tree plantations (sample unit = 1000m²) and greenhouses (sample unit = 50m²). Values are shown as means ± SE. Bars with common letters do not differ significantly according to Tukey's pairwise comparison tests (p-value > 0.05).

Fig. 4 .

 4 Fig. 4. Mean percentage of species / transect in each life form, considering crop types separately. For this analysis, only the data from 2003 to 2010 were used: cereal plots (sample unit =200m²), palm tree plantations (sample unit = 1000m²) and greenhouses (sample unit = 50m²). Values are shown as means ± SE. Bars with common letters do not differ significantly according to Tukey's pairwise comparison tests (p-value > 0.05).

Fig. 5 .

 5 Fig. 5. Mean number of weed species in each year from 1990 to 2010 considering crop types separately: in the cereal plots (dot; sample unit =200m²), in the greenhouses (triangle; sample unit = 50m²) and in the palm plantations (cross; sample unit = 1000m²). Analyses could only be run on the data from 2003 to 2010. See results of analyses in Table 2. Values are shown as means ± SE. Points with common letters do not differ significantly according to Tukey's pairwise comparison tests (p-value > 0.05).

Fig. 6 .

 6 Fig. 6. Percentage of species in a particular chorotype in each year from 1990 to 2010, considering crop types separately: in the cereal plots (sample unit =200m²), in the greenhouses (sample unit = 50m²) and in the palm plantations (sample unit = 1000m²). Analyses were run on the data from 2003 to 2010, on the two most common chorotypes for each crop type. See results of analyses in Table 2. Values are shown as means ± SE. Points with common letters do not differ significantly according to Tukey's pairwise comparison tests (p-value > 0.05).

Fig. 7 .

 7 Fig. 7. Correspondence Analysis run on cereal plot data (71 samples × 104 species). Only the most correlated species are shown. Samples are grouped according to the year of study from 2003 to 2010. From 1990 to 2003, there is only one point (in grey) / year.

Fig. 8 .

 8 Fig. 8. Correspondence Analysis run on greenhouse data (139 samples × 31 species). Only the most correlated species are shown. Samples are grouped according to the year of study from 2003 to 2010. From 1990 to 2003, there is only one point (in grey) / year.

Fig. 9 .

 9 Fig. 9. Correspondence Analysis run on palm tree plantation data (72 samples × 119 species). Only the most correlated species are shown. Samples are grouped according to the year of study from 2003 to 2010. From 1990 to 2003, there is only one point (in grey) / year.

Table 2 .

 2 . Results of the Mixed-effects Generalized Linear Models analysing differences in numbers of weed species (species richness lines) or percentages of weed species (other lines) / transect along time. Each crop type was analyzed separately. Two models were compared for analysis: one with time included as an explanatory variable and one without. For these analyses, only the data from 2003 to 2010 were used. If the model with time had a lower AIC

		Cereal plots	Greenhouses	Palm plantations
	1990-91	×	×	×
	1991-92		×
	1992-93		×	×
	1993-94	×	×	×
	1994-95	×	×
	1995-96		×
	1996-97		×
	1997-98		×	×
	1998-99	×	×	×
	1999-00	×	×	×
	2000-01	×	×	×
	2001-02	×	×	×
	2002-03	×	×	×

Acknowledgements

We thank Mr Douadi Gousmi (Director of Benabdellah Experimental Station in Ouargla).

We also thank all the people who helped with our work, in particular Daniel Pavon (IMBE botanist) who provided invaluable assistance in identifying plants and seedlings and who helped us compiling plant information. We are grateful to the numerous dedicated reviewers who helped to greatly improve the manuscript.


Emmanuel J P Douzery
email: douzery@isem.univ-montp2.fr

Elizabeth A Snell

Eric Bapteste

Frédéric Delsuc

Hervé Philippe
email: philippe@umontreal.

The timing of eukaryotic evolution: does a relaxed molecular clock reconcile proteins and fossils?

The use of nucleotide and amino acid sequences allows improved understanding of the timing of evolutionary events of life on earth. Molecular estimates of divergence times are however controversial and are generally much more ancient than suggested by the fossil record. The limited number of genes and species explored, and pervasive variations in evolutionary rates are the most likely sources of such discrepancies. Here we compared concatenated amino acid sequences of 129 proteins from 36 eukaryotes to determine the divergence times of several major clades, including animals, fungi, plants, and various protists. Due to significant variations of their evolutionary rates and to handle the uncertainty of the fossil record, we used a Bayesian relaxed molecular clock simultaneously calibrated by 6 paleontological constraints. We show that, according to 95% credibility intervals, the eukaryotic kingdoms diversified 950-1,259 million years ago (Mya), animals diverged from choanoflagellates 761-957 Mya, and the debated age of the split between protostomes and deuterostomes occured 642-761 Mya. The divergence times appeared to be robust with respect to prior assumptions and paleontological calibrations. Interestingly, these relaxed clock time estimates are much more recent than those obtained under the assumption of a global molecular clock, yet bilaterian diversification appears to be about 100 million years more ancient than the Cambrian boundary.

Manuscript information: 30 text pages (including references and figure legends) ; two Figures ; two Tables ; five Supporting Informations (S1 to S5). Characters count:

All text characters plus spaces = 40,231 Fig. 1 (2-columns, 12 cm high = 360 x 12) = 4,320 Fig. 2 (1-column, 5 cm high = 180 x 5) = 900 Table 1 (1-column, 6 lines high = 60 x 6) = 360 Table 2 (

INTRODUCTION

Reconstructing and dating the early evolutionary history of eukaryotes has proven a challenging question due to the scarcity of the fossil record, especially for protists [START_REF] Knoll | Life on a young planet: the first three billion years of evolution on Earth[END_REF]. In the Archean eon, biomarker compounds characteristic of possibly extinct stem eukaryotes are found ~2,700 million years ago (Mya) [START_REF] Brocks | [END_REF]. Morphological clues for eukaryotes seem to occur in the Paleoproterozoic, but convincing microfossil evidence appears ~1,500 Mya in Lower Mesoproterozoic successions (3). Fossils that can be interpreted as crown eukaryotes are found in Upper Mesoproterozoic / Lower Neoproterozoic rocks, around 1,200 Mya for red algae (4), and 1,000 Mya for stramenopiles (5). Chlorophytes and testate amoebae are found in Middle to Late Neoproterozoic formations at ~750 Mya (5,6). Special attention has been attracted to the origin of the animal phyla. Bilaterian metazoans first appeared at the end of the Late Proterozoic (~600 Mya at the Doushantuo Formation) (7) and were already quite diversified ~530 Mya (8), suggesting a rapid diversification, known as the Cambrian explosion.

The advent of molecular data started a new era by providing an independent approach to address the history of life on Earth, and thus creating a closer connection between geology, paleontology, and biology (9). The constancy of molecular evolutionary rate over time, i.e. the molecular clock hypothesis, was suggested early on [START_REF] Zuckerkandl | Evolving Genes and Proteins[END_REF]. Calibration of this clock based on taxa with a rich fossil record allowed the inference of divergence times for other living organisms [START_REF] Bromham | [END_REF]. Yet, DNA and protein age estimates of land plants, fungi, and animal phyla possibly appeared several hundred million years (Myr) older than indicated by paleontology (12).

A major concern in molecular dating is that independent studies based on multiple genes yielded highly variable estimates. For example, the protostome / deuterostome divergence time ranges from 573 to 1,200 Mya, albeit either close to or drastically more ancient than the Cambrian explosion [START_REF] Feng | Proc. Natl. Acad. Sci[END_REF][START_REF] Ayala | Proc. Natl. Acad. Sci[END_REF][START_REF] Gu | [END_REF](16)(17)(18)(19)(20)(21)(22)[START_REF] Peterson | Proc. Natl. Acad. Sci[END_REF]. The observed discrepancy between clocks and rocks might be due to an extended period of Precambrian metazoan diversification without any trace in the fossil record [START_REF] Cooper | [END_REF], but also to biases in molecular dating methods. Among the limitations of the latter approaches, the non-clocklike behavior of genes and proteins (25), the limited power of the relative rate tests used to detect lineage-specific rate variation (26), and the existence of a systematic bias toward overestimation of evolutionary time scales (27) are the most confounding. For example, incorrectly assuming a clocklike property of the sequence data involves higher (22,28,29) or lower (30) divergence times. Consequently, dating methods have been developed to relax the molecular clock assumption by allowing rate variations along branches of a phylogenetic tree (28,(31)(32)(33).

To estimate the timing of eukaryotic evolution while reducing the impact of stochastic error due to the finite length of the sequences, we used a large data set of 129 combined nuclear proteins for 36 animals, fungi, plants, and protists. Despite the expectation that potential lineage-specific evolutionary rate variations from one protein to another might be counterbalanced in large-scale analysis (34), extensive among species variations in the rate of amino acid replacements were identified. We thus took advantage of a relaxed molecular clock approach developed in a Bayesian framework, which approximates divergence times by a Markov chain Monte Carlo (MCMC) numerical method (32,33). Such a relaxed clock method has three major advantages: (i) the biologically unverified hypothesis of a constant evolutionary rate all over the history of eukaryotes is unnecessary, (ii) the incorporation of prior constraints on divergence times is preferred over the use of fixed time points in order to handle the inherent uncertainties of paleontological data (fossils never match exactly on nodes of phylogenetic trees), and (iii) several independent calibrations can be used simultaneously. Since divergence times based on this relaxed protein clock were much more recent than those based on a global (constant) clock, computer simulations were used to explore the reliability of the Bayesian dating method.

MATERIAL & METHODS

The assembly of the data set of 129 proteins from 36 eukaryotes and its partitioned maximum likelihood (ML) analysis are provided elsewhere (35). Due to ambiguous alignment, two third of the sites were conservatively removed to ensure the sitewise homology of the remaining 30,399 positions at this deep evolutionary time scale. Alignments are available from HP upon request. We verified orthology through the analysis of each gene with standard phylogenetic methods, and we discarded several proteins, including some that were frequently used for studying eukaryotic evolution (tubulins, actin, and HSP70). Our reasonable taxon coverage in this large multiprotein framework involves that some genes are lacking for some taxa, leading to ~25% of missing character states. The best topology (Fig. 1) was rooted by Dictyostelium, as suggested by a derived gene fusion separating opisthokonts (animals, choanoflagellates, fungi) from bikonts (eukaryotes ancestrally with two cilia) (36,37). However, because of the potential homoplasy of gene fusion and fission characters (38), a more classical rooting along the kinetoplastid branch was also considered.

Relaxed molecular clock analyses

Divergence times were estimated under a Bayesian relaxed molecular clock (33), constrained by six independent fossil references. To explicitly incorporate the inherent uncertainty of paleontological estimates, each primary calibration was defined as a prior time constraint, with lower and / or upper bounds. In a conservative approach, we used the stratigraphic range of the geological period to which key fossils for the node under focus Myr at the ingroup root node, 0.001 (SD = 0.001) for the parameter (ν) that controls the degree of rate autocorrelation per Myr along the descending branches of the tree. The highest possible time between tips and root was 4,500 Mya. The MCMC were started from random values. After conservatively discarding 100,000 generations as the burn-in, they were run for 10 million generations. A sample was taken each 100 generations. Posterior divergence times were identical for different MCMC started independently, evidencing convergence of the Bayesian dating procedure. The uncertainty on divergence time estimates was characterized by 95% credibility intervals, calculated after sorting the 100,000 MCMC samples, and then reporting the 2,501 st and 97,500 th values.

Simulation studies

Simulation studies were conducted to evaluate (i) the ability of the global and relaxed clock approaches to accurately estimate divergence times when rates do correlate or not along tree branches, and (ii) the impact of missing data. We generated 10 matrices of 36 taxa and 30,399 amino acid positions (i.e., 1,094,364 character states) under PSeq-Gen [START_REF] Grassly | [END_REF], by using empirical amino acid frequencies, and a Gamma distribution parameter (α = 0.70) for rate heterogeneity among sites. Each branch length was set equal to the product of its rate (i.e., the average of the rates at the nodes that begin and end it) by its time duration (i.e., the age difference between its beginning and ending nodes) as measured by DIVTIME.

For each of these 10 simulated data sets, we applied the six calibration constraints, and computed the Bayesian divergence times as described in the previous section. Then, we removed each character corresponding to the missing data of the original data set (281,517 amino acids, i.e., ~25% of the data), and again calculated the divergence times. Age estimates computed with and without missing data were then compared for each node.

Moreover, a second set of 10 complete matrices was generated under a global clock model. The same topology, amino acid frequencies, Gamma parameter, and distribution of node times were used, but a single (constant) rate of 0.025 amino acid replacement per 100 sites was enforced for branch length entries. Autocorrelated-rate and constant-rate matrices were analyzed under ESTBRANCHES. Then, MCMC runs were conducted under DIVTIME to estimate divergence times, assuming either rate autocorrelation along branches (mean and SD for ν ≠ 0) or a perfect clock (mean and SD for ν = 0).

RESULTS

A molecular time scale for the main eukaryote clades

The phylogeny resulting from the ML analysis of the 30,399 unambiguously aligned amino acid positions (Fig. 1) is in agreement with current opinion about eukaryotic relationships (44). The log-likelihoods (lnL) of the phylogenetic tree were lnL UNCONSTRAINED = -779,284 without a molecular clock constraint, and lnL CLOCK = -783,020 under the global clock and a rooting by Dictyostelium. A likelihood ratio test significantly rejected the hypothesis of a clocklike behavior of our data (P < 0.0001). Thus, the rate of amino acid replacements extensively changed among the eukaryote lineages here analyzed. For example, trypanosomes and nematodes evolved two times and three times faster than mammals, respectively (data not shown).

The relaxed molecular clock based on the 129 proteins was calibrated by six time constraints spread over the phylogenetic tree and defined within green plants, animals, and fungi. Credibility intervals at 95% indicated that the basal split between the major eukaryotic kingdoms documented here occurred 950-1,259 Mya (mean: 1,085), followed by their diversification in less than 200 Myr (Fig. 1). Plantae originated 892-1,162 Mya (mean: 1,010), and red algae branched off 825-1,061 Mya (mean: 928) whereas stem land plants separated from green algae 662-812 Mya (mean: 729). This suggests that the endosymbiosis of a freeliving cyanobacterium that led to primary plastids (surrounded by two membranes) occurred between 825 to 1,162 Mya (Fig. 1). Plastids surrounded by four membranes originated from secondary endosymbiosis, in which a red alga was engulfed and retained by a flagellate protist. This event probably happened along the branch leading to stramenopiles + alveolates, i.e., between 767-1,072 Mya, shortly after the primary endosymbiosis (Fig. 1).

Among opisthokonts, most of the inferred divergence times were more recent than previous molecular multigene estimates. The divergence between animals and fungi was for example estimated 872-1,127 Mya (mean: 983 Mya) whereas others provide estimates up to 1,600 Mya (45). The divergence between choanoflagellates and animals was estimated during the Late Proterozoic, between 761 and 957 Mya (mean: 849). The major transition during animal evolution from unicellular to multicellular state might have thus occurred at least 200 Myr before the Cambrian explosion. Interestingly, the divergence time between protostomes and deuterostomes (Fig. 1) was 642-761 Mya (mean: 695). This value is much lower than most of the previously published ones (12,17,46), but is close to recent estimates (22,[START_REF] Peterson | Proc. Natl. Acad. Sci[END_REF], reducing to about 40 Myr the minimum gap between the protein relaxed clock timing and the first indisputable bilaterian fossil (7). Within fungi, the divergence between ascomycetes and basidiomycetes is estimated around 629-837 Mya (mean: 727), whereas a recent study using a global clock method on a similar number of proteins, but a lower number of taxa and calibration points, yielded 1,200 Mya (45).

Global versus relaxed clock estimates

Our molecular time scale for the major eukaryotic clades is on average 1.5 times more recent than the one based on a global clock (12). To evaluate the accuracy and precision of dates estimated under linear versus autocorrelated rates extrapolations beyond calibrations, we conducted simulations. 1).

Impact of missing data on divergence times

Despite the expectation that the large size of our dataset would have likely buffered the impact of missing data (47), our dating estimates might have been seriously influenced by the absence of about 25% of amino acids. This point was evaluated on simulated protein sequences, by comparing divergence times inferred from complete and gapped alignments (Supporting Information S1). The presence of missing data leads to slightly under-or overestimated divergence ages, depending on the node under focus, and induce a mean difference of -0.2 % ± 2.8 over all nodes (extreme values: -9.8 to +8.7%). On average, the relative difference is respectively +1.7 % and -2.5 % for ages which are over-and underestimated by missing data. This corresponds to variations of ca. +17 Myr to -25 Myr for a 1,000 Myr old node. For the root age, these values are indeed twice less than the actual standard deviation measured by the Bayesian approach. Interestingly, the random introduction of an additional 25% of missing data-yielding a data set with 50% of missing amino acids-again provides very similar dates relative to the original data set (data not shown). The relaxed clock method is therefore efficient, even when the alignment contains partial sequences.

Sensitivity of posterior divergence ages upon prior assumptions

Since Bayesian posteriors might depend on priors, we evaluated the sensitivity of our dating results upon a priori hypotheses. Different expected number of time units between tree root and tips were tested by comparing posterior ages obtained using different Gammadistributed priors: 540 ± 270 Mya (the Cambrian boundary), 1,000 ± 500 Mya (earliest body fossils of stramenopiles ( 5)), 2,000 ± 1,000 Mya (a recent molecular estimate for eukaryotes origin ( 12)), and 3,000 ± 1,500 Mya (a deep value for stem eukaryotes (2)). Posterior ages were plotted against the four independent root priors for each of the 34 tree nodes. This resulted in 34 regression lines with a mean slope of 0.016 (range: -0.002 ± 0.0005 to 0.067 ± 0.007). Such an average impact over all nodes indicated that a 1,000 Myr variation of the root prior will involve a posterior divergence age variation of only 16 Myr, which is far less than the order of magnitude of the uncertainty on the divergence ages themselves. Posterior divergence times thus appeared only slightly sensitive to huge changes in prior specification of root age. Furthermore, comparison of a priori and a posteriori age estimates for calibration nodes did not reveal any systematic trend (Supporting Information S2), suggesting that paleontological uncertainties did not bias our calculations since they were incorporated in the dating procedure (31,33).

Changes in the reference topology may potentially affect molecular dating (39).

Since the root location is uncertain among eukaryotes (44), we recalculated divergence times after moving the root from the branch leading to the cellular slime mold to the one leading to kinetoplastids (Fig. 1). Divergence times in opisthokonts, especially among animals, were not sensitive to this modification (e.g., we obtained 627-738 Mya for bilaterians). The kinetoplastid rooting actually involved deeper divergences among Plantae (899-1,191 Mya).

The main clades of eukaryotes also diverged earlier (1,042-1,412 Mya: Supporting Information S3), likely because such an alternative root location involves a pectinate tree shape. Thus, we confirm that the above-mentioned observations of Middle-to Late-Proterozoic cladogeneses among most eukaryotes were not the result of a rooting artifact.

Molecular clock approaches are also sensitive to the choice of calibrations (48). For example, divergence time estimates under the single vertebrate calibration were on average 1.4 times deeper than those measured with the six calibrations (Supporting Information S4).

We systematically studied the effect of calibration selection by analyzing all the possibilities of choosing one or combining 2, 3, 4, 5 or 6 calibrations. Whereas divergence times averaged over all possibilities of using a given number of simultaneous calibrations seems not to vary much with those numbers, the minimum and maximum estimates may dramatically vary depending upon the choice of a given combination of fossil references (Table 2). However, a reduction of the uncertainty on divergence times was clearly associated to the increase of the number of simultaneous calibrations. Importantly, our datings were robust against the random exclusion of one or two calibration points.

DISCUSSION

Robustness of the relaxed clock estimates

Molecular clocks are highly controversial, dates from various sources being often so different that some researchers doubt that molecular dating of evolutionary events is possible (49). Three points are of crucial importance: (i) the size of the molecular sample (i.e.

stochastic errors due to a limited number of species and/or genes), (ii) the efficiency of the dating method (i.e. systematic errors in the inference), and (iii) the quality of the fossil calibrations (i.e. amplification of paleontological uncertainties). We tried to reduce the impact of all these potential issues.

Our divergence time estimates should not be seriously affected by the stochastic error since we used a very large data set (129 genes, 30,399 amino acid positions) corresponding to proteins involved in different functional pathways like translation, transcription, cytoskeleton, and metabolism. Contrary to all previous studies with such a large sample of genes, a reasonable taxonomic diversity (36 species) was sampled. Here, half of the total number of proteins contains at most three missing taxa over 36, and 25% of the amino acids are missing in our final alignment. Nevertheless, our dating estimates are robust upon missing data (Supporting Information S1). Indeed, a potentially reduced phylogenetic (and/or dating) accuracy associated with missing data is primarily caused by the occurrence of too few complete characters rather than too many missing positions (47).

Presently, since the shortest concatenated sequence includes 7,500 amino acids, our complete alignment retains a strong phylogenetic and dating signal.

Second, our relaxed molecular clock datings (Fig. 1) appear to be robust to the specification of prior distributions on model parameters. This suggests that much posterior dating information regarding measure of divergence times among eukaryotes is attributable to the 129 concatenated proteins rather than to prior specifications. Moreover, the relaxed clock approach retrieves the true divergence times (Table 1) when substitution rates are distributed according to the autocorrelated-rate model developed by Thorne et al. (32).

Third, the concomitant use of six calibration constraints chosen among plants, fungi, and animals (Figure 1) also helped to more accurately measure the absolute amino acid replacement rates in these three distinct clades, and stabilized the Bayesian divergence times on the whole tree. An alternative procedure (22,28) is to use several independent calibrations, and, for each gene, to average the corresponding divergence estimates.

However, such an approach is "inferior to a simultaneous analysis of multiple genes under the constraints of multiple calibrations" (22). Accordingly, we have shown that the use of simultaneous calibrations reduces the uncertainty on age estimates (Table 2).

nuclear genes (22). Second, simulation studies (Table 1) indicate that, at least in the present case, the autocorrelated rates model does not underestimate divergence times, neither on clocklike nor on relaxed clock sequences.

Comparison with other relaxed molecular clock datings

Recently, two relaxed clock dating studies-among metazoans [START_REF] Peterson | Proc. Natl. Acad. Sci[END_REF] and photosynthetic eukaryotes (50)-yielded results that are, at first sight, incongruent with our estimates. The use of seven concatenated proteins and the non parametric rate smoothing r8s approach ( 31) indicate a more recent occurrence of the last common bilaterian ancestor, between 556-592 Mya ( 23) instead of 642-761 here. Importantly, these estimates were obtained without taking into account among site rate variation. In the same paper, calculations conducted by modeling such a rate variation with a Gamma distribution yield a confidence interval of 636-678 Mya for the bilaterian diversification, i.e., an estimate overlapping with ours.

On the basis of a tree derived from 5 concatenated plastid proteins, calibrated by 5 intervals under the r8s approach, Yoon et al. (50) suggested a red-green algae split between 1,591-1,757 Mya, instead of 825-1,061 or 899-1,191 Mya here depending on the root location. Our reanalysis of this data set under the DIVTIME relaxed molecular clock approach suggests a major incompatibility between "green" and "red" calibrations. Under the ML topology inferred from the 5 plastid proteins, and enforcing the angiosperms, spermatophytes, and land plants nodes between 90-130, 290-320, and 432-476 Mya (50) would involve a red algae diversification between 392-761 Mya (mean 551 Ma), and a

Chondrus / Palmaria split between 165-381 Mya. This is incompatible with the respective (4,50). More generally, molecular dating still needs to be refined by improved sampling of genes and especially of species, and improved methods with more realistic models of rate variation.

Towards a partial reconciliation of clocks and rocks

Deeper molecular dates are actually expected since the fossil record documents only the first appearance of a morphologically recognizable (crown) group, and not the actual time of genetic divergence (51). Previous molecular datings yielded molecular ages always more ancient than paleontological ones, but often to a large extent (12,[START_REF] Gu | [END_REF]20,46). Interestingly, our molecular dating (Fig. 1) does not appear to be biased in such a systematic way.

Molecular ages are more recent than the fossil record for red algae but more ancient for animals. We dated the divergence between green plants and red algae between 825-1,061

Mya (Fig. 1) whereas a fossil of 1,200 Mya has been interpreted as belonging to a specific subgroup of red algae, i.e. bangiophytes (4). Molecular dating with a kinetoplastid rooting would be in better agreement with this fossil (95% credibility interval: 899-1191 ; Fig. 1), but such a rooting might represent a tree reconstruction artifact (36). Actually, single fossils may also bear their load of uncertainty in taxonomy identification and geological ages [START_REF] Bromham | [END_REF].

Illustrating the latter case, the age of the Bangiomorpha assemblage is in fact bounded between 1,267 and 723 Mya (52). For bilaterians, the discrepancy between our molecular estimates (642-761 Mya) and the fossil record ( 7) is reduced to a minimum of 40 Myr.

According to protein data, the divergence of the main stem groups of metazoans (Deuterostomia, Ecdysozoa and Lophotrochozoa) would have occurred in the Late Neoproterozoic, but their diversification would date close to the Precambrian-Cambrian boundary. This suggests that the stem lineages would have been organisms with a poor fossil record because of their small size, taphonomic problems, or undistinguishable morphology.

The use of a relaxed molecular clock on our large data set estimated the diversification time of most extant eukaryotic kingdoms at ~1,100 Mya in the mid-Proterozoic (Fig. 1). Even if cytoskeletal and ecological prerequisites for eukaryote diversification were already established some 1,500 Mya (3), the postulated anoxic and sulfidic redox status of oceans between ~2,000 to 1,000 Mya might have limited the rise of photosynthetic eukaryotes (53). Around 1,000 Mya, the diversification of plants may have been facilitated by the primary endosymbiotic event that led to plastids. Moreover, it has been proposed that a major increase in atmospheric oxygen level occurred between 1,000 and 640 Mya (54). In association with the Neoproterozoic appearance of more fully oxic oceans (53), this major transition towards a more oxygenic environment could possibly have accelerated the diversification of the eukaryotic organisms, thanks to the aerobic respiration provided by mitochondria [START_REF] Philippe | Evolutionary relationships among Protozoa[END_REF]. Close to the Meso-/ Neo-Proterozoic transition, choanoflagellates diverged from their closest animal relatives marking the dawn of cellular cooperation in animals. The major eukaryotic kingdoms therefore originated much more recently than previously thought (12), but did not seem to have been affected by the postulated global glaciations ("snowball Earth" hypothesis) between 750 and 600 Mya [START_REF] Lubick | [END_REF]. The present study illustrates how geological and biological records of organismal evolution might be partially reconciliated by using molecular relaxed clock datings from large genomic comparisons.

  1-column, 8 lines high = 60 x 8

Figure 1 .

 1 Figure 1. Divergence time estimates (in million years ago [Mya]) among

Figure 2 .

 2 Figure 2. Relationship between node rates and node times. For each node of

  Table 1 recapitulates the age estimates of two selected nodes-bilaterians and our tree root-calculated under either global or relaxed clock assumptions onto data simulated with or without constant rate of evolution. Relaxed clock timings were accurate, though slightly deeper when estimated on constant-rate data.

They were also less precise-as attested by 95% credibility interval widths-relative to global clock estimates, because of the greater number of parameters involved

(22)

. However, constant-rate extrapolations conducted with all other conditions being identical yield less accurate, i.e., 1.2 to 1.4 deeper divergence times when actual rates do vary. The simulations therefore indicate that the autocorrelation model of substitution rate evolution does not underestimate ages on clocklike data, whereas a linear extrapolation on rate-variable data markedly overestimates divergence times (Table

Table 1 . Relaxed and global clock average age estimates (in Mya) for bilaterians and the tree root

 1 . Divergence times are computed on 30,399 amino acid long data for 36 taxa, simulated under variable (autocorrelated) or constant rates of evolution, by using either a relaxed or a global molecular clock. Expected (simulated) ages are 695 Mya for bilaterians, and 1,085 Mya for the root. The 95% credibility interval widths are given in parentheses. All values are averages given with ± the standard-error on 10 replicates.

		CLOCK	SIMULATIONS
			Variable rates	Constant rate
	BILATERIA	Relaxed	692 ± 6 (119 ± 4)	706 ± 7 (52 ± 6)
		Global	864 ± 9 (22 ± 0.4)	700 ± 7 (35 ± 5)
	ROOT	Relaxed	1,083 ± 15 (308 ± 14)	1,122 ± 11 (89 ± 8)
		Global	1,541 ± 14 (41 ± 1)	1,079 ± 8 (54 ± 7)

Table 2 . Posterior divergence ages for bilaterians and tree root computed on real data under a relaxed molecular clock calibrated by different combinations of one to six time constraints.

 2 Mean ages, range of ages, and mean 95% credibility interval widths (δ CredI) are computed on the 6, 15, 20, 15, 6, or 1 possibilities of incorporation of respectively C = 1, 2, 3, 4, 5, or 6 simultaneous calibrations.

		BILATERIANS	ROOT
	C	Ages (Range)	δ CredI	Ages (Range)	δ CredI
	1	694 (452-999)	330 ± 201	1,140 (751-1,617)	578 ± 307
	2	708 (536-999)	199 ± 97	1,147 (885-1,621)	379 ± 135
	3	709 (560-863)	157 ± 73	1,132 (908-1,361)	329 ± 82
	4	706 (642-845)	136 ± 56	1,117 (1,029-1,298)	312 ± 68
	5	699 (669-757)	123 ± 32	1,099 (1,045-1,174)	301 ± 39
	6	695 (-)	119	1,085 (-)	309

ACKNOWLEDGEMENTS

We thank Henner Brinkmann, Wilfried de Jong, Christophe Douady, Eric Fontanillas, Peter Holland, Emmanuelle Javaux, Franz Lang, Nicolas Lartillot, David Moreira, and two anonymous referees for helpful suggestions. Debashish Bhattacharya kindly gave us the alignment of 5 plastid genes. This work has been supported by the "ACI Informatique-

Variable-rate and constant-rate molecular clocks

The greatest discrepancy between age estimates without and with the clock assumption is observed for the Saccharomyces / Candida split-here estimated at 168-308 (mean: 235) Mya (Fig. 1) compared to 841 Mya (45). These two fungi appear to be fast evolving (data not shown), and such high rates are probably not accommodated by global clock based methods. When we constrained a perfect molecular clock under the Bayesian approach, this divergence was estimated at 439-468 Mya (mean: 454). On average, the divergence times among eukaryotes were 1.7 times deeper with a global clock than with a relaxed clock (Supporting Information S5). The greatest contrast (ratio > 2.0) between global and relaxed clock ages is observed for alveolates, nematodes, and kinetoplastids, i.e., the fastest-evolving taxa for the 129 proteins here evaluated. A similar result was obtained for animals under various models of clock relaxation (22,28), and for mitochondrial genomes of mammals (29).

How to evaluate which method, either a global or a relaxed clock, is better describing the data at hands? Estimating a time scale by the mean of age extrapolation beyond animal, fungi, and plant calibrations may possibly lead to biased estimates. Whereas interpolation of date estimates are certainly less biased because evolutionary rates are bounded by zero, extrapolation could conduct to rates increasing on average when moving far back in time. Thus, increased rates-faster earlier and slower since the Neoproterozoic / Cambrian transition-may result into apparently too recent estimates for deeper nodes. Two elements suggest that our analyses have not been affected by such a trend. First, we plotted the posterior evolutionary rate of each node against its mean posterior age estimate as measured on our concatenated data set (Fig. 2). No trend towards a rate increase far back in time is observed, and the rate for the last common bilaterian ancestor is close to the mean rate measured along the eukaryotic phylogeny. In particular, our 129 proteins did not record the Late Precambrian episodic burst of evolutionary rates detected on 22 mitochondrial and

SUPPORTING INFORMATION S1. Sensitivity of Bayesian dating estimates to missing data.

The histogram of the divergence time estimates of the nodes computed with and without missing data is reported. The percentage of variation of divergence ages after introducing 25% random gaps is given, and defined as 100 x [ 1 -(age of a given node computed with missing data) / (age of the same node computed with complete data) ]. It has been computed for 280 nodes = 28 nodes that are not under paleontological constraints x 10 replicates of simulations. Interesting to note, the greatest difference introduced by missing data is always less than 10%, meaning that the potential error on the age of the eukaryotic root due to missing data is inferior to ~100 Myr.

SUPPORTING INFORMATION S2. Divergence time estimates under a Bayesian relaxed molecular clock for the six calibration nodes.

The age of the six calibration nodes were a priori forced to fall into a certain time range, and after the MCMC analysis, posterior estimates of these calibration nodes are reported. Some posterior divergence times lay in the middle of the prior time interval (insects, arthropods), whereas others lay close to one end, either the recent (land plants, deuterostomes) or ancient (flowering plants) bound. Examination of the a posteriori calibrations thus did not reveal a systematic trend of posterior estimates to be bound upwards or downwards. The reference topology was rooted by Dictyostelium. Geological times containing the posterior estimates are also given.

The posterior Gamma distributions used to relax the molecular clock assumption were characterized by the following mean values, ± 1 standard deviation: ν = 0.00099 ± 0.00030 (prior: 0.001 ± 0.001) for the parameter that controls the degree of rate autocorrelation per million years along the descending branches of the tree, 1,085 ± 79 Mya (prior: 1,000 ± 500 Mya) for the root age, and 0.03775 ± 0.00804 (prior: 0.03641 ± 0.0182) amino acid replacements per 100 sites per million years at the ingroup root.

All absolute ages of the geological periods are taken from the 1999 Geological Time Scale of the Geological Society of America (www.geosociety.org/science/timescale/timescl.pdf).

Calibration Prior

SUPPORTING INFORMATION S3

.

Divergence time estimates (in million years ago [Mya]

) for the major groups of eukaryotes, based on a relaxed molecular clock applied to 30,399 amino acid positions, and calibrated by the simultaneous use of six independent paleontological references. The Kinetoplastidae rooted the tree (leading to a sister-group relationship between Stramenopiles + Alveolates, and a clade containing Plantae associated to Dictyostelium + Opisthokonta), but this outgroup was pruned from dating analyses. See the caption of Figure 1 for the conventions used.

SUPPORTING INFORMATION S4. Molecular divergence time estimates of some eukaryote groups obtained under a Bayesian relaxed molecular clock calibrated by either a single vertebrate constraint (Actinopterygii vs. Mammalia) or the six simultaneous constraints, and comparison with the global clock calibrated by the single vertebrate point.

The reference topology was rooted by the amoeba Dictyostelium, and the ingroup topology matched the one depicted on Figure 1. Mean posterior divergence times are given in million years ago, and 95% credibility intervals are provided between parentheses. On the 29 nodes not constrained by prior calibrations, divergence ages based on a single calibration within the slow-evolving vertebrates are 1.4 times deeper than estimates based on the synergetic use of the six calibrations (listed in SUPPORTING INFORMATION S2). This is an indication of the importance of checking dating estimates against multiple independent calibrations. Combining the global clock with the vertebrate calibration involves on average 2.9 deeper age estimates relative to the timings derived from the relaxed clock calibrated by six points. The reference topology was rooted by the amoeba Dictyostelium, and the ingroup topology matched the one depicted on Figure 1. Mean posterior divergence times are given in million years ago, and 95% credibility intervals are provided between parentheses. On the 28 nodes not constrained by prior calibrations, global clock estimates are 1.7 times deeper than relaxed clock estimates.

Number of calibrations

Relaxed