Joana Patrício

João Magalhães Neto

Heliana Teixeira

Fuensanta Salas

João Carlos Marques

The robustness of ecological indicators to detect long-term changes in the macrobenthos of estuarine systems

Keywords: ecological indicators, benthic macrofauna, intertidal, ecological status, estuaries, Portugal

Accurate and reliable benthic quality indicators are in great demand following the recent developments and the strict time schedule for implementing the European Water Framework Directive. The Mondego estuary has experienced a progressive deterioration during the 1990's, followed by a partial ecological recovery due to restoration measures in 1997/1998. We have used the estuary as a model system to test the performance and robustness of a set of ecological indicators in highlighting the changes in the ecological state of intertidal areas. Over a period of 17 years (1985-2002), we calculated Margalef, Shannon-Wiener, Berger-Parker, Taxonomic Distinctness measures, AZTI's Marine Biotic Index, Infaunal Trophic Index, and Eco-Exergy based indices and tested differences across periods characterised by different anthropogenic disturbance. We combined temporal data within three periods: before, during and after disturbance, based on ACCEPTED MANUSCRIPT 2 progressive information on the changes in the extended type of anthropogenic disturbance. Indices were then compared with biological and abiotic descriptors (macroalgae, macrophytes, benthic macrofauna, nutrients concentration, sediment grain size and total organic carbon). We found great disparity in the indicators ability to capture temporal changes, showing distinct performances at each site. At the Z. noltii site, only Margalef, Total Taxonomic Distinctness and the thermodynamically based indices captured temporal changes, despite giving higher values during the disturbance period. At the bare sediment site, Taxonomic Distinctness, ITI, Shannon-Wiener, Berger-Parker, AMBI and the TBI were able to distinguish between periods, in agreement with the differences observed analysing the macrobenthic assemblages. Furthermore, Taxonomic Distinctness was not robust enough to detect any temporal or spatial change. We thus suggest further research to understand the behaviour of ecological indicators, in view of their crucial importance for the management and protection of marine coastal areas.

Introduction

Background

Ecological indicators are commonly used to provide synoptic information about the state of ecosystems. Most often they summarise information on the structure and/or functioning of ecosystems using a particular biological aspect or component of these ecosystems, such as diversity, productivity or environmental conditions [START_REF] Salas | Ecological indicators in coastal and estuarine environmental quality assessment. A user-friendly guide for practitioners[END_REF]. In general, the application of ecological indicators is not exempt from criticism.

Aggregation of different metrics results in oversimplification of the ecosystem under observation. Moreover, the choice of which ecological indicators to use should follow criteria based on the question we want to address and the theoretical background on which the indicator was created. Before being able to use an indicator it is thus necessary to establish its applicability and robustness for a set of conditions, by testing its consistency with the ecological expectation.

The implementation of the European Water Framework Directive (WFD) (EC, 2000), as well as the recent European Marine Strategy Directive 2008/56/EC (EC, 2008) require the use of well-founded ecological indicators, which has contributed greatly to develop this research field and stimulate analyses to establish indicators' applicability to a variety of ecosystems. Indeed, since 2000, numerous ecological quality assessment tools have been proposed by the European scientific community (e.g. [START_REF] Borja | Implementation of the European water framework directive from the Basque country (northern Spain): a methodological approach[END_REF][START_REF] Rosenberg | Marine quality assessment by use of benthic species-abundance distribution: a proposed new protocol within the European Union Water Framework Directive[END_REF][START_REF] Bald | Assessing reference conditions and physico-chemical status according to the European Water Framework Directive: A casestudy from the Basque Country (Northern Spain)[END_REF][START_REF] Simboura | A synthesis of the biological quality elements for the implementation of the European Water Framework Directive in the Mediterranean ecoregion: The case of Saronikos Gulf[END_REF][START_REF] Devlin | Setting nutrient thresholds to support an ecological assessment based on nutrient enrichment, potential primary production and undesirable disturbance[END_REF]Muxica et al., 2007;[START_REF] Pinto | Review and evaluation of estuarine biotic indices to assess benthic condition[END_REF]. The WFD suggests that the water body ecological quality should be primarily determined by biological quality elements and then by hydromorphological and physicochemical quality elements.

The composition and structure of benthic macrofauna is one of the indicated biological quality elements to be used in transitional and coastal waters for quality status assessment. Benthic communities may detect the level of stress locally and integrate the recent history of stress, constituting a sort of memory for the system. Macrofauna are dominated by species with different mobility, life-cycle and tolerance to stress [START_REF] Dauer | Abundance-biomass comparison ABC method: effects of an estuarine gradient, anoxic/hypoxic events and contaminated sediment[END_REF], which covers the WFD demand of integrating differently sensitive species. The response of macrobenthic communities to several types of stress is well studied, based on multivariate analyses that take into account variations in species diversity and their relative abundance between perturbed and control sites (e.g. [START_REF] Pearson | Macrobenthic succession in relation to organic enrichment and pollution of the marine environment[END_REF][START_REF] Warwick | Comparing the severity of disturbance: a meta-analysis of marine macrobenthic community data[END_REF], Gray 2002). Based on this knowledge it is possible to determine a priori the ecological indicator behaviour and thus its appropriateness to detect changes in the variable of interest. The aim of our study was to test the performance and robustness of a large set of ecological indicators based on the composition and structure of macrobenthic communities.

A plethora of methodologies with hundreds of indices, metrics and evaluation tools are presently available [START_REF] Borja | Assessing the environmental quality status in estuarine and coastal systems: comparing methodologies and índices[END_REF]; e.g. see summary in Diaz et al., 2004). An ecologically parsimonious approach dictates that investigation should place greater emphasis on evaluating the suitability of indices that already exist prior to developing new ones (Díaz et al., 2004). We thus selected a large set of ecological indicators of different types (Margalef, Shannon-Wiener, Berger-Parker, Taxonomic Distinctness measures, AZTI's Marine Biotic Index (AMBI), Word Infaunal Trophic Index (ITI), and Eco-Exergy based indices) and compared the status of two sites "Z. noltii site" and a "bare sediment" in the Mondego estuary (western coast of Portugal; 40º08'N, 8º50'W) over a period of seventeen years. During this period, the estuary experienced changes in the type and severity of anthropogenic stress, which changed macrophyte distribution and macrofauna communities. We asked whether different ecological indicators were able to detect changes and whether these changes were in agreement with the differences observed analysing species composition of macrobenthic assemblages. Specifically, we tested for differences in the macrofauna composition and biological indicators between the periods during and after the major anthropogenic changes. We also assessed which ecological indicators were more robust in detecting these changes.

Synopsis of human induced changes in the Mondego estuary ecosystem

The Mondego estuary is a relatively small (860 ha) warm-temperate intertidal system. At present, the entire catchment area contributes a large loading of nutrients and several chemicals into the system. Moreover, the estuary supports industrial activities, saltextraction, and aquaculture farms, and is located close to a tourism centre (Figueira da Foz). The estuary's terminal part is 7 km long and is 2-3 km across at its widest part, consisting of two arms, separated by the Murraceira Island (Fig. 1). The northern arm is deeper (5-10 m during high tide), constituting the main navigation channel and the harbour location. The southern arm is shallower (2-4 m during high tide), and is characterised by large areas of intertidal flats exposed during low tide (almost 75% of the total area of the southern arm).

A detailed description of the system can be found in [START_REF] Cardoso | Dynamic changes in seagrass assemblages under eutrophication and implications for recovery[END_REF] and Marques et al. (2007). A comprehensive study on the environmental quality of the Mondego estuary was carried out during the last two decades [START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF] and covered several quality elements (water quality, hydraulics and sediments dynamics, benthic intertidal and subtidal communities, Z. noltii beds and green macroalgal cover). Such analyses allowed detection increased eutrophication symptoms and ecological deterioration in the South arm since early 1990s, when the communication between the two arms was interrupted, increasing water residence time and nutrient concentration [START_REF] Lillebø | Management of a shallow temperate estuary to control eutrophication: the effect of hydrodynamics on the system nutrient loading[END_REF].

Since then, macroalgal blooms became more frequent and intense (mostly Ulva spp.) (e.g. Martins et al., 2001;[START_REF] Cardoso | Dynamic changes in seagrass assemblages under eutrophication and implications for recovery[END_REF][START_REF] Pardal | Assessing environmental quality: a novel approach[END_REF][START_REF] Leston | The response of primary producer assemblages to mitigation measures to reduce eutrophication in a temperate estuary[END_REF] and have affected Z. noltii beds, reducing the area covered by the seagrass and the bodiversity and community biomass in the remaining meadow [START_REF] Cardoso | The impact of extreme flooding events and anthropogenic stressors on the macrobenthic communities' Dynamics[END_REF][START_REF] Lillebø | Management of a shallow temperate estuary to control eutrophication: the effect of hydrodynamics on the system nutrient loading[END_REF]. Before the interruption, the Z. noltii covered 15 ha (in 1986), which decreased to 1.6 ha in 1993, when macroalgal blooms increased up to 415 g m 2 (Fig. 2). In 1997 Z.

noltii bed reached its lowest value (0.02 ha). Through time, such modifications promoted the emergence of a simplified trophic structure, which has been analysed in the existing literature (e.g. Dolbeth et al., 2007;[START_REF] Martins | Seasonal variation of short-term survival of Zostera noltii transplants in a declining meadow in Portugal[END_REF][START_REF] Lopes | Influence of macroalgal mats on abundance and distribution of dunlin Calidris alpina in estuaries: a long termapproach[END_REF][START_REF] Patrício | Mass balanced models of the food web in three areas along a gradient of eutrophication symptoms in the south arm of the Mondego estuary (Portugal)[END_REF]. Since 1998, the system has started recovering, after the decrease in freshwater discharge from the Pranto River and the limited re-establishment of communication between the two arms [START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF][START_REF] Leston | The response of primary producer assemblages to mitigation measures to reduce eutrophication in a temperate estuary[END_REF][START_REF] Cardoso | The impact of extreme flooding events and anthropogenic stressors on the macrobenthic communities' Dynamics[END_REF]. Z. noltii meadow has began to recover and in 2005, it occupied an area of about 4.2 ha.

Material and Methods

Description of the database used

Intertidal benthic communities were first sampled in 1985 (during winter) and 1986 (during summer). Subsequently, the community was regularly monitored from January 1993 to January 1996, and again from January 1999 to December 2002 once every two weeks from 1993 to 1996, and monthly from 1999 to 2002.

In all cases, sampling was conducted at two sites in the South arm (Fig. 1), within Z. noltii meadow, which persisted during early 1990s, and on an area where the Z. noltii has disappeared, now characterised by bare sediment. Each time and at each site, six replicate cores were taken to a depth of 20 cm (13 cm inner diameter and then pooled together for macrofauna analyses. Sediment was washed through a 500 µm mesh sieve and animals preserved in 4% buffered formalin. Animals were identified to species level and subsequently dried at 70ºC for 72h to estimate biomass as dry weight (DW) and ash free dry weight (AFDW), after combusting samples for 8h at 450ºC.

Hydrological changes and variation in macroalgal and macrophyte coverage and biomass

Hydrologic changes in the Mondego estuary, following human interventions, have been assessed using technical reports, maps and aerial photographs. Samples of macrophytes and macroalgae were taken once every two weeks from January 1993 to January 1996, and subsequently, from January 1999 to December 2002, with a monthly periodicity. Six cores (13 cm diameter to a depth of 20 cm) were randomly taken from each site (Z. noltii site and bare sediment site), during low tide. Macroalgal and macrophyte species were identified to the species' level, dried at 70ºC for 48 hours and then combusted at 450ºC for 8 hours to estimate biomass as AFDW. Additionally, the area of Z. noltii meadows was measured through aerial photographs and field campaigns [START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF].

Ecological indicators

Description and calculation

The following ecological indicators were calculated using the macrofauna database, for each site and time. Replicates within site and time were pooled.

1. The AZTI's Marine Biotic Index, AMBI, [START_REF] Borja | Marine Biotic Index to establish the ecological quality of softbottom benthos within European estuarine and coastal environments[END_REF] is based on the proportion of five ecological groups to which the benthic species are allocated. The soft bottom macrofauna was divided into five groups according to their sensitivity as a function of an increasing stress (see Table 1 for details). The AMBI 4.0 software, freely available at http://www.azti.es, was used to calculate this index, using a species list of July 2006 and following [START_REF] Borja | Guidelines for the use of AMBI (AZTI's Marine Biotic Index) in the assessment of the benthic ecological quality[END_REF] guidelines. AMBI ranges from 0 to 7, and a better quality trend would be indicated by a decrease on the index value (see Table 1). 2. To estimate Taxonomic Diversity indices -Taxonomic Diversity, ∆ ; Taxonomic Distinctness, * ∆ and Total Taxonomic Distinctness, s + -a hierarchical Linnean classification was used as a proxy for cladograms representing the relatedness of individual species. For each location, a composite taxonomy was compiled, five taxonomic levels were considered (species, genus, family, order, class and phylum) [START_REF] Clarke | Non-parametric multivariate analyses of changes in community structure[END_REF]Warwick, 1998, 2001a) and the indices were calculated using PRIMER 5

(Software package from Plymouth Marine Laboratory, UK) (Table 1). Higher values of these indices suggest that, on average, the species in the assemblage are not closely related, being higher biodiversity.

3. The Berger-Parker, Margalef and Shannon-Wiener indices were also calculated (see Table 1 for details). The Berger-Parker takes values from 0 to 1, with higher values corresponding to lower diversity. An increase on Margalef index is generally accepted as an indication of quality improvement. The log 2 was used to compute the Shannon-Wiener; therefore the results were expressed in bits/individual and higher values correspond to higher diversity.

4. The Infaunal Trophic Index, ITI, was calculated following [START_REF] Codling | Development of a biotic index for the assessment of pollution status of marine benthic communities[END_REF].

For each sample, the taxa contributing to the top 80% of the abundance have been allocated to trophic groups. In addition, in samples with ten or less identified taxa, all were taken into consideration for calculation purposes. Once generated, the list of the unclassified taxa was compared with the lists provided by [START_REF] Word | The Infaunal Trophic Index[END_REF][START_REF] Word | The Infaunal Trophic Index, a functional approach to benthic community analyses[END_REF]. Table 2 shows the Mondego estuary species' trophic classification. Values of the index vary from 0 to 100 with low values indicating degraded conditions. ITI values were used to classify areas of seabed into either "normal" (ITI values 100-60), "modified" (60-30) or "degraded" (30-0) [START_REF] Bascom | Establishing boundaries between normal, changed and degraded areas[END_REF].

5. Exergy expresses the work capacity relative to the same system but at thermodynamic equilibrium with its environment, which means that all the components would be inorganic at the highest possible oxidation state and homogeneously distributed at random in the system (no gradients). Of course exergy cannot be measured precisely, as the total internal energy content of a body cannot be measured exactly. Nevertheless, with regard to far from equilibrium systems (e.g. ecosystems), a reasonable approximation can be obtained if exergy is calculated only from the chemical potentials (which are in fact extremely dominant in ecosystems) [START_REF] Jørgensen | Principles of pollution abatement[END_REF]. This approximate measure can be called Eco-Exergy index [START_REF] Marques | Analysis of the properties of exergy and biodiversity along an estuarine gradient of eutrophication[END_REF]. If the total biomass in the system remains constant through time, then the variation of Eco-Exergy will be a function of only the quality of the biomass (structural complexity) or, in other words, of the information embedded in the biomass, which may be called Specific Eco-Exergy. Both indices were calculated using macrofauna biomass data and -values given by [START_REF] Jørgensen | Calculations of exergy for organisms[END_REF] according to the algorithm shown in Table 1.

Criteria used to assess the ecological indicators' success/failure

According to [START_REF] Borja | Assessing the environmental quality status in estuarine and coastal systems: comparing methodologies and índices[END_REF] indicate low ecological quality status (Table 1).

2) Determining the level of agreement provided by an index with best professional judgment. In our case study, the indices scores were compare with the existing knowledge about the ecological changes of the system (e.g. Marques et al., 1993[START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF][START_REF] Pardal | Impact of eutrophication on the life cycle, population dynamics and production of Ampithoe valida (Amphipoda) along an estuarine spatial gradient (Mondego estuary, Portugal)[END_REF][START_REF] Pardal | Assessing environmental quality: a novel approach[END_REF]Dolbeth et al., 2007;[START_REF] Cardoso | Dynamic changes in seagrass assemblages under eutrophication and implications for recovery[END_REF][START_REF] Cardoso | The impact of extreme flooding events and anthropogenic stressors on the macrobenthic communities' Dynamics[END_REF][START_REF] Lillebø | Management of a shallow temperate estuary to control eutrophication: the effect of hydrodynamics on the system nutrient loading[END_REF][START_REF] Lillebø | Apllying quality status criteria to a temperate estuary before and after mitigation measures to reduce eutrophication symptoms[END_REF]; and the present study) according to ultimate causal stressors that had induced such changes (e.g. hydrology and eutrophication) and to other potential stressors (proximate causes) generated by changes in hydrology and eutrophication (e.g. changes in seagrass cover and biomass (Z. noltii site); changes in the algal mats biomass (bare sediment site); nutrients levels; sediment grain size and total organic carbon (TOC). The "before the interruption-1985-1986" period was left out from this analysis due to it small number of replicates. Data were checked for normality using the Kolmogorov-Smirnov test and homogeneity of variances was assessed by the Levene's test. Data were transformed when homoscedasticity assumption was not met. The analyses were performed using the software package Minitab version 12.2. In each site, we expected no differences during the years of interruption, which strongly increased eutrophication effects and sedimentation, thus affecting Z. noltii and the post-management years.

Statistical analysis 2.4.1. Univariate analysis

Multivariate analysis

In order to examine the benthic macroinvertebrate communities' temporal dynamics in Mondego estuary South arm, a multivariate analysis was performed on biological data from the Z. noltii site and bare sediment site using the PRIMER 5.2.6 © software (Software package from Plymouth Marine Laboratory, UK). Data of abundance were square root transformed prior to analysis in order to scale down the effects on the ordination of highly abundant species [START_REF] Clarke | Change in marine communities: an approach to statistical analysis and interpretation[END_REF]. The Bray Curtis similarity matrix was calculated and used to generate a 2-dimensional plot with the non-metric multidimensional scaling (nMDS) technique [START_REF] Clarke | Non-parametric multivariate analyses of changes in community structure[END_REF]. The one-way ANOSIM test was performed to test for differences among time periods. Species contributing 70% to differences between periods were identified using similarity percentages (SIMPER; [START_REF] Clarke | Non-parametric multivariate analyses of changes in community structure[END_REF].

Results

Changes in macrofauna assemblages

In both sites, there was a significant decline in total macrofauna abundance in the "after the experimental opening" period (1-way ANOVA, Z. noltii site: F 1,58 = 12.59, p < 0.01; bare sediment site: F 1,56 = 24.06, p < 0.01) (in seagrass bed, from 1993 to 1995:

35181 -249966 ind m -2 ; from 2000 to 2002: 5785 -120375 ind m -2 ; in the bare sediment site, from 1993 to 1995: 6515 -167554 ind m -2 ; from 2000 to 2002: 3430 -28324 ind m -2 ) (Fig. 3A). Whereas the total biomass, in the Z. noltii site, was significantly higher "during the interruption" period (1-way ANOVA, F 1,58 = 7.63, p < 0.01), ranging from 58 to 158 g m -

2

. In the bare sediment site, the total biomass was significantly higher after the implementation of preliminary management measures (range 20-75 g m -2 ) (1-way ANOVA, F 1,56 = 20.59, p < 0.01) (Fig. 3B). In both sites, more species were present "during the interruption" period compared to "after the experimental opening" period. (Fig. 3C) (1-way ANOVA, Z.noltii site: F 1,58 = 77.23, p < 0.01; bare sediment site: F 1,56 = 15.65, p < 0.01).

Generally, the Z. noltii site supported more species than the bare sediment site (in the Z.noltii site, from 1993 to 1995: 14-40 species; from 2000 to 2002: 9-25 species; in the bare sediment site, from 1993 to 1995: 7-28 species; from 2000 to 2002: 7-20 species).

In both sampling sites, the results of the one-way ANOSIM test showed significant differences (p < 0.001) between the two periods. Nonetheless, the Z. noltii site had higher Global R value than the bare sediment site (Global R = 0.42 and 0.32, respectively). The nMDS showed that the "after the experimental opening" period is closer to the reference points, particularly in the Z. noltii site (Fig. 4).

The SIMPER analysis showed that in the Z. noltii site, Hediste diversicolor, Hydrobia ulvae, Scrobicularia plana, Cyathura carinata, Tubificidae, Alkmaria romijni, Heterosmastus filiformis, Capitela capitata, Cerastoderma edule and Chaetozone setosa were the species responsible for the separation between periods (Table 3). During the period of works and closure of the channel, the abundance of some polychaete species (A. romijni, H. filiformis, C. capitata and C. setose ), of Tubificidae and H. ulvae dramatically increased as compared to the period before disturbance (1986) and also to the period following mitigation measures (1998-2002) (Table 3).

Ecological indicators

In both sites, Margalef (Fig. 5A) exhibited significant differences between the two time periods considered (1-way ANOVA, Z. noltii site: F 1,58 = 58.70, p < 0.01; bare sediment site: F 1,56 = 12.29, p < 0.01), presenting unexpected higher values "during the interruption" period. Overall, Margalef was always higher in the Z. noltii site than the bare sediment site. The qualitative measure based on the presence/absence of species, TTD (Fig. 5B), showed similar behaviour (1-way ANOVA, Z. noltii site: F 1,58 = 76.82, p < 0.01; bare sediment site: F 1,56 = 15.12, p < 0.01).

In contrast to the diversity indices of Margalef and TTD, there were no difference in the Taxonomic Distinctness between periods (Fig. 5C), (p > 0.05), and the index was not larger in the Z. noltii compared to the bare sediment site.

The Taxonomic Diversity (Fig. 5D) and ITI (Fig. 5E) were unable to detect temporal changes in the Z. noltii site, but, as theoretically expected, exhibited lower values "during the interruption" period at the bare sediment site (1-way ANOVA, F 1,56 = 10.67, p < 0.01 and F 1,56 = 6.49, p < 0.01, respectively). Both indices were lower at Z. noltii site, differing from the behaviour that we could predict based on the authors classification and our own knowledge regarding the system. According to [START_REF] Bascom | Establishing boundaries between normal, changed and degraded areas[END_REF], ITI classified the benthic community in the bare sediment site community as a "slightly modified" to "normal".

The Shannon-Wiener (Fig. 5F) values were higher in the bare sediment site, because of the dominance of H. ulvae in the Z. noltii site (almost 90% of the total abundance).

While no significant differences were detected between the pre-and post-management periods in the reference site (1-way ANOVA, F 1,58 = 0.14, p = 0.71), in the bare sediment site, the values were significantly higher "after the experimental opening" period (1-way ANOVA, F 1,56 = 9.68, p < 0.01). As expected, the Berger-Parker dominance (Fig. 5G)

showed the opposite behaviour, with higher values in the Z. noltii site, again due to H.

ulvae's dominance. In the bare sediment site, the dominance significantly decreased after the implementation of the management plan (1-way ANOVA, F 1,56 = 12.82, p < 0.01).

A significantly higher AMBI value (Fig. 5H) was found "during the interruption" period at the reference site (1-way ANOVA, F 1,58 = 5.91, p < 0.02), indicating worse ecological conditions during this time period. In the bare sediment site, no differences regarding the pre-or post restoration periods were detected (1-way ANOVA, F 1,56 = 2.53, p > 0.05). This index showed similar values in both sites.

In the Z. noltii site, Eco-Exergy and Specific Eco-Exergy (Fig. 5I andJ) were significantly higher "during the interruption" period (1-way ANOVA, Eco-Exergy: F 1,58 = 8.66, p < 0.01; Specific Eco-Exergy : F 1,58 = 56.85, p < 0.01). In the bare sediment site, these indices showed the opposite behaviour, with significantly higher values "after the experimental opening" period (1-way ANOVA, Eco-Exergy: F 1,56 = 25.24, p < 0.01; Specific Eco-Exergy : F 1,56 = 4.8, p < 0.03), indicating higher biomass structure and information embedded in it. The indices had higher values in the Z. noltii site than in the bare sediment site, following the behaviour that we would predict based on a priori established criteria.

Discussion

Ecological indicators performance

Our analyses on the long-term changes in the composition of macrofauna assemblages confirmed the previous studies (e.g. Dolbeth et al., 2007;[START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF][START_REF] Cardoso | The impact of extreme flooding events and anthropogenic stressors on the macrobenthic communities' Dynamics[END_REF] and highlighted that the system seems to be recovering from the environmental changes occurred during the 1990s. At both sites, during the 1993-1995, the total dissolved inorganic nitrogen in the water column was higher than "after the experimental opening period" [START_REF] Leston | The response of primary producer assemblages to mitigation measures to reduce eutrophication in a temperate estuary[END_REF], the macroalgae assemblages were dominated by Ulvaceans, the Z. noltii cover and biomass reached their lowest values [START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF] and the abundance of the opportunistic species increased. After 1997/1998, the macrofauna and macrophyte assemblages showed a partial recovery towards the reference situation.

The selected ecological indicators are based upon community characteristics that are typically part of, or inherent to, the diversity of definitions of biotic integrity that include elements of species diversity, abundance, energy-flow-food web structure, maintenance of complexity and self-organization [START_REF] Borja | Assessing the environmental quality status in estuarine and coastal systems: comparing methodologies and índices[END_REF]. Therefore, according to the classification proposed, we would expect significantly higher values for Margalef, Shannon-Wiener, Taxonomic Distinctness measures, ITI, Eco-Exergy and Specific Eco-Exergy; AMBI and Berger-Parker should present significant lower values during the posmanagement period, at both sites. In the Z. noltii site, only AMBI reflected the ecological expectation, of a community with higher abundance of opportunist species "during the interruption" period. Margalef, TTD and the TBI indicated higher diversity and biomass structure/information "during the interruption" period, indicating that this community was apparently able to obtain some benefit from the existing conditions (higher nutrient input, lower water residence time). In the bare sediment site, with the exception of Margalef and TTD, all ecological indicators able to detect significant differences between time periods, performed as initially projected, indicating a better ecological status "after the experimental opening" period.

Furthermore, even if hydrology and eutrophication (both water column related), were probably the major drivers "shaping" South arm environmental conditions, the two study sites were differently affected by distinct stressors resulting from these drivers. The bare sediment site, characterised by sandy sediments, which have not supported rooted macrophytes for more than 15 years, lower organic matter content (mean 3.0%±1.14) and lower water-flows (0.8-1.2 m s -1 ) has been regularly covered by extensive algal mats (e.g. [START_REF] Lillebø | Management of a shallow temperate estuary to control eutrophication: the effect of hydrodynamics on the system nutrient loading[END_REF]Dolbeth et al., 2007;[START_REF] Marques | Analysis of the properties of exergy and biodiversity along an estuarine gradient of eutrophication[END_REF][START_REF] Marques | Impact of eutrophication and river management within a framework of ecosystem theories[END_REF][START_REF] Marques | Monitoring the Mondego estuary. Anthropogenic changes and their impact on ecological quality[END_REF][START_REF] Cardoso | The impact of extreme flooding events and anthropogenic stressors on the macrobenthic communities' Dynamics[END_REF]. These macroalgal blooms have well-documented effects, causing dramatic changes in sediment chemistry, nutrient release, rates of oxygen consumption and accumulation of metal sulfides and H 2 S in sediment pore waters. These shifts are profound changes in habitat quality for invertebrate assemblages (e.g. [START_REF] Pearson | Macrobenthic succession in relation to organic enrichment and pollution of the marine environment[END_REF][START_REF] Raffaelli | Ecological impact of green macroalgal blooms[END_REF][START_REF] Cloern | Our evolving conceptual model of the coastal eutrophication problem[END_REF]Patrício & Marques, 2006). On the other hand, the Z. noltii site, represents a distinct complex and heterogeneous habitat with particular physical and chemical conditions, characterised by muddy sediments covered with rooted macrophytes, higher organic matter content (mean 6.2%±1.76), higher waterflow velocity (1.2-1.4 m s -1 ) (e.g. Dolbeth et al., 2007) and supporting a higher diversity of species (Cardoso, 2004). Knowing that the direct and indirect responses to changes in hydrology and nutrient enrichment are expressed in varying degrees in different habitats [START_REF] Cloern | Our evolving conceptual model of the coastal eutrophication problem[END_REF], and taking into consideration the effect of site-specific operative stressors, we have considered Z. noltii site as representing a better ecological condition. Therefore, we would expect higher values for Margalef, Shannon-Wiener, Taxonomic Distinctness measures, ITI, Eco-Exergy and Specific Eco-Exergy in the reference site; AMBI and Berger-Parker would present the opposite behaviour. In fact, only Margalef, TTD and TBI followed the behaviour that we could predict based on our knowledge of the system (Fig. 4). Taxonomic Distinctness and AMBI presented similar values at both sites, indicating that these two metrics may be of limited efficacy in understanding the environmental conditions being investigated. AMBI values were dominated by H. ulvae abundance. As numerous estuarine species, this gastropod belongs to ecological group III, which comprises species tolerant to an excess of organic matter enrichment, a widespread stressor in transitional systems. Such explanation is supported by considering that Shannon-Wiener and Berger-Parker indices indicated that the Z. noltii site was of lower quality than the bare sediment.

These indices are based on measures of equitability and dominance, and can vary unexpectedly if one species becomes dominant, as is the case of H. ulvae. Indeed, excluding H. ulvae from the dataset resulted in opposite results. H. ulvae is a very mobile species that may migrate during low tide towards areas more protected and with high food availability, as is the case of Z. noltii site. Therefore, more attention should be paid to the ecology and distribution of species to decrease the probability of giving an erroneous classification based on a species with high mobility. Furthermore, it should be noticed that in estuaries species dominance occurs more often than on coastal systems, even in relatively pristine conditions, due the large natural variability of environmental parameters.

According to [START_REF] Constanza | Ecosystem health: new goals for environmental management[END_REF] and Elliot and Quintino (2007), in these ecosystems, the "positive effect" of natural environmental variability is shown as the ability of some organisms, which can tolerate the adverse and variable environmental conditions, to capitalise on the lack of inter-specific competition and thus achieve high population densities. Therefore, caution should be taken when interpreting equitability and dominance results based in estuarine macrobenthic assemblages.

Following [START_REF] Word | The Infaunal Trophic Index, a functional approach to benthic community analyses[END_REF], ITI showed a higher value in the bare sediment site, where the sediment organic matter content was lower. However, in the Z. noltii site the sediment organic matter has a natural source (macrophyte primary production). For that reason, the metric is not correctly capturing the environmental conditions being investigated.

Moreover, ITI is difficult to apply in environmental studies once determining the organisms' diet may become reasonably problematic. Generally, the real diet (studied observing stomach content) is difficult to establish, and can vary from one population to another among the same taxonomic entity (e.g. H. filiformis was classified by [START_REF] Brown | The effects of crude oil pollution on marine organisms: a literature review in the South African context: conclusions and recommendations[END_REF] as subsurface deposit feeder and as surface detritus feeder by [START_REF] Word | The Infaunal Trophic Index, a functional approach to benthic community analyses[END_REF]).

According to our study, it seems that the ecological indicators response was not always temporally and spatially consistent (Table 4). These results are quite interesting, raising further questions.

It is well known, that estuaries are particularly challenging due to strong spatial, seasonal and inter-annual variations of environmental characteristics that influence benthic communities [START_REF] Chainho | Influence of seasonal variability in benthic invertebrate community structure on the use of biotic índices to assess the ecological status of a Portuguese estuary[END_REF]. Apart from salinity regime, plant presence/absence, sediment grain size and organic matter content, the hydrological conditions, such as flood and drought periods, may have great impact in communities and thus on diversity measures [START_REF] Pearson | Macrobenthic succession in relation to organic enrichment and pollution of the marine environment[END_REF][START_REF] Gray | Marine biodiversity: patterns, threats and conservation needs[END_REF][START_REF] Raffaelli | Ecological impact of green macroalgal blooms[END_REF][START_REF] Cloern | Our evolving conceptual model of the coastal eutrophication problem[END_REF][START_REF] Salas | Ecological indicators in coastal and estuarine environmental quality assessment. A user-friendly guide for practitioners[END_REF][START_REF] Chainho | Influence of seasonal variability in benthic invertebrate community structure on the use of biotic índices to assess the ecological status of a Portuguese estuary[END_REF]Dolbeth et al., 2007;[START_REF] Cardoso | The impact of extreme flooding events and anthropogenic stressors on the macrobenthic communities' Dynamics[END_REF][START_REF] Teixeira | Ecological indices tracking distinct impacts along disturbance-recovery gradients in a temperate NE Atlantic Estuary -guidance on reference values[END_REF].

In summary, we may say that although hydrology and eutrophication were the major drivers, other ultimate causal factors, probably also played a role during the study period, making it difficult to unambiguously distinguish the anthropogenic disturbance effects from the natural temporal wide-scale variability. In any case, the ecological indicators' results

were constrained to the type and magnitude of the drivers' combined effect that occurred at either or both sites during the study period. Therefore, caution must be taken interpreting values isolated from the ecological community context.

Scientific implications and limitations of the study

This paper laid out some of the history of anthropogenic modifications and sketched a brief image of the ecological changes over the considered timeframe (derived partially from previously reported literature). In addition, many readers have pondered comparisons among all of the new indices that seem to continuously emerge (e.g. Díaz et al., 2004;[START_REF] Labrune | Characterization of the ecological quality of the coastal Gulf of Lions (NW Mediterranean). A comparative approach based on three biotic indices[END_REF][START_REF] Aubry | The use of environmental integrative indicators to assess seabed disturbance in estuaries and coasts: application to the Humber Estuary, UK[END_REF][START_REF] Quintino | The derivation, performance and role of univariate and multivariate indicators of benthic change: case studies at different spatial scales[END_REF][START_REF] Dauvin | The ecological quality status of the Bay of Seine and the Seine estuary: use of biotic indices[END_REF]Blanchet et al., 2008). The present study compared several ecological indicators (all using benthic macrofauna) to the image of change that can be cast through a set of measurements on alterations in water quality, hydraulic/sediments, benthic invertebrate and plant/algal communities, and spatial distribution. Even if this broad comparison has been made in the past, given the European initiatives to use indices to monitor estuarine system health, we hope this study could be a useful and directly applicable demonstration of the care which must be taken in selecting appropriate indices. Indeed, it is very challenging to see such different predictions from the diverse indices, given that the same data were used (albeit in different ways, but that was the aim). Notwithstanding, it is obvious for us, that the perfect way to compare indices would be to (ideally) set up an ecological situation, in as controlled a manner as possible, in which one knew the "precise, accurate and totally clear answer" already and where strong putative gradients were selected for validation testing (e.g. [START_REF] Borja | Marine Biotic Index to establish the ecological quality of softbottom benthos within European estuarine and coastal environments[END_REF][START_REF] Quintino | The derivation, performance and role of univariate and multivariate indicators of benthic change: case studies at different spatial scales[END_REF]. Nevertheless, the available ecological tools will have to work in "existent" ecosystems and available datasets. In our case study the ecological indicators diagnostic ability was only tested in two sites: a "Z. noltii site" and a "bare sediment site" during two different time periods, characterised by different levels of impact. We did not cover the five categories featured within the WFD (bad, poor, moderate, good and high). Therefore, this study only addresses two conditions and for that reason provides a limited assessment.

Clearly, with other case-studies, further steps will have to be done to fulfil the WFD requirements. The second phase of the European WFD Intercalibration exercise will certainly be critical, particularly for transitional systems. Even so, it is our conviction that initially the proposed ecological indicators should be able to give the correct assessment in this case, particularly, when some of these metrics are being selected and combined into multimetric indices used to assess the ecological quality status all over Europe with obvious legal implications. We truly agree that without a sound and obvious ecological foundation, an index will not be policy-relevant and therefore difficult to use to make policy choices [START_REF] Borja | Assessing the environmental quality status in estuarine and coastal systems: comparing methodologies and índices[END_REF], thus the results must be justifiable and scientifically defensible. Experiences gather from different ecosystems, distinct geographical ranges, various hydrological regimes, and diverse pressure types are vital and urgent.

Conclusions

Relevant conclusions can be taken from our study: (1) Taxonomic Distinctness did not detect temporal and spatial changes; (2) ITI did not detect spatial changes largely because estuarine macrofauna shows diet plasticity and is able to adapt its trophic role;

(3) The equitability and dominance measures were totally influenced by the presence of a common species with high mobility, (4) The best distinction between sites was given by Margalef,TTD and TBI,(5) The indicators ability to capture temporal changes was different in each site, revealing habitat-dependency and (6) although some ecological indicators have been effective in differentiating the two time periods, not all have behaved consistently with ecological expectation. These considerations are crucial for management, once an indicator is expected to respond to a particular pressure and lead to a human response to reduce environmental impacts through the systems' driving forces, with all the social and economic repercussions of that action. I. Species very sensitive to organic enrichment and present under unpolluted conditions; II. Species indifferent to enrichment, always in low densities with nonsignificant variations with time; III. Species tolerant to excess of organic matter enrichment; IV. Second-order opportunist species, mainly small sized polychaete; V. First-order opportunist species, essentially deposit-feeders Normal 0-1.2 Slightly polluted 1.2-3.2 Moderately polluted 3.2-5.0 Highly polluted 5.0-6.0 Azoic 6.0-7.0

Diversity

Average Taxonomic Diversity [START_REF] Warwick | New 'biodiversity' measures reveal a decrease in taxonomic distinctness with increasing stress[END_REF])

[ ] ( ) [ ] 2 / 1 / - ∑ ∑ < = n n j i j x i x ij ω ∆
Average Taxonomic Distinctness [START_REF] Warwick | New 'biodiversity' measures reveal a decrease in taxonomic distinctness with increasing stress[END_REF] [

] [ ]

∑ ∑ < ∑ ∑ < = ∆ j x i x j i j i j x i x ij δ ω / *
Total Taxonomic Distinctness [START_REF] Warwick | New 'biodiversity' measures reveal a decrease in taxonomic distinctness with increasing stress[END_REF]

( ) ( ) [ ] [ ] ∑ - ∑ ≠ = + ∆ 1 s j i ij i s ω i x : abundance of the i th of s species observed; ( ) ∑ = i i
x n : total number of individuals in the sample; ij ω : "distinctness weight" given to the path length linking species i and j in the taxonomy Higher value, higher biodiversity -on average species in the assemblage are not closely (phylogenetic) related 70% to the difference between periods in the Z. noltii and in the bare sediment sites. for the three time periods: before the perturbation under assessment (1985 and 1986), period without connection between arms (1993, 1994 and 1995) and after experimental mitigation measures (1999)(2000)(2001)(2002).

  , index validation should ideally include (1) testing of the index using an independent data set, different from the index development data set; (2) setting a priori correct classification criteria and/or (3) presentation of a strong a posteriori justification for use based upon best professional judgment. The first requirement was fulfilled using our own data set from two intertidal sites in the Mondego estuary. In order to assess the indices success or failure to discriminate the ecological quality status, two main criteria were used: 1) The classification proposed by the index authors. High values of Margalef, Shannon-Wiener, Taxonomic distinctness indices, ITI and the thermodynamically based indices (TBI) indicate high ecological quality status. High values of Berger-Parker and AMBI

  One-way ANOVA, with time as fixed factor (2 levels: "during the interruption of flow between the North and South arms -1993 to 1995" and "after the experimental opening of the connection between the two arms of the estuary-2000 to 2002") and 30 and 29 replicates for Z. noltii site and bare sediment site, respectively, was used to test for temporal differences in each of the two sites (with respect to total abundance, total biomass, number of species, Margalef, Shannon-Wiener, Berger-Parker, Biotic Coefficient (AMBI), Taxonomic Distinctness indices, ITI, Eco-Exergy and Specific Eco-Exergy).

  Blanchet et al. (2008), testing five biotic indices (AMBI, BENTIX, Shannon-Wiener diversity, BQI and BOPA) in two semi-enclosed, sheltered coastal ecosystems and in one transitional water body situated along the Western French coast, also found that these five indices rarely agreed with each other, describing very different pictures of the overall ecological quality status of the three study areas. A similar conclusion was made by[START_REF] Chainho | Influence of seasonal variability in benthic invertebrate community structure on the use of biotic índices to assess the ecological status of a Portuguese estuary[END_REF] studying the Mondego estuary subtidal assemblages, using four indices(Margalef, Shannon- Wiener, AMBI and W-Statistic).

Figures captions Figure 1 .

 captions1 Figures captions

Figure 2 .

 2 Figure 2. Temporal changes of Z. noltii meadows' area (ha) and biomass (g AFWD m -2), green macroalgae biomass (g AFWD m -2 ) and the connection width between the two estuarine arms. Black areas correspond to Z. noltii site, bare sediment site and connection between the two arms. NA-north arm, SA-south arm.

Figure 3 .

 3 Figure 3. Macrofauna A. total abundance, B. total biomass and C. number of species in the Z. noltii site and in the bare sediment site, during the study period.

Figure 4 .

 4 Figure 4. A Two-dimensional, non-metric, multidimensional scaling (nMDS) ordination based on abundance data. A. Z. noltii site and B. bare sediment site,

  Each point represents a different time period at the same site. Stress values are shown for each nMDS plot indicating the similarity between the multi-dimensional representation and the 2-dimension plot.

Figure 5 .

 5 Figure 5. Temporal and spatial variation of different ecological indicators applied to data of the intertidal communities of the South arm of the Mondego estuary, during the study period.

  Figure

  Figure

FigureI.

  Figure

Table 1 .

 1 Díaz, R.J., Solan, M., Valente, R.M., 2004. A review of approaches for classifying benthic habitats and evaluating habitat quality. Journal of Environmental Management 73, 165-181. Dixon, I.M.T., Moore, P.G., 1997. A comparative study of the tubes and feeding behaviour of eight species of corophioid Amphipoda and their bearing on phylogenetic relationships within the Corophioidea. Philosophical transactions of the Royal Society of London 352, 93-112. Ecological indicators: algorithms and ecological classification.

	Dolbeth, M., Cardoso, P.G., Ferreira, S.M., Verdelhos, T., Raffaelli, D., Pardal, M.A., 2007.
	Anthropogenic and natural disturbance effects on a macrobenthic estuarine community over
	a 10-year period. Marine Pollution Bulletin 54, 576-585.
	EC, 2000. Directive 2000/60/EC of the European Parliament and of the Council establishing a
	framework for community action within the field of water policy PE-CONS 3639/1/00, 72 p.
	! " ! #		$ %	&	&
	' (	)	# & *	$ %	+ . Official Journal of the European Union, 22p.
	Elliot, M., Quintino, V., 2007. The Estuarine Quality Paradox, Environmental Homeostasis and the
	difficulty of detecting anthropogenic stress in natural stressed áreas. Marine Pollution Bulletin

Table 4 .

 4 Summary of the ecological indicators trends obtained for the Z. noltii and bare sediment sites, considering "during the interruption" and "after the experimental opening" (after) periods. In bold are the values that are compliant with a priori established classification criteria.

	728
		Before the interruption	During the interruption	After the experimental opening
	Zostera noltii site
	Hediste diversicolor	1189	19	523
	Tubificidae	0	2944	714
	Alkmaria romijni	204	993	443
	Heterosmastus filiformis	265	550	357
	Capitela capitata	114	241	92
	Chaetozone setosa	0	1756	31
	Hydrobia ulvae	1960	80691	41862
	Scrobicularia plana	596	1023	724
	Melita palmata	83	336	75
	Cyathura carinata	271	217	173
	Oligochaete	0	220	100
	Streblospio shrubsolii	60	349	138
	Cerastoderma edule	13	519	74
	Bare sediment site
	Hediste diversicolor	523	50	265
	Hydrobia ulvae	1166	27192	5573
	Alkmaria romijni	93	6076	1483
	Cyathura carinata	638	2043	2177
	Scrobicularia plana	319	2347	2763
	Amphicteis gunneri	106	0	0

Acknowledgments

The present study was carried out in the scope of WADI (INCO-CT-2005-015226), EXTREMIS (III/36/2008) and RECONNECT (PTDC/MAR/64627/2006) research projects. It was supported by the FCT (Portuguese National Board of Scientific Research) through three grants (SFRH/BPD/26604/2006, SFRH/BPD/41701/2007, SFRH/BD/24430/2005).

This work had the institutional support of INAG-National Water Institute (Portugal). Authors wish to thank two anonymous referees for their critical review of the first version of the manuscript and for their useful comments. We would also like to sincerely thank F. Rossi for helping us to improve the quality of the manuscript. Special thanks to all who assisted us, at any moment, during field and laboratory work.

Berger-Parker Index [START_REF] Berger | Diversity of planktibuc Foraminifera in deep sea sediments[END_REF] N n D max = nmax: number of individuals of the most abundant species; N: total number of individuals

Low diversity 1 High diversity 0

Margalef Index [START_REF] Margalef | Perspectives in Ecological Theory[END_REF])

S :number of species found; N : total number of individuals

Low diversity

Low value

Shannon-Wiener Index

Ecological strategies

Infaunal Trophic Index [START_REF] Word | The Infaunal Trophic Index[END_REF])


Didier Fraix-Burnet
email: fraix@obs.ujf-grenoble.fr

ASTROCLADISTICS: MULTIVARIATE EVOLUTIONARY ANALYSIS IN ASTROPHYSICS

Keywords: astrophysics, cladistics, multivariate, evolution, galaxies, continuous characters

The Hubble tuning fork diagram, based on morphology and established in the 1930s, has always been the preferred scheme for classification of galaxies. However, the current large amount of data up to higher and higher redshifts asks for more sophisticated statistical approaches like multivariate analyses. Clustering analyses are still very confidential, and do not take into account the unavoidable characteristics in our Universe: evolution. Assuming branching evolution of galaxies as a 'transmission with modification', we have shown that the concepts and tools of phylogenetic systematics (cladistics) can be heuristically transposed to the case of galaxies. This approach that we call "astrocladistics", has now successfully been applied on several samples of galaxies and globular clusters. Maximum parsimony and distance-based approaches are the most popular methods to produce phylogenetic trees and, like most other studies, we had to discretize our variables. However, since astrophysical data are intrinsically continuous, we are contributing to the growing need for applying phylogenetic methods to continuous characters. classification;

INTRODUCTION

The extragalactic nature of galaxies was discovered by Edwin Hubble less than 90 years ago. He also discovered the expansion of the Universe and established in 1936 the famous tuning fork diagram or Hubble diagram based on the idea that elliptical galaxies should evolve into flattened systems such as ordinary or barred spiral galaxies. Since then, telescopes and the associated detectors have made remarkable progress. We are now able to study galaxies in great detail, identifying individual stars, gas and dust clouds, as well as different stellar populations. Imagery brings very fine structural details, and spectroscopy provides the kinematical, physical and chemical conditions of the observed entities at different locations within the galaxy. For more distant objects, information is scarcer, but deep systematic sky surveys gather millions of spectra for millions of galaxies at various redshifts.

Like paleontologists, we observe objects from the distant past (galaxies at high redshift), and like evolutionary biologists, we want to understand their relationships with nearby galaxies, like our own Milky Way. Strangely enough, the Hubble classification, depicted in the Hubble diagram, is still frequently used as a support to describe galaxy evolution, even though it ignores all observables except morphology (Fraix-Burnet et al 2006a, Hernandez and Cervantes-Sodi 2006). Can the enormous amount of very detailed observations and galaxy diversity be solely depicted by few large families characterized only by their global shape? There obviously must be a better way to exploit the data gathered by very large telescopes and their sophisticated detectors, and to account for the complicated physical processes that lead to the diversification of galaxy.

Abandoning the one-parameter classification approach and using all available descriptors means taking a methodological step equivalent to the one biologists took after Adanson and Jussieu in the 18th century. One basic tool, the Principal Component Analysis, is relatively well-known in astrophysics (e.g. [START_REF] Cabanac | Classification and redshift estimation by principal component analysis[END_REF][START_REF] Recio-Blanco | Multivariate analysis of globular cluster horizontal branch morphology: searching for the second parameter[END_REF]. However, only a very few attempts to apply multivariate clustering methods have been made very recently (Chattopadhyay and [START_REF] Chattopadhyay | Objective classification of spiral galaxies having extended rotation curves beyond the optical radius[END_REF]Chattopadhyay , 2007;;Chattopadhyay et al 2007[START_REF] Chattopadhyay | Globular Clusters in the Milky Way and Dwarf Galaxies -A Distribution-Free Statistical Comparison[END_REF][START_REF] Chattopadhyay | Horizontal Branch Morphology of Globular Clusters: A Multivariate Statistical Analysis[END_REF]. Sophisticated statistical tools are used in some areas of astrophysics and are developing steadily, but multivariate analysis and clustering techniques have not much penetrated the community.

A supplementary difficulty is that evolution, an unavoidable fact, is not correctly taken into account in most classification methods. By mixing together objects at different stages of evolution, most of the physical significance and usefulness of the classification is lost. In practice, the evolution of galaxies is often limited to the evolution of the properties of the entire population as a function of redshift [START_REF] Bell | Galaxy Assembly. Planets to Cosmology: Essential Science in Hubble's Final Years[END_REF]. Since environment (the expanding Universe) and galaxy properties are so much intricate, this kind of study is relevant to a first approximation. However, recent observations have revealed that galaxies of all kinds do not evolve perfectly in parallel, as illustrated for instance by the so-called downsizing effect which shows that large galaxies formed their stars earlier than small ones (e.g. [START_REF] Neistein | Natural downsizing in hierarchical galaxy formation[END_REF]). New observational instruments now brings multivariate information at different stages of evolution, and in various evolutive environments. In this multivariate context, we believe that the notion of "evolution", easy to understand for a single parameter, is advantageously replaced by "diversification".

In this conference, I present our efforts to implement a method to reconstruct the "galactogenesis", based on tools and concepts largely developed in evolutionary biology and bioinformatics. I first present cladistics and some of our results. Then I replace this approach among multivariate statistical methods and explain why we also need these other tools to explore the parameter space. Finally, I present one of our developments that makes the astrophysical problematics to join current bioinformatics and mathematics studies about the use of continuous characters in reconstructing phylogenies.

Multivariate clustering methods compare objects with a given measure and then gather them according to a proximity criterion. Distance analyses are based on the overall similarity derived from the values of the parameters describing the objects. The choice of the most adequate distance measure for the data under study is not unique and remains difficult to justify a priori. The way objects are subsequently grouped together is also not uniquely defined. Cladistics uses a specific measure that is based on characters (a trait, a descriptor, an observable, or a property, that can be given at least two states characterizing the evolutionary stages of the object for that character) and compares objects in their evolutionary relationships [START_REF] Wiley | The Compleat Cladist: A Primer of Phylogenetic Procedures[END_REF]. Here, the "distance" is an evolutionary cost. Groupings are then made on the basis of shared or inherited characteristics, and are most conveniently represented on an evolutionary tree.

Character-based methods like cladistics are better suited to the study of complex objects in evolution, even though the relative evolutionary costs of the different characters is not easy to assess. Distance-based methods are generally faster and often produce comparable results, but the overall similarity is not always adequate to compare evolving objects. In any case, one has to choose a multivariate method, and the results are generally somewhat different depending on this choice (e.g. [START_REF] Buchanan | Phenetics, cladistics, and the search for the Alaskan ancestors of the Paleoindians: a reassessment of relationships among the Clovis, Nenana, and Denali archaeological complexes[END_REF]). However, the main goal is to reveal a hidden structure in the data sample, and the relevance of the method is mainly provided by the interpretation and usefulness of the result.

Multivariate evolutionary classification in astrophysics has been pioneered by the author (Fraix-Burnet et al 2006a, 2006c, Fraix-Burnet 2009). Called astrocladistics, it is based on cladistics that is heavily developed in evolutionary biology. Astrocladistics has been first applied to galaxies (Fraix-Burnet et al 2006b) because they can be shown to follow a transmission with modification process when they are transformed through assembling, internal evolution, interaction, merger or stripping. For each transformation event, stars, gas and dust are transmitted to the new object with some modification of their properties. Cladistics has also been applied to globular clusters (Fraix-Burnet et al 2009), where interactions and mergers are probably rare. These are thus simpler stellar systems, even though we have firm evidence that internal evolution can create another generation of stars and that globular clusters can lose mass. Basically, the properties of a globular cluster strongly depend on the environment in which it formed (chemical composition and dynamics), and also on the internal evolution which includes at least the aging of its stellar populations. Since galaxies and globular clusters form in a very evolving environment (Universe, dark matter haloes, galaxy clusters, chemical and dynamical environment), the basic properties of different objects are related to each other by some evolutionary pattern.

In both cases, there is a kind of transmission with modification process, which justifies a priori the use of cladistics. It must be clear that this is not a "descent with modification" in the sense that there is no replication. But evolution does nevertheless create diversity. We are dealing with phylogeny (relationships between species), not with genealogy (relationships between individuals). Since a multivariate classification of galaxies is not yet available, we assume that each object represents a "species" that will have to be defined later on. The word "progenitor", often used for galaxy evolution, is to be understood A cladistic analysis works as follows. One first builds a matrix with values of the four parameters for all galaxies. In contrast to multivariate distance methods, undocumented values are not a problem in cladistic analyses. The values for each parameter are discretized into 30 bins representing supposedly evolutionary states. Discretization of continuous variables is quite a complex problem, especially in the evolutionary context (see Section 4). Here, we took equal-width bins. The choice of the number of bins cannot be made in a simple objective way. A priori, one could consider a compromise between an adequate sampling of continuous variables and the uncertainties on the measurements. The first constraint is given by the software (32 in this case). The second one would a priori give a lower limit of something like total range/uncertainty, but Shannon's theorem would multiply this by 2. Hence, 30 bins would account for about 7% measurement errors. Even so, border effects always imply that some objects could belong to a bin or its neighbor, a process that add some more artificial noise. The best way to avoid this effect is to make several analyses with different number of bins and check that the result does not depend on this number. From our experience in astrocladistics, we know that the results are very generally identical between 20 and 30, and can differ with 10 bins. For lower number of bins, the results are very dependent on this number.

Then, all possible arrangements of galaxies on a tree-like structure are con- structed, and using the discretized matrix, the total number of state changes is computed for each tree. The most parsimonious tree is finally selected. If several such trees are found, then a consensus (strict or majority rule) tree is built. The whole procedure is computerized since the number of arrangements is very large. The result (depicted on a cladogram) is a diversification scenario that should be confronted to other knowledge and parameters.

Figure 1 shows the cladogram obtained for globular clusters of our Galaxy (Fraix-Burnet et al 2009). Three groups are identified. The first one (in blue) has on average the lower ratio Fe/H that measures the proportion of heavy atomic elements that are processed within stars. This group is consequently considered as more primitive.

Figure 2 shows bivariate diagrams for the 4 parameters used in the analysis: logTe, that measures the temperature of stars that are at a specific point in their evolution, Fe/H, MV that is the total visible intensity (magnitude) and roughly indicates the mass of the globular cluster, and Age that can be measured quite precisely because all stars of a given globular cluster are formed nearly at the same time. However Age is not an intrinsic property discriminating evolutionary groups since it evolves in the same way for all. But we gave it a half weight to arrange the objects within each group. Looking at plots like Figure 2 and using other parameters (orbital elements, kinematics, more refined chemical abundances...), we are able to infer that each group formed during a particular stage of the assembly history of our Galaxy. The blue group is the older one. It formed during the dissipationless collapse of the protogalaxy. They are located mainly in the outer halo. The red group belongs to the inner halo and the corresponding clusters formed at a later stage during the dissipational phase of Galactic collapse, which continued in the halo after the formation of the thick disc and its globular clusters. These clusters were very massive before "star evaporation" took place. The latter group (green) formed during an intermediate and relatively short period and comprises clusters of the disk of our Galaxy.

STRUCTURE OF THE PARAMETER SPACE

The spectacular astrophysical results obtained so far were possible because the characters are pertinent from the evolutionary and physical points of view. However, they are too few to obtain a detailed classification of tens of thousands of galaxies. Since astrophysics is a science of observation, not of experimentation, it is difficult to obtain the intrinsic physical and chemical properties of galaxies or globular clusters especially if distant. The spectra carries all the information we gather, so that much data is available. But this information is generally not straightfully pertinent, either because it is not directly connected to physical parameters, or because there are redundancies or incompatibilities that are paticularly annoying in a cladistic analysis. Being as multivariate as possible ensures objectivity, but some care is required in selecting the parameters. Different multivariate tools are thus needed to explore the structure of the data.

The first and obvious approach is to use Principal Component Analysis (PCA) to reduce the dimensionality of the parameter space and avoid redundancy. It is possible to perform a cladistic analysis directly with the main components. However, the correlations revealed by such analyses are not necessarily intrinsic nor due to scaling effects. They can be generated by evolution. For instance, in Figure 3, the parameter logs, that measures the central velocity dispersion, is not physically linked to mgbfe, that measures to the global metallic composition of the stars, but they appear "correlated" through their respective evolutions. They possibly characterize two independent facets of galaxy diversification. Such information would be lost with the use of the PCA components.

Alternatively, it is possible to select the parameters that appear to have more weight in the PCA. In this approach, we both keep the evolutionary information and gain a precious knowledge on the parameters themselves, which are physical in contrast to the principal components. Parameters for subsequent analyses are thus selected according to their loadings and their physical meaning. At the end, we eliminate true redundancies (parameters that measure the same quantity) and undiscriminating parameters.

There are other techniques that surely could be interesting and that we have not yet explored. Rather, we are investigating approaches more specific to the problem of obtaining a phylogeny with continuous characters (see Section 4). It appears that the structure of the parameter space must be at least grossly understood before embarking on sophisticated clustering tools, especially if the Figure 3: Projection on the two first Principal Components of parameters available for a typical sample of galaxies. evolutionary information is to be retrieved. One difficulty in multivariate clustering of continuous and evolutive parameters, certainly not specific to astrophysics, is the cosmic variance that makes the groups to largely overlap in the parameter space.

Once the parameters are selected, it is instructive to compare results obtained with several multivariate clustering techniques in several different conditions (sub-samples, subsets of parameters). This provides additional insights on the parameters, and the different groupings can be confronted and analysed. Similarities reinforce the results, and mismatches can be interpreted in light of specificities of the techniques. In this way, we have obtained very convincing results on samples of hundreds of galaxies (Fraix-Burnet, Chattopadhyay, Chattopadhyay, Davoust and Thuillard, in prep). We demonstrate that multivariate clustering techniques and cladistics are the right tools to use. They yield results that are largely consistent, the cladogram providing the additional evolutionary relationshups between the groups. These groups cannot be guessed by any a priori splitting of the data, but we find that their properties happen to delineate different assembling histories, like for the globular clusters above.

CONTINUOUS CHARACTERS AND SPLIT NETWORKS

Cladistics is a difficult concept and its use in the case of continuous variables (morphometric data) is under intense investigation in the bioinformatics community (Gonzàles-Rosé et al 2008). In particular, the formal relationship between distance-based (clustering, continuous variables) and character-based (cladistics, discrete values) is a very interesting topic in itself for which astrophysics brings an entirely new kind of data [START_REF] Thuillard | Phylogenetic Applications of the Minimum Contradiction Approach on Continuous Characters[END_REF]. Maximum parsimony and distance-based approaches are the most popular methods to produce phylogenetic trees. While most studies in biology use discrete characters, there is a growing need for applying phylogenetic methods to continuous characters. Examples of continuous data include gene expressions [START_REF] Planet | Systematic analysis of DNA microarray data: ordering and interpreting patterns of gene expression[END_REF], gene frequencies [START_REF] Edwards | Reconstruction of evolutionary trees. Phenetic and Phylogenetic Classification[END_REF]1967), phenotypic characters [START_REF] Oakley | Independent contrasts succeed where ancestor reconstruction fails in a known bacteriophage phylogeny[END_REF] or some morphologic characters [START_REF] Macleod | Morphology, Shape and Phylogeny[END_REF]González-José et al. 2008). The simplest method to deal with continuous characters using maximal parsimony consists of discretizing the variables into a number of states small enough to be processed by the software.

Distance-based methods are applied to both discrete and continuous input data. Compared to character-based approaches, distance-based approaches are quite fast and furnish in many instances quite reasonable results. As pointed out by [START_REF] Felsenstein | Inferring phylogenies[END_REF], the amount of information that is lost when using a distance-based algorithm compared to a character-based approach is often surprisingly small. The use of continuous characters in distance-based methods may at first glance be less problematic than in character-based methods as algorithms like the Neighbour-Joining work identically on discrete or continuous characters. But also here it is often not easy to determine if the data can be described by a tree. In Thuillard and Fraix-Burnet (2009), we investigate the conditions for a set of continuous characters to describe a split network or a Xtree (split networks are generalized trees). We show that a set of m continuous characters can be described by a split network or a weighted X-tree if the mdimensional space representation of the taxa state values is on an orthogonal convex hull. In that case, there exists an order of the taxa, called perfect order, for which the distance matrix satisfies the Kalmanson inequalities for each character.

In practice, identifying a priori characters that comply to these conditions is difficult. For complex objects in evolution, this requires some good knowledge of the evolution of the characters together with some ideas about the correct phylogeny or at least a rough evolutionary classification. In astrophysics, the study of galaxy evolution has yet not reached this point. However, we can show that the approach presented in this work, based on the Minimum Contradiction Analysis [START_REF] Thuillard | Minimizing contradictions on circular order of phylogenic trees[END_REF][START_REF] Thuillard | Minimum contradiction matrices in whole genome phylogenies[END_REF]) is extremely valuable even in cases with very little a priori hints. Here, we illustrate how this method can be used in practice, in particular to discover structuring characters and objects (more can be found in Thuillard and Fraix-Burnet 2009). We have taken from Ogando et al (2008) a sample of 100 galaxies described by some observables and derived quantities.

In this first example, three variables have been chosen: logs, that is the logarithm of the central velocity dispersion, Kabs that is the absolute magnitude in the K band and measures the total quantity of near-infrared K light emitted mainly by the stars, and ldiam that is the diameter of the galaxy. They are all expected to be correlated with the total mass of the galaxy. In addition, Kabs and ldiam both strongly depend on the total number of stars. Figure 4 shows the three variables logs, Kabs, Idiam after ordering of the galaxies with the Minimum Contradiction Analysis that provides the best order (the closest to the perfect order). The two last characters fulfil quite well the conditions for perfect order and the associated distance matrices are well ordered. The first character is in very first approximation well ordered. The distance matrix is well approximated by an expression characterizing a line tree structure that do not show particular substructures.

As expected, the variables Kabs and ldiam are well correlated as can be seen from the ordered taxa corresponding to the best order. They are essentially responsible for the two groups seen on the distance matrix of Fig. 4 and thus for the line tree structure. The variable logs is only in first approximation correlated to the other variables. It is much noisier and has no obvious extremum.

In the second example, we have chosen two variables possibly very little correlated and with evolution function probably not far from the case yielding a convex hull. D/B is the ratio between the size of the disk to that of the bulge. This is a quantitative and objective measure of the morphology of galaxies, spiral (that are disky) systems having a high D/B and elliptical ones having low D/B. The OIII variable is an intensity measure of a region of the spectrum where the OIII emission or absorption line is present. A strong emission line indicates the presence of interstellar gas and thus probably star formation. It is a priori not directly physically linked to the morphology.

Despite some large deviations to perfect order for a number of galaxies, the level of contradiction is quite low as can be seen in Figure 5. The distance matrix can be described to a good approximation by a split network. A typical split network is shown in Figure 6. Obviously, split networks can represent more complex relationships than trees, but they are very difficult to use for inferring phylogenetic hypotheses, even more than reticulograms that take into account hybridization or horizontal transfer. In Thuillard and Fraix-Burnet (2009), we propose a simple a posteriori discretization of the variables that, in favorable cases, simplifies a split network into a X-tree. This is depicted in Figure 7 where the values of the variables are plotted in the best order of Figure 5.

Using a threshold value for the two variables OIII and D/B furnishing the discrete characters OIII large, OIII small, D/B large and D/B small. Compared to the peak in the D/B curve, the large peak for OIII is shifted to the right by roughly 5-10 galaxies. This shift generates some structure in the distance matrix (Figure 5). The galaxies can then be grouped into 3 different groups. The first group has large D/B values and large OIII except for 2 or 3 galaxies that show small OIII, the second group has small D/B and small OIII values while the group 3 contains all galaxies with small values of D/B and large OIII. This grouping can be represented on the tree of Figure 8.

The first group is made of disky galaxies, which have generally a high star formation rate. There are only 2 or 3 galaxies of this group that have a small OIII. The second and third groups have more bulgy galaxies. In the latter, some galaxies show an OIII line in emission as strong as in group 1. These galaxies seem to contradict the long known observation, that elliptical galaxies have less gas, and thus form less stars than spiral galaxies. This observation is believed to be due to historical reasons in the evolution of galaxies that are not at all clear. Is it due to a sweeping of the gas in the elliptical galaxies by the intergalactic medium? Or is it due to the way elliptical galaxies formed, Figure 8: The "phylogenetic" tree obtained from the discretization shown in Figure 7.

exhausting nearly all their gas in very efficient star formation episodes? Even if our sample is relatively small, OIII and D/B are obviously not tightly correlated since the first group contains some disky galaxies with low star formation and group 3 has bulgy galaxies with OIII in emission. Anyhow, this example is only for illustration, and does not pretend to generality. At least, it shows that since the formation scenarios for galaxies are complicated, they cannot be studied with two parameters only, like the present morphologies (like D/B) and the OIII signature.

As a conclusion of this study, the level of contradiction of the minimum contradiction matrix furnishes an objective measure of the deviations to a tree or network structure. This measure can be used as a criterion to select the most important characters. The minimum contradiction matrix can also be of great help to a posteriori discretize continuous characters so that the resulting discrete characters can be well described by a tree or a split network. Quite interestingly, while discretization of continuous characters is often problematic, the posterior discretization after ordering of the taxa can help removing contradictions from a split network or tree structure. We believe there is here a great hope not only for astrophysics but for many evolutionary and multivariate studies.

COMMENTS AND CONCLUSION

It is now clear that cladistics can be applied and be useful to the study of galaxy diversification. Many difficulties, conceptual and practical, have been solved,. Significant astrophysical results have been obtained and will be extended to larger samples of galaxies and globular clusters. However, many paths remain in the exploration of this new and large field of research.

There are difficulties that seem to be intrinsic to astrophysics. Most notably, we have millions of objects, but a few tens of descriptors. Of course the situation will improve with time, in particular with integral-field spectroscopy (spectra of detailed regions in the galaxy, e.g. Ensellem et al. 2007). Spectra might not be currently employed at their full capacity of description. Here also, sophisticated statistics can help. But it is not clear whether these data would lead to the discrimination of hundreds of classes. Perhaps this is an erroneous target, perhaps galaxies cannot be classified with such a refinement. Nonetheless, this is already a matter of using multivariate clustering methods and interpreting their results usefully. This is probably more than a conceptual question because clustering with continuous data and large intraspecific variance is a very complicated problem in itself. Sophisticated statistical tools must be used, but the question of characterizing the groups in this context requires a different culture that is not yet widespread in astrophysics. We are convinced that improvements can be made here.

Cladistics is supposed to identify clades, that are evolutionary groups, whereas the concept of "species" is not defined at all in astrophysics, and we have even not converged toward groupings based on multivariate analyses. Understanding the formation and evolution of galaxies, depicting a global scenario for galaxy diversification, is a major challenge for the astrophysics of XXIst century.

Figure 1 :

 1 Figure 1: Cladogram of Galactic Globular Clusters.

Figure 2 :

 2 Figure 2: Scatterplots for the 4 parameters used in the analysis, colors corresponding to groups defined on the cladogram of Figure 1.

Figure 4 :

 4 Figure 4: Left: Character values (top: logs, middle: Kabs, bottom: ldiam) vs taxa arranged according to the best order. Right: Distance matrix on the 3 characters.

Figure 5 :

 5 Figure 5: Minimum contradiction matrix for the characters D/B and OIII.

Figure 6 :

 6 Figure 6: An example of a split network.

Figure 7 :

 7 Figure 7: Behavior of the two variables D/B and OIII arranged according to the best order represented in Figure 5. The a posteriori discretization is shown, and yields the tree in Figure 8.