Prof Michele Usseglio-Poltera

Prof Young-Seuk Tackx

Dr Nam Park

Prof So

Olivier

Thao Sokunthia

Hap Sophorn

Ith Saveng

MEAS Senghun Phauk Sophany

Kheam Sokha

Nicholas Souter

Nicole Hommet

Linda Jalabert

Dominique Pantalacci

Sigrid Schuermans

Marianne Van Den Hove

Veerle Sylvie Bauwens

Gaël Grenouillet

Sébastien Brosse

Loïc Tudesque

Mathieu Chevalier

Pengbun Ngor

Ratha Chea

Ratha Seng

Kong Heng

Bunyeth Chan

Phanith Long

Sophea Iv

Sophany Soeun

Chuanbo Guo

Chao Wang

Mingli Lin

Kevin Cilleros

Marine Cambon

Jade Bruxaux

Lucie Kuczynski

Niels De Troyer

Wout Van Echelpoel

Sacha Gobeyn

Jana Van Butsel

Naomi Damanik Ambarita

Marie Anne

Eurie Forio

Selamawit Tuan Long Ho

Negassa Chawaka

Natalia Carolina

Rubén Donoso Pantoja

Jerves Cobo

Daniel Mercado

Thi Hanh

Tien Nguyen

Tri Trương

Trịnh Từ

Stijn Bruneel

Nathalie Claire Paracueles

Channy

Boets, P., Lek, S, Goethals, P. R Sor

Optimizing the reliability of classification tree models in predicting alien mollusc occurrence: a hindcasting-and forecasting-based approach

Overall aims: Freshwater tropical and temperate river systems are known to support different biotic communities. In this study, I investigated benthic macroinvertebrate community composition and diversity and its spatial and temporal variation both in tropical Asian and temperate European regions. I also examined the influences of physical-chemical water quality variables on community composition, variations and diversity, and modelled the occurrence of selected species.

Locations: Tropical Asia: the Lower Mekong Basin (LMB), covering an area of 609,000 km 2 ; Temperate Europe: Western Europe, Flemish rivers (Belgium), covering an area of 13,787 km 2 .

Materials and Methods:

For the LMB, data collected from 2004 to 2008 were used, and median values of this period were analysed. For Flemish rivers, data collected from 1991 to 2010 were used. The data were divided into 4 periods: D1: 1991-1995, D2: 1996-2000, D3: 2001-2005 and D4: 2006-2010. The medians of each period were used for detailed spatial analyses. Multivariate analyses were applied to relate community composition and diversity to physical-chemical variables. Five modelling techniques namely Logistic Regression (LR), Random Forest (RF), Support Vector Machine (SVM), Artificial Neural Network (ANN) and Classification Tree (CT) were used to model the occurrence of selected species.

Main results:

Community composition variations, diversity and relationship with environmental variables

From the LMB, 299 macroinvertebrate taxa belonging to 196 genera and 90 families were identified: 131 insects, 98 molluscs, 38 crustaceans, and 32 annelids. These are the largest numbers ever reported for the LMB. Alien taxa were not assessed due to the lack of distribution data and the limited taxonomic knowledge. Insects were mostly found at the upstream parts, which were characterized by high altitude, clear water with high concentration of dissolved oxygen (DO). Molluscs, crustaceans and annelids dominated the downstream parts, which were characterized by a large surface area of watersheds, deep and wide rivers and high water temperature. These diverse compositions found at different spatial scales resulted in an increased local diversity (Shannon-Weiner, H') from the upstream (H' = 1.9) to the downstream (H' = 2.7) parts, and in a great amount of total variation, or beta (β) diversity (BDTotal = 0.8, on a 0-to-1 scale). When computed separately for each component community (i.e. insects, molluscs, crustaceans and annelids), molluscs and insects had a higher variation and their LCBDs greatly contributed to β diversity of global communities, whereas crustaceans and Summary -English vi annelids had a lower variation and contributed less to the β diversity. A high degree of uniqueness in community composition in this tropical system mostly occurred along the mainstream of the rivers, which are highly associated with anthropic disturbance.

From Flemish rivers, 207 macroinvertebrate taxa belonging to 145 families were identified: 131 insects, 34 molluscs, 21 crustaceans and 21 annelids. Seventy-three alien macroinvertebrate taxa were collected. From the past (D1) to the recent period (D4), most taxa and their abundance were linked to high values of DO and low values of other water quality variables. However, Chironomidae thummi-plumosus, Naididae and Asellidae were found to be negatively associated with DO concentration, but positively related to the concentration of ammonium (NH4 + ), phosphate, chemical oxygen demand (COD) and nitrate (NO3 -). Local diversity increased from D1 (H' = 1.6) to D3 (H' = 1.9), but decreased in D4 (H' = 1.6), which is due to a tremendous increase in overall abundance of macroinvertebrates in the latest period when the water quality was substantially improved. This may indicate that water quality improvement is favourable for some dominant taxa, and thus increases the homogenization of the communities and subsequently reduces the local diversity. Beta diversity of global communities was moderate, but significantly increased from the past (BDTotal = 0.50) to the most recent period (BDTotal = 0.59). When the water quality was poor (i.e. in the 1990s), a significant seasonal diversity (α and β) was observed; the diversity was low in Spring compared to Summer and Autumn. However, the seasonal difference was not found in the 2000s during which the rehabilitation program for improving water quality in most of parts of Flanders was step by step implemented. When computed separately for each component community, molluscs and insects had a high β diversity which always increased from D1 to D4, compared to β diversity of annelids (the lowest) and crustaceans which fluctuated between D1 and D4. LCBDs of the four components significantly and highly contributed to the β diversity of global communities. A high degree of uniqueness in community composition in this temperate system was more related to sites situated in the main harbour watercourses and in brackish polders, where high values of EC and pH were recorded.

Species occurrence modelling

For the LMB, the occurrence of 199 species could be predicted using LR, RF, SVM and ANN.

The four modelling techniques yielded significantly different performances (p<0.01), of which ANN yielded the highest performance and was found to better predict the occurrence of rare species. For Flemish rivers, the occurrence of alien molluscs and their co-occurrence with native molluscs were predicted using CT. Based on field data from D1 to D4, the CT models Summary -English vii were only able to reliably predict the "co-occurrence", but not the sole occurrence of "alien" molluscs. The co-occurrence was mainly dependent on sinuosity and a set of chemical water quality variables (e.g. NH4 + , NO3 -, COD, pH). When the CT models were optimized by incorporating field and cloned data (i.e. a dataset obtained by independently duplicating the field data points by k different individuals) via hindcasting and forecasting models, the occurrence of alien molluscs was correctly predicted with a low error rate. This result corresponds to field observations, where alien mollusc occurrence has been observed over the last two decades in Flanders.

Main conclusion and remarks:

The environmental conditions of the two systems certainly favour a different macroinvertebrate community composition, and thus lead dissimilar variation. The LMB was found to support higher diversity compared to Flemish rivers. This could be due to the fact that most invertebrates from the LMB were identified to the species level, while the invertebrates from Flemish rivers were identified only to family or genus level.

However, these findings indeed are not revealing the real composition and diversity for the LMB because it is very large, compared to Flemish rivers, but has been hardly studied. The number of reported taxa from the LMB is thus most likely to be underestimated. On the other hand, Flemish rivers have been extensively investigated and monitored regularly. Nevertheless, the two systems showed some similarities: molluscs and insects had a higher total variation, compared to crustaceans and annelids. Moreover, a high degree of uniqueness in community composition of the two systems mostly occurred at sites where a high level of anthropic disturbance was observed. Among the measured environmental variables, altitude, surface area of watersheds, river width and depth, Secchi depth, DO, EC and water temperature were the key factors affecting macroinvertebrate community composition and diversity in the LMB; whilst in Flemish rivers, DO, EC, pH, NH4 + , NO3 -, COD, phosphate and sinuosity were the key variables. These findings can provide useful information and insights which could be used to support management, conservation and restoration planning in each system. Among modelling techniques applied in the LMB, ANN performed the best, and yielded better results when predicting the occurrence of rare species. The prediction of the occurrence of alien molluscs in Flemish rivers was successfully optimized using CT models. Whether in the past or recent periods, the results of this optimization correspond to field observations. To test for transferability, the successfully optimized models are suggested to be validated using data collected outside Flanders.

Financial means are very crucial to make this scientific mission possible. I am therefore grateful to the Erasmus Mundus LOTUS Unlimited project (European Commission) and the Belmont forum TLSCC projects for providing the chance to do this research. And also, I would like to thank the Mekong River Commission (MRC) and the Flemish Environment Agency (VMM) for the opportunity to study their samples.

My thankful appreciation owe to Prof Sovan LEK and Mrs Sithan LEK for their hospitality, kindness, care, and delicious food. They consider me like their own child, and their house feels like a home for me in Toulouse. I do enjoy gardening (especially during the weekend!), shopping and travelling with them. Many thanks to my second parents and their family in Belgium and my second mum (Judy Santmire) in the USA for their warm welcome, treats and talks during my stressing moments, and to brothers and sisters in churches (in Cambodia, France, and Belgium) for their supportive prayers. Highest thanks be to God for he has given me this wonderful scientific mission and for his wisdom and blessing so that I can complete this assignment with perseverance and joy.

Grateful thanks is extended to my beloved parents, brother and sister. It is blessing to be born and grow up in such a warm and lovely family.

General Introduction 1.1 Background to the study

Invertebrates, defined as organisms without backbones, are the majority of the global animal species, with an estimation of ~98% [START_REF] May | How many species are there on Earth?[END_REF][START_REF] Chapman | Numbers of living species in Australia and the World[END_REF]). The invertebrates that can be seen without the aid of a microscope are referred to "macroinvertebrates", which mainly comprise insects (an estimation of ~73%), arachnids (~9%), nematodes (~7%), molluscs (~3%), crustaceans (~2%), flatworms (~1%) and annelids (~0.4%) [START_REF] Chapman | Numbers of living species in Australia and the World[END_REF]IUCN 2014). In freshwater ecosystems, insects, molluscs, crustaceans and annelids are the most commonly found and most studied, compared to the other groups of invertebrates (Collier and Lill 2008;Arab et al. 2009;Pérez-Quintero 2011;[START_REF] Szöcs | Effects of anthropogenic salinization on biological traits and community composition of stream macroinvertebrates[END_REF]Sor et al. 2017a).

Macroinvertebrates are a key component of freshwater ecosystems (Palmer et al. 1997;Bogan 2008). Firstly, they are considered as ecosystem engineers [START_REF] Jones | Organisms as ecosystem engineers[END_REF][START_REF] Chowdhury | Ecosystem engineering by mussels supports biodiversity and water clarity in a heavily polluted lake in Dhaka, Bangladesh[END_REF]. For example, crustaceans and insects are responsible for regulating decomposition, shredding detritus and bioturbation. Molluscs also contribute to bioturbation, sediment formation and filtering of water, while most annelids regulate decomposition and autotrophs although some of them also promote bioturbation and sediment formation (Palmer et al. 1997).

Secondly, macroinvertebrates are key networks of food chains that are important in maintaining freshwater and terrestrial food webs (Fig. 1.1). Molluscs, crabs, shrimps and other benthos (e.g.

annelids and insects) are the food sources for reptiles (e.g. turtles), numerous carnivorous fish species (Chea et al. 2016), and some terrestrial predators including birds [START_REF] Poulsen | Distribution and Ecology of Some Important Riverine Fish Species of the Mekong River Basin[END_REF]).

Moreover, freshwater macroinvertebrates are generally used as bioindicators and bioassessment. For instance, the presence of mayflies, caddisflies and stoneflies (insects) can indicate good water quality conditions, while the presence of annelids, clams, pouch snails, water bugs and crayfish can indicate a moderate to high level of pollution in the water (Feld and Hering 2007;Collier and Lill 2008;Królak and Korycińska 2008;Wang et al. 2012). Due to these facts, macroinvertebrates have been used for bioassessment, which is useful to support management, restoration and conservation planning in freshwater ecosystems [START_REF] Heino | Assessing physical surrogates for biodiversity: do tributary and stream type classifications reflect macroinvertebrate assemblage diversity in running waters?[END_REF]Kudthalang and Thanee 2010). Structure and spatial patterns of lotic macroinvertebrate communities are known to organize along a longitudinal downstream river/stream gradients (Vannote et al. 1980;Friberg et al. 2010). In the upstream parts, the communities are mainly characterized by a high abundance of insects (Collier and Lill 2008;Arab et al. 2009;[START_REF] Jiang | Longitudinal patterns of macroinvertebrate communities in relation to environmental factors in a Tibetan-Plateau river system[END_REF], whereas the downstream communities are characterized by a high abundance of molluscs, crustaceans and worms (Arscott et al. 2005;Collier and Lill 2008;Pérez-Quintero 2011). For pristine or undisturbed rivers, structural and functional macroinvertebrate communities are adapted to conform to the most suitable positions, depending on the available energy flow and environmental variability, throughout the river's continuum (RCC, Vannote et al. 1980). According to RCC, communities in the headwaters, collecting energy from decayed leaves, needles and stems, are mostly composed of shredders, collectors and less grazers/scrapers. The mid-reach, being strongly exposed to sunlight, supports more grazers/scrapers and collectors, and the lower-reach, having a low photosynthesis production in the rivers (due to high turbidity and surface film) and high energy inputs (mostly from upstream sources), is home to numerous collectors. However, in large/floodplain rivers, which receive a high level of disturbance, the RCC cannot be applied to address the biological systems. This is because differences between biotic community composition in these rivers are determined by spatial and temporal heterogeneity along the rivers (Sedell et al. 1989), and by natural and human-derived disturbances [START_REF] Clarke | Macroinvertebrate diversity in headwater streams: a review[END_REF][START_REF] Muñoz | Bridging levels of pharmaceuticals in river water with biological community structure in the Llobregat River basin (Northern Spain)[END_REF].

Distribution, composition and diversity patterns of macroinvertebrates greatly vary depending on studied climatic regions (e.g. tropical vs temperate), zoogeographic regions (e.g. Palaearctic, Nearctic, Neotropical, Afrotropical, Oriental and Australasian) and geographic regions (e.g.

Asia, Europe, America and Africa) [START_REF] Martin | Global diversity of oligochaetous clitellates ("Oligochaeta"; Clitellata) in freshwater[END_REF][START_REF] Yeo | Global diversity of crabs (Crustacea: Decapoda: Brachyura) in freshwater[END_REF]Bogan 2008;Ferrington 2008). The variation in composition and diversity found from each climatic and zoogeographic region may reveal the different favourable environmental conditions for diverse taxonomic groups to live on (Dudgeon et al. 2006;[START_REF] Boulton | Are tropical streams ecologically different from temperate streams?[END_REF]. However, the ecological processes in these systems appear to be driven by more or less the same variables, e.g. drought, disturbance, nutrient concentration and trophic structure [START_REF] Boulton | Are tropical streams ecologically different from temperate streams?[END_REF][START_REF] Dudgeon | Tropical Stream Ecology[END_REF]).

Nevertheless, geographic regions that lie in the tropical zone harbour a higher biodiversity, at least for most invertebrate taxa, than those lie in the southern or northern temperate zone (Sodhi et al. 2004;[START_REF] Boulton | Are tropical streams ecologically different from temperate streams?[END_REF]). However, stream invertebrate ecology in many parts of tropical region (e.g. South America, Africa, and Asia) remains little investigated, whereas stream invertebrates in temperate regions (e.g. North America, Europe, Australia and New Zealand)

have been well studied [START_REF] Dudgeon | Tropical Stream Ecology[END_REF]Boyero et al. 2009).

Macroinvertebrates in Asian and European rivers: a general overview

Taxonomic and ecological knowledge on freshwater macroinvertebrates in Asia, as mentioned earlier, is still limited [START_REF] Boulton | Are tropical streams ecologically different from temperate streams?[END_REF]Boyero et al. 2009). Most research in the tropical Asia is largely restricted to a few geographic regions including the Hong Kong and Peninsular Malaysia [START_REF] Resh | Multinational, freshwater biomonitoring programs in the developing world: lessons learned from African and Southeast Asian river surveys[END_REF]Boyero et al. 2009;[START_REF] Leung | Scales of spatiotemporal variability in macroinvertebrate abundance and diversity in monsoonal streams: detecting environmental change[END_REF]Al-Shami et al. 2013). A wider range of studies has also been revealed from the northern subtropical Asian rivers, e.g.

Yangtze and the Upper Mekong River or the so-called Lancang River in China. Research topics from these river basins include species distribution, spatio-temporal patterns and species records [START_REF] Nieser | A new subgenus and six new species of Nepomorpha (Insecta: Heteroptera) from Yunnan, China[END_REF][START_REF] Shao | Macroinvertebrate community structure in Three-Gorges Reservoir, China[END_REF][START_REF] Qi | Convergent, parallel and correlated evolution of trophic morphologies in the subfamily schizothoracinae from the qinghai-tibetan plateau[END_REF]. However, most of the research findings are reported in Chinese and are not publically available; only a few are accessible, e.g.

the benthic macroinvertebrates as indicators of ecological status in Yangtze River [START_REF] Pan | An exploratory analysis of benthic macroinvertebrates as indicators of the ecological status of the Upper Yellow and Yangtze Rivers[END_REF]) and the seasonal variability of metazooplankton (including crustaceans) communities and new mollusc species records from the Upper Mekong Basin (Du et al. 2011;[START_REF] Wu | Seasonal variability of water quality and metazooplankton community structure in Xiaowan Reservoir of the upper Mekong River[END_REF].

For the Lower Mekong Basin, more investigations have been recently conducted in Thai streams [START_REF] Boonsoong | An implementation plan for using biological indicators to improve assessment of water quality in Thailand[END_REF]Kudthalang and Thanee 2010;[START_REF] Phaphong | Development of a benthic macroinvertebrate biotic index to evaluate wetland health in Northeastern Thailand[END_REF]David and Boonsoong 2014) and recently also some Philippine streams have been investigated [START_REF] Tampus | Water quality assessment using macroinvertebrates and physico-chemical parameters in the riverine system of Iligan City, Philippines[END_REF][START_REF] Sinco | Riverine biota as indicators of water quality in tropical Cagayan de Oro River, Philippines[END_REF][START_REF] Fajardo | Aquatic macroinvertebrate diversity and physicochemical characteristics of freshwater bodies in Tubay, Agusan Del Norte, Philippines[END_REF][START_REF] Magbanua | Water physicochemistry and benthic macroinvertebrate communities in a tropical reservoir: the role of water level fluctuations and water depth[END_REF][START_REF] Forio | Assessment and analysis of ecological quality, macroinvertebrate communities and diversity in rivers of a multifunctional tropical island[END_REF]. Most of these studies are related to species diversity, description of new species and using benthos to assess water quality in river systems [START_REF] Parnrong | New records of Behningiidae, Potamanthidae, and Prosopistomatidae (Ephemeroptera) from Thailand[END_REF][START_REF] Sangpradub | New records of Heptageniid mayflies Asionurus and Thalerosphyrus (Ephemeroptera: Heptageniidae) from Northeastern Thailand[END_REF][START_REF] Flores | Macroinvertebrate composition, diversity and richness in relation to the water quality status of Mananga River, Cebu, Philippines[END_REF]. Macroinvertebrates from other geographic areas including Myanmar, Laos, Cambodia and Vietnam, which mainly share the Lower Mekong Basin, remain very scarce.

On the contrary, knowledge on macroinvertebrates from river systems in Europe has been extensively studied (Boyero et al. 2009). Since the adoption of the European Water Framework Directive (WFD) (European Commission 2000), freshwater macroinvertebrates have become the central focus [START_REF] Pollard | The European water framework directive: a new era in the management of aquatic ecosystem health?[END_REF][START_REF] Hering | The European Water Framework Directive at the age of 10: a critical review of the achievements with recommendations for the future[END_REF]. Macroinvertebrates from hundreds of streams have been studied and used to assess water quality [START_REF] Buffagni | A Europe-wide system for assessing the quality of rivers using macroinvertebrates: The AQEM Project and its importance for southern Europe (with special emphasis on Italy)[END_REF][START_REF] Verdonschot | Testing the European stream typology of the Water Framework Directive for macroinvertebrates[END_REF]. Within 10 years of the implementation, ~1,900 papers resulted from research projects associated with WFD [START_REF] Hering | The European Water Framework Directive at the age of 10: a critical review of the achievements with recommendations for the future[END_REF]. This results in a very well documented knowledge on freshwater macroinvertebrates and their ecological applications for Europe. Furthermore, a diverse assessment methods have been developed [START_REF] Birk | Three hundred ways to assess Europe's surface waters: An almost complete overview of biological methods to implement the Water Framework Directive[END_REF], some of which have applied a predictive modelling framework that is based on macroinvertebrates or use environmental variables to predict future distribution, occurrence and abundance of particular taxa (Goethals et al. 2007;[START_REF] Everaert | Development and assessment of ecological models in the context of the European Water Framework Directive: Key issues for trainers in datadriven modeling approaches[END_REF][START_REF] Boets | Evaluation and comparison of data-driven and knowledge-supported Bayesian Belief Networks to assess the habitat suitability for alien macroinvertebrates[END_REF].

When macroinvertebrate composition and diversity are related to measured environmental variables, key factors driving spatio-temporal changes have been known to be more or less the same regardless of geographic or climatic regions. For instance, macroinvertebrate communities in river basins from southern China (Pearl, Yangtze and Qiangtang rivers), from northern Portugal (the Olo, Corgo, Pinhao and Tua rivers) and from Susquehanna River (New York, North America) have been reported to be influenced by land use types including anthropogenic disturbance (Allan 2004;[START_REF] Bruns | Macroinvertebrate response to land cover, habitat, and water chemistry in a mining-impacted river ecosystem: A GIS watershed analysis[END_REF]Cortes et al. 2011;Cortes et al. 2013).

Another example can be found from European Mediterranean (Evrotas River, Greece) and Asian streams (Peninsular Malaysia) that stream size (e.g. width and depth), dissolved oxygen and pH were the key factors influencing macroinvertebrate composition and variation (Al-Shami et al. 2013;Salmah et al. 2014;[START_REF] Karaouzas | Spatial scale effects on Chironomidae diversity and distribution in a Mediterranean River Basin[END_REF]. These indicate that similar ecological processes can be expected from different ecological systems.

Modelling techniques and applications

Various modelling techniques have been widely and increasingly implemented in ecological systems (Lek et al. 1996;Park et al. 2003;[START_REF] Schröder | Statistical classification of terrestrial and marine ecosystems for environmental planning[END_REF]Lencioni et al. 2007;Guo et al. 2015). The techniques applied are generally used to explain and predict the relationship between the occurrence or abundance of studied species and environmental variables (Goethals et al. 2007;Boets et al. 2013). Utilization of modelling techniques to combine both explaining and predicting such relationships is also commonly applied (Roura-Pascual et al. 2009;Call et al. 2016). Applications of predictive models have provided knowledge and improved the understanding of the ecology and behaviour of studied taxa, which could be used to support decision making, management and conservation planning. For instance, many previous studies have used predictive models to predict the occurrence and distributional areas of plants, herbs, macroinvertebrates and fish (Thuiller et al. 2005;Roura-Pascual et al. 2009;Vicente et al. 2011;Boets et al. 2013;Chen et al. 2015;Guo et al. 2015).

However, the application of predictive models has been suggested to be carefully taken into account because they can have a wide variation in performance (Segurado and Araujo 2004;Elith et al. 2006;Guisan et al. 2007). Some models even yield contrasting predictions of habitat suitability (e.g. Guisan et al. 2007;Evangelista et al. 2008;Roura-Pascual et al. 2009).

Furthermore, predictive models are sensitive to parameterization and selection criteria during the modelling process (Araújo and Guisan 2006;Elith et al. 2006), and thus can result in an uncertainty of current or past/future projections of species distributions (Svenning et al. 2008;Buisson et al. 2010). Due to this fact, when calibrating and validating predictive models, carefully taking into account the data characteristics (e.g. sample size, species prevalence or environmental predictors), parameterization and selection criteria are usually recommended (Luoto et al. 2006;Dormann et al. 2008).

Research problem, aims and objectives

The Lower Mekong Basin (LMB), which includes portions of Thailand, Laos, Cambodia and Vietnam, is characterized by a long and large floodplain (Eastham et al. 2008) and is known for its high biodiversity (Sodhi et al. 2004). However, the knowledge of macroinvertebrates in the LMB is poorly investigated. Given that this river basin is being impacted by various anthropogenic disturbances such as agricultural activities, aquaculture, urbanization and mining (Sodhi et al. 2004;Nhan et al. 2007;Köhler et al. 2012), there is an urgent need to study the patterns of spatial organization, community structure and variations of macroinvertebrates in this basin and their relation to environmental factors. Up to date, only a few studies (except for those conducted in Thailand) have been published on the basin, e.g. community structure and composition of littoral invertebrates in the Mekong delta (Wilby et al. 2006) and the diversity and distribution of crustaceans and molluscs in the Indo-Burma region (Cumberlidge et al. 2011;Köhler et al. 2012). Yet, no attempt has been made to examine the large spatial patterns, community structures, variations (i.e. β diversity) of macroinvertebrate communities and their relation to key environmental variables nor the application of predictive modelling in this hardly studied basin.

On the other hand, river systems in Europe as well as in Flanders suffered from severe water quality degradation in previous decades (e.g. from 1980s to 1990s). During these periods, some native species were reported to disappear (Bernauer and Jansen 2006) and only those that were able to withstand the water quality degradation remained. At the same time, most European river systems have been exposed to a number of alien macroinvertebrate species (Leuven et al. 2009;Boets et al. 2016). From the 2000s until now, the water quality of European rivers has been greatly improved. This water quality improvement does not only promote the occurrence and abundance of native species, but also favours the alien species to spread widely, which consequently may lead to changes in community composition. As such, investigation spatiotemporal changes in community composition, variations and predicting the occurrence of alien species across Flemish rivers, which have been poorly studied, will provide insights into the ecology of overall communities and of studied alien species. Results from this investigation can be used to support management and conservation planning.

The aims and specific objective (or questions) of the present study are:

Aim 1. Investigating general patterns of macroinvertebrate communities and their relation to environmental variables in the two systems, i.e. the LMB and Flemish rivers.

1.a. Investigating patterns of spatio-temporal variation in macroinvertebrate assemblages/communities.

1.b. Analysing the variability of macroinvertebrate composition among the assemblages/communities, and determining key indicator/important taxa (the most representative taxa/taxa with high among-site variance).

1.c. Identifying the important environmental variables that are associated with the particular macroinvertebrate assemblages/communities. Aim 2. Determining the total variation in macroinvertebrate communities (i.e. total β diversity) and the key determinants in the two systems.

2.a. Is there a moderate or a large amount of total β diversity?

2.b. What are the taxa that contribute most to the total β diversity?

2.c. What are the environmental conditions and component communities (e.g.

annelids, crustaceans, insects and molluscs) that significantly influence the total β diversity? Aim 3. Predicting macroinvertebrate species occurrence and analysing the performance of modelling techniques applied in the LMB.

3.a. Predicting the occurrence of macroinvertebrate species and comparing the performance of the applied techniques based on a complete prevalence range (i.e. 0.0-1.0), and different prevalence ranges (i.e. at a 0.1 interval).

3.b. Analysing how the species prevalence affects the behaviour of modelling techniques' performance.

Aim 4. Predicting the occurrence of alien species and their co-existence with native species and identifying the key determining variables in Flemish rivers over the past two decades (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010).

4.a. Identifying key determining physical-chemical variables associated with the occurrence of alien species (i.e. alien molluscs) and with the co-occurrence of alien and native species, using a classification tree modelling technique.

4.b. Optimising the reliability of classification tree models in predicting alien mollusc occurrence.

Materials and Methods

Case study in the LMB and dataset

The LMB

The Mekong River Basin is divided into the Upper Mekong Basin (UMB) and the Lower Mekong Basin (LMB). The UMB on the Tibetan plateau in China is composed of narrow, deep gorges and small, short tributaries, whereas the LMB stretches from Yunnan province in South

China to the delta in Vietnam and it covers approximately 70% of the total length of the whole basin (Eastham et al. 2008). The LMB consists of a large floodplain and long, broad tributaries and it drains more than 76% of the Mekong basin. The climate of the LMB is dominated by a tropical monsoon rainfall system, which is characterized by a dry (November -April) and a wet (May -October) season generated by the northeast monsoon and the south-west monsoon, respectively. The most intensive rainfall falls from July to September, while the lowest precipitation is observed between January and April (Adamson et al. 2009). The annual rainfall of the LMB varies from 1,000 -1,600 mm in the driest regions to 2,000 -3,000 mm in the wettest regions (Hoanh et al. 2003). A higher precipitation is found in the eastern mountainous regions of Laos and in northeast Thailand (Eastham et al. 2008).

The largest floodplain water body of the LMB is the Tonle Sap Lake (TSL) in Cambodia (Adamson et al., 2009), which is the largest freshwater lake in Southeast Asia (Sarkkula et al. 2003). The TSL is connected to the Mekong through the Tonle Sap River, and thus creating an exceptional hydrological cycle. In the wet season, the TSL receives excess water from the Mekong River and expands its surface area from 2,500 km 2 to 15,000 km 2 . In the dry season when the rain ceases and water levels drop in the Mekong, a reverse flow occurs; the drained water from the TSL flows to the Mekong delta (Arias et al. 2011). The Mekong delta is characterized by a number of man-made canals, which are mostly used for domestic and agricultural activities (Kummu et al. 2008).

Data collection and processing

Benthic macroinvertebrates were sampled at 60 sampling sites along the main channel of the LMB and its tributaries by the Mekong River Commission (MRC) (Fig. 2.1). This sampling was carried out once a year in March during the dry season from 2004 to 2008. At each sampling site, macroinvertebrates were sampled from three locations in the benthic zone: near the left and right banks, and in the middle of the rivers. At each location, a minimum of three samples (where inter-sample variability is low, e.g. tributaries) to a maximum of five samples (where inter-sample variability is higher, e.g. the main channel and the delta) were collected using a

Petersen grab sampler which has a sampling area of 0.025 m 2 . With the grab sampler, four subsamples were taken and pooled to give a single sampling unit covering a total area of 0.1 m 2 .

In total, between nine (3 samples × 3 locations) and fifteen (5 samples × 3 locations) pooled samples were collected at each sampling site. Each pooled sample was rinsed using a sieve (0.3 mm mesh size). In the field, the samples were sorted and then preserved by adding 10%

formaldehyde to obtain a final concentration of about 5%. In the laboratory, they were identified to the lowest taxonomic level possible and counted using a compound microscope (40 -1,200 magnification) or a dissecting microscope (16 -56 magnification). Macroinvertebrate abundance data per sampling unit was averaged across all samples (between 9 and 15 samples) collected from each sampling site. At the sampling site, geographic coordinates and altitude were determined with a GPS (Garmin GPS 12XL). River width was measured in the field using a Newcon Optik LRB 7x50 laser rangefinder. Other physical-chemical variables were measured at the three locations where macroinvertebrates were sampled. River depth was measured using a line metre. With a handheld water quality probe (YSI 556MP5), water temperature, dissolved oxygen, pH and electrical conductivity were measured at the surface (0.1-0.5 m) and at a depth of 3.5 m or at a maximum depth of the river (wherever less than 3.5 m) and then the average value was recorded for each location. Water transparency was measured with a Secchi disc by lowering it into the water and recording the depth at which it was no longer visible. The physical-chemical data of each sampling site was the averaged value across the three sampling locations. Distance from the sea and the surface area of watersheds drained at each sampling site was determined using a Geographic Information System (ArcGIS 10.0, ESRI). Geographic data (ArcGIS shapefiles) about the LMB (river networks, basin boundaries, land covers, and subcatchments derived from topographical maps) was provided by the MRC.

In total, 108 samples were collected from the 60 sampling sites. In 2008, 3 sampling sites were sampled further away from their original sampling coordinates, and thus were considered as different sampling sites (see Appendix T1). Therefore, a total of 63 sampling sites were taken into account in the analyses. Because of unequal sampling efforts (i.e. unequal and different number of samples at each site during the 5-year sampling period) and missing values of environmental variables, we used median values from the collected data to represent each site in the analyses, as suggested by McCluskey and Lalkhen (2007). These median values were used in all of the analyses corresponding to the case study of the LMB.

Case study in Flemish rivers and dataset

Flanders

Flanders (northern Belgium) is located in Northwest Europe and its Northwestern part is bordered by the North Sea (Fig. 2.2A). Flanders has a total area of 13,522 km 2 , and is considered as one of the most densely populated regions in Europe (477 inhabitants/km 2 in 2015, https://en.wikipedia.org/wiki/). Flanders is classified as a lowland area, which is divided into different rivers basins (Fig. 2.2B). This region is influenced by a temperate oceanic climate, as same as most of northwestern European countries are (e.g. UK, France, Luxembourg, Netherland and Denmark) [START_REF] Peel | Updated world map of the Köppen-Geiger climate classification[END_REF]). Flanders has a dense watercourse network including navigable canals. Agriculture, industry and residential areas are the main land use types of Flanders and its landscape is characterized by highly fragmented and complex mosaic of land use types [START_REF] Poelmans | Detecting and modelling spatial patterns of urban sprawl in highly fragmented areas: a case study in the Flanders-Brussels region[END_REF]. This fragmentation and complexity may have put a high pressure on habitat quality and biodiversity in Flanders.

Data collection and processing

The Flemish Environment Agency (VMM) has collected biological and environmental data in Flanders since 1989. The monitoring sites include all types of watercourses from all river basins. Every three year from the beginning, a fixed set of sampling locations was sampled.

Most of the sampling locations were only sporadically sampled, and thus results in a large dataset of more than 11,000 biological samples collected at more than 2500 sites spread over different water bodies (Fig. 2.2C). In this monitoring program, the sampling protocol was entirely based on the method as described by Gabriels et al. (2010). Macroinvertebrates were collected using a standard handnet, which is made of a metal frame (0.2 m by 0.3 m) to which a conical net is attached with a mesh size of 300 µm. The kick sampling was made along the watercourses at a stretch of approximately 10-20 m. Each sample was collected for three minutes for small watercourses (less than 2 m wide) or five minutes for larger rivers. At sampling sites where the kick sampling method was not possible, artificial substrates were used.

Three replicates of artificial substrates, which consisted of polypropylene nets (5 litres) filled with bricks of different sizes, were left in the water for a period of at least three weeks after which they were retrieved. Leaving this period enables species to colonize the substrates. The different sampling efforts of the two sampling approaches (the kick and artificial substrate sampling) may have repercussion on the diversity of sampled invertebrates. However, according to Gabriels et al. (2010), the two approaches are standardized semi-quantitative methods and are similar in terms of sampled macroinvertebrate abundance. In the laboratory, macroinvertebrates in the VMM database were identified to the level (family or genus) needed for the calculation of the biotic water quality index.

Electrical conductivity (EC), pH and dissolved oxygen (DO) were measured in the field with a hand-held probe (Cond 315i, oxi 330, wtw, Germany and 826 pH mobile, Metrohm, Switzerland). All additional chemical variables, i.e. ammonium (NH4 + ), chemical oxygen demand (COD), biological oxygen demand (BOD), total phosphorus (Pt), nitrate (NO3 -), nitrite (NO2 -), Kjeldahl nitrogen, orthophosphate (oPO4), were retrieved from the monitoring dataset compiled by the VMM and which is online accessible (www.vmm.be). Nutrient analysis was performed spectrophotometrically in accordance to ISO 17025. GIS software (version 9.3.1) applied on the Flemish Hydrographic Atlas was used to determine the slope and the sinuosity of a watercourse at a different height in between two points (1000 m apart) and on a stretch of Data from 1991 to 2010 was used for the analyses. Based on the preliminary data mining, the overall communities had a temporal change (especially a somewhat different community composition for the late 2000s) which could be grouped based on a five-year interval.

Therefore, the data were then divided into 4 periods. Each period consisted of samples from a five-year sampling effort (i.e. D1: 1991-1995, D2: 1996-2000, D3: 2001-2005 and D4: 2006-2010). This division can provide useful information on changes in community composition for each period. To analyse spatial variation in the community composition, the median values were used to represent each site for each period. This is because, as mentioned earlier, only a fixed set of sampling locations was sampled regularly whereas most of the other sampling locations were sporadically sampled. For the modelling purposes (see the detailed in the "Modelling" section below), all collected samples were used.

Statistical analyses and modelling approaches

I analysed only four groups of macroinvertebrates in this study. This is because they are the most commonly studied animals and are generally used as bioindicators and assessment in freshwater ecosystems (Feld and Hering2007;Collier and Lill 2008;Wang et al. 2012). The four groups included annelids, crustaceans, insects and molluscs. These four groups were designated as component communities in following paragraphs onwards. All the applied statistical analyses and modelling approaches were performed using functions of packages in the R language program (R Core Team 2013).

Communities clustering and diversity measures

Samples were clustered based on the Bray-Curtis dissimilarity of macroinvertebrate abundance data by using Ward's hierarchical method. The Bray-Curtis dissimilarity distance (Legendre and Legendre 2012) between the macroinvertebrate samples was calculated using the Hellinger transformation in the package vegan of R (Rao 1995).

The macroinvertebrate indicator taxa in each assemblage were determined using the Indicator Value (IndVal, Dufrene and Legendre 1997) with the package labdsv of R (Roberts 2013). The Indicator Value of a taxon is an index ranging from 0 to 1, indicating the least to most important taxa occurring in a group of sites. A value of 1 is obtained when every individual of the taxon is found only in the group and when it occurs at all sites of that group. A high number of taxa with significant Indicator Values may provide information on the habitat they prefer to share.

Taxa having Indicator Values with a p-value ≤0.01 were retained as the most important taxa representing a given assemblage (consisting of a group of sites).

Macroinvertebrate richness, abundance and Shannon-Wiener diversity (H') were calculated for each sampled site, cluster and group. To quantify beta (β) diversity, the community composition data were first Hellinger-transformed (Legendre and Gallagher 2001;Legendre and Legendre 2012). For Hellinger-transformed data, the total variance, or total β diversity (BDTotal), of a community composition data table is an index between 0 and 1, and it can be partitioned into local contribution (LCBD) and species contribution (SCBD) indices. An LCBD value is an index showing the degree of uniqueness in taxonomic composition in each site, computed as the relative contribution of a site to BDTotal, so that the LCBD indices sum to 1, whereas an SCBD index shows the relative degree of variation of a taxon across all sites. The BDTotal, LCBD and SCBD indices were computed using the function "beta.div" available in the adespatial package in R (Dray et al. 2016). The Hellinger transformation was used because the corresponding Hellinger distance is one of the dissimilarity functions admissible for β diversity analyses (Legendre and Gallagher 2001;Legendre and De Cáceres 2013); it does not give high weights to the rare species. In addition to LCBD, Hellinger-transformed data also allow researchers to compute SCBD indices; this is not allowed by most other admissible dissimilarity functions (Legendre and Gallagher 2001;Legendre and De Cáceres 2013). SCBD indices that were higher than the mean of SCBD values identified the taxa that were the most important contributors to BDTotal. In the following paragraphs, BDTotal, LCBD and SCBD designate the indices of the global macroinvertebrate communities, whereas BDATotal, BDCTotal, BDMTotal, BDITotal, and LCBDA, LCBDC, LCBDM and LCBDI designate the BDTotal and LCBD indices for annelid, crustacean, mollusc and insect communities, respectively.

Comparative analyses

Descriptive statistics were used to describe and summarize the information of the collected data.

These included minima, maxima, mean, range, standard deviation (sd), sample size, and percentage. In most cases, mean and standard errors were used to indicate significant differences in macroinvertebrate composition, environmental conditions and model performances between/among groups. Where applicable and appropriate, a one-way ANOVA or a Kruskal-Wallis test was used to test for significant differences between/among unmatched groups (≥3 groups). One-way ANOVA was applied when residuals of the models were normal (Shapiro-Wilk test, p > 0.05, and homoscedastic (Bartlett's test, p > 0.05)); otherwise, the nonparametric test (Kruskal-Wallis) was used. For matched groups (dependent samples), a multifactor ANOVA and a Friedman test were used when the data was normal and not-normal distributed, respectively.

Regression and Multivariate analyses

Simple and multiple regression models were used to access the influence of independent variables (e.g. environments) on response variables (e.g. communities, diversity measures and indices). To identify the strength of the regression models, the stepwise selection with the Akaike Information Criterion (AIC) was applied. The models having the lowest AIC and highest adjusted R 2 were considered to have the strongest influence on the response variables.

Linear Discriminant Analysis (LDA) was performed, using the package ade4 of R (Chessel 2006), to assess which measured environmental variables best accounted for the differences among the macroinvertebrate assemblages grouped by the hierarchical clustering. Before performing the LDA, environmental variables were tested for multivariate homogeneity of within-group covariance (Borcard et al. 2011). The contribution of each variable to the discrimination among assemblages was represented by the standardized factorial coefficient, projected as an arrow on the LDA plot.

With a complex community data, as in the case of Flemish river data, Redundancy Analysis (RDA) was conducted on the Hellinger-transformed abundance data and environmental factors.

RDA is powerful tool for the analysis of community composition data tables (Legendre and Legendre 2012). The RDA model was first tested at global scale to detect for its significance, and afterwards, the forward selection method was carried out in order to select the most importing factors associating with the community composition.

Model development, validation and performance

Five commonly used modelling techniques were applied in this study.

Species occurrence prediction in the LMB

Macroinvertebrate species occurrence was predicted using four modelling techniques: LR, RF, SVM and ANN. The presence/absence of each species was used as the response variable and the measured environmental variables such as altitude, river width, river depth, distance from the sea, water temperature, DO, EC, pH and Secchi depth were used as the input predictors. For every prediction, species that occurred only in one instance (a prevalence of ~0.02) were removed from the data. This is because Leave-One-Out cross-validation (LOO) was used to validate the models, due to small sample size. With LOO, it is not feasible to split data into training and validation sets for species that have only one occurrence instance.

To evaluate the four modelling techniques, Cohen's Kappa Statistic (Kappa), area under receiver operating characteristic curve (AUC), and error rate (ER) of the prediction were used.

The overall performance differences among the four modelling techniques were assessed using a Friedman test (nonnormal distribution). P-values ≤0.05 were considered to indicate significant differences. Each performance measure of the four modelling techniques was then regressed against the species prevalence using two types of models: a linear and a quadratic model. Mean performance values and standard errors were compared to identify which technique better performed based on the complete prevalence and based on different prevalence ranges.

Species occurrence prediction in Flemish rivers

A CT model was used to predict the occurrence of alien and native molluscs and their cooccurrence, and to identify the determining physical-chemical variables. The CT was implemented in this system because it is well suited for analyses of complex ecological data (De'ath and Fabricius 2000). Moreover, it is relatively simple to implement, easy to interpret, and it tolerates missing values during both the training and testing cases (Therneau and Atkinson 1997;De'ath and Fabricius 2000). The prediction was made for each period. The response variable of each sampling site of each period was categorized as: a "native" site (i.e.

a site having only native molluscs present), an "alien" site (i.e. a site having only alien molluscs present) and a "co-occurrence" site (i.e. a site having both alien and native molluscs present).

Due to a limited number of occurrence instances of most alien mollusc genera, all alien genera were merged to form one categorical variable. This provided a higher number of instances for the predictive models and thus a better and more robust development of the model. The results of these predictions could reveal common environmental conditions that most of the alien molluscs prefer. In the same way, all native genera were also merged to form one categorical variable. Physical and chemical water quality variables (i.e. BOD and Kjeldahl nitrogen) that had missing values for more than 5% of the total samples were removed from the analyses.

Therefore, each period consisted of one response categorical variable (native/alien/cooccurrence) and 11 predictor variables. The summary of the physical and chemical water quality variables and of the response variable is shown in Table 2.1. During the model development, all CT trees were initially pruned by setting a complex parameter at cp =0.05. Where the tree had only a root, the cp was lowered to a level (e.g. cp =0.04, 0.03, 0.02, 0.01) that at least two terminal nodes were produced. In practice, the first few splits mostly provide a very informative division of the data (Therneau and Atkinson 1997). These criteria were set in order to make the trees easily interpretable and comparable in terms of the number of variables and complexity. For each period, a three-fold cross-validation was used to train and validate the models. To across the nine models. To identify which variables significantly determine the preferred conditions of each class, the importance of each variable was compared based on the standard error. The same procedures and criteria were applied for the modelling of data of each period.

Optimizing alien mollusc occurrence prediction

Due to the small number of occurrences, the CT models were not able to predict sites where only alien molluscs occurred (see Results section 3.3.2). Subsequently, the CT models were optimized by combining the field and cloned data, and testing two types of resampling approaches (stratified and random) under same and different parameterization settings.

The combined data was prepared as follows. First, a cloned dataset was developed by independently duplicating the original field data of the less frequently occurring classes (mostly the "alien" and "co-occurrence" class) by k different individuals. Then the cloned data was combined with the original field data to obtain the same sample size for each class (Fig. 2.3A).

One may argue to use an equal-stratified dataset of field observations by randomly selecting a number of observations of the more frequently occurring classes in an equal number as there are of the less occurring class. By following this routine, valuable information of many of the unselected field observations will be excluded. That would in this case lead to a small sample size for the model development, which consequently leads to a less reliable and less robust models. On the contrary, using the cloned data to combine with all field observations yields an equal data distribution and sufficient samples for each period. Incorporating the cloned data into the models will increase the maximum likelihood of the prior class distribution, making the models more robust (Lele et al. 2007).

Model optimization was carried out in two scenarios: a hindcasting and forecasting. For the hindcasting, only data period D4 was used to calibrate the models. Then the calibrated models were used to hindcast (validate) the response variable based on the environmental predictors of the data from period D4, D3, D2 and D1. For the forecasting, only data period D1 were employed to calibrate the models. Then the calibrated models were used to forecast (validate)

the response variable based on the environmental predictors of the data from period D1, D2, D3 and D4. During the pruning phase, the complex parameter of all calibrated models was set in the same way for both scenarios, following the procedure described in an earlier section.

The calibration and validation process was based on a three-fold cross-validation (CV), following the procedure described earlier. When splitting the dataset for the 3 fold CV, a random and a stratified resampling approach were applied. The random resampling was made by shuffling and randomly splitting the whole dataset into three folds. From this approach, the sample size of each class was not equally stratified across the three folds. To allow a reliable error estimation of the hindcasting and forecasting models for each period, three replicates of the 3 fold CV were made by reshuffling and randomly splitting the whole dataset again into 3 new folds. For the stratified resampling, the three classes of the response variable and the environmental predictors of the whole dataset were first allocated into three separated data subsets (DS, i.e. DS1, DS2, DS3). Each DS corresponds to the data of each class of the response variable and its corresponding environmental predictors. Then, each DS was randomly divided into 3 smaller data subsets (i.e. DS1a, DS1b, DS1c; DS2a, DS2b, DS2c; DS3a, DS3b, DS3c).

Thus, each of the smaller data subsets of the same class had exactly the same number of samples (n, e.g. nDS1a = nDS1b = nDS1c) or 1 sample more or less than the others (e.g. nDS1a ± 1).

Finally, I recombined each smaller data subset to obtain new three equal-stratified data subsets, i.e. Fold1 (nDS1a + nDS2a + nDS3a), Fold2 (nDS1b + nDS2b + nDS3b) and Fold3 (nDS1c + nDS2c + nDS3c), which were later used to calibrate and validate the models. The schematic diagram showing this resampling procedure is depicted in Fig. 2.3B. To be consistent, I also made three replicates of the 3 fold CV for this resampling approach, by reshuffling the original data and then following the same procedure. The calibration and validation of the models was performed separately for each resampling approach.

For each scenario (hindcasting and forecasting), nine models (3 models of 3 fold CV × 3 replicates) were built for each resampling approach in each period. In this way, a performance value of the nine different models was calculated. From the field data, I built 36 models ([3 models of 3 fold CV hindcasting × 3 replicates × 2 resampling types] + [3 models of 3 fold CV forecasting × 3 replicates × 2 sampling types]) for each period. In the same way, 36 models were built based on the combined data for each period. Therefore, I totally built 144 models over the four periods for each data type (field and combined data). Kappa and CCI were used to evaluate the model performance.

Main Results

Macroinvertebrate communities and diversity patterns in the LMB

Overall macroinvertebrate communities

In total, 299 taxa, 196 genera, 90 families, 23 orders and 5 clades, accounting for a total of 21,810 individuals (inds), were identified; of which, 131 taxa (44%) and 6,481 inds (30%) were insects, 98 taxa (33%) and 10,603 inds (49%) were molluscs, 38 taxa (13%) and 2,054 inds (9%) were crustaceans, and 32 (10%) and 2,672 inds (12%) were annelids. The most commonly found insect orders were Diptera (37 taxa), Ephemeroptera (32), Odonata ( 22) and Trichoptera (20). For molluscs, most taxa belonged to the order Unionida (18), Veneroida (15) and Caenogastropoda (50); whereas for crustaceans, most taxa belonged to the order Decapoda (18)

and Amphipoda (9). Annelids were mainly represented by the order Haplotaxida (15).

At the family level, insects were mainly represented by Chironomidae (Diptera; 24 taxa) and Gomphidae (Odonata; 14 taxa). Molluscs were characterized by Unionidae (Bivalvia; 18 taxa), Corbiculidae (Bivalvia; 14 taxa), Viviparidae (Gastropoda; 12 taxa) and Stenothyridae (Gastropoda; 9 taxa). Crustaceans and annelids were represented by Palaemonidae (10 taxa)

and Naididae (15 taxa), respectively. Regardless the taxonomic groups, 36 macroinvertebrate families were represented by only one species (see Appendix T2).

Over all of 299 species reported, 20 were common (present in ≥16 samples or ≥25% occurrence), 106 were uncommon (present between 4-15 samples or between 6-24% occurrence) and 173 were rare (present in ≤3 samples or ≤5% occurrence). Alien species were not accessed due to the lack of distribution data and the limited taxonomical knowledge of each species. The most widely distributed species belonged to two insects: Ablabesmyia sp. (73% occurrence) and Polypedilum sp. (70%) and one was a mollusc, Corbicula tenuis (67%). In addition to being widely distributed, these 3 taxa were among the top 10 most abundant. Of the total individuals, Ablabesmyia sp. accounted for 2.9%, Polypedilum sp. for 3.8%, whereas the 3 most abundant species, Corbicula leviuscula, Limnoperna siamensis and Corbicula tenuis, accounted for 8.4%, 6.1% and 5.8%, respectively. The information of each species occurrence is provided in the Appendix T2.

Spatial community patterns and their relationship with environmental factors

Community clusters and indicator taxa

Based on the dissimilarity of macroinvertebrate abundance data and cluster analysis, the 63 sampling sites could be distinguished into four clusters (Fig. The clusters represent four different macroinvertebrate assemblages, which were characterized by different environmental conditions, macroinvertebrate richness, abundance and diversity (p<0.01) (Fig. 3.2, Table 3.1). An increasing richness, abundance and diversity of macroinvertebrates was found from the upstream (i.e. tributaries) to the downstream assemblage (the delta). The richness and abundance of molluscs, crustaceans and annelids generally increased from the up to the downstream assemblage, while insects dominated the upstream and tributary assemblages. 12.5 (7.91) -IIb2 24.38 (15.76) WT: water temperature, DO: dissolved oxygen, EC: electrical conductivity, SAW: the surface area of watersheds, Mix_evg.dec: mixed evergreen and deciduous forests. The number of samples [n] in each cluster is indicated between square brackets. * indicates ANOVA and Kruskal-Wallis Test for significant differences among clusters at p<0.05. Superscripts (IIa, IIb1, IIb2) indicate significant pair-wise comparisons between the corresponding cluster (each column) and superscript-labeled clusters (i.e. IIa, IIb1, IIb2) at p<0.05. The number of indicator species followed the overall trend of macroinvertebrate richness and abundance: the delta assemblage (I) were represented by 53 indicator species, most of which were molluscs, annelids and crustaceans. The upstream assemblage along the main channel were presented by 14 indicator taxa, most of which were insects. The in-between assemblage (IIb1, between the delta and main upstream assemblages) and the tributaries (IIb2) were represented by two different indicator taxa. The detailed information on indicator taxa for each assemblage is provided in Table 3.2.

Relationship between community clusters and environmental factors

The

Macroinvertebrate diversity and its relation to environmental factors

Diversity and variation in important taxa

The diversity of macroinvertebrates remained relatively high across the LMB. The mean richness and abundance at each site was 23 species (range: 6-74 species) and 346 inds (range:

13-2009 inds). Alpha (α) diversity (H') at each site was 2.1 (range: 0.8-3.3). Beta diversity of the global macroinvertebrate communities (i.e. the communities that include all component communities: molluscs, crustaceans, annelids and insects), measured as the total variance, was exceptionally high, at BDTotal =0.80 on a 0-to-1 scale. When β diversity of each component community was computed separately, the total variance of mollusc communities was the highest (BDMTotal =0.78), followed by insect (BDITotal =0.74) and annelid communities (BDATotal =0.72). Crustacean communities had the lowest total variation (BDCTotal =0.38).

A total of 60 macroinvertebrates were identified as the important species (i.e. the species that had SCBD indices larger than the mean SCBD (0.003). Among them, 29 species belonged to insects, 18 to molluscs, 7 to annelids and 6 to crustaceans (Table 3 6) and the order Veneroida (5) and Unionida (4). Crustaceans and annelids were dominated by the order Decapoda (6 taxa) and Amphipoda (4), and Rhynchobdellida (7)

and Arhynchobdellida (5), respectively.

The most occurring families of insects were Corixidae (Hemiptera, 9 taxa), Coenagrionidae (Odonata, 8), of molluscs was Planorbidae (Hygrophila, 11) and of annelids was Glossiphoniidae (Arhynchobdellida, 5). Crustaceans were identified only to family level, and that did not reveal which family dominated the communities. In total, 73 alien taxa have been collected across river system in Flanders (see Boets et al. 2016).

Spatio-temporal community composition and environmental factors

Spatio-temporal community composition

From the past to the recent period, there was a gradual increase in macroinvertebrate richness and abundance (Fig. 3.4, Table 3.4). The average richness and abundance per site for the four periods (D1, D2, D3 and D4) were: 14±0.1 taxa and 380±32 inds, 14±0.2 taxa and 279±20 inds, 16+0.2 taxa and 456±24 inds, and 16±0.2 taxa and 1009±70 inds, respectively. The descriptive information on the richness and abundance of the global macroinvertebrate communities and of each component were provided in Table 3.4. In the first two periods in the 1990s, significantly seasonal differences (i.e. Spring vs Summer vs Autumn) was observed for macroinvertebrate abundance (Friedman chi-squared F =14.8; p<0.001) in the early 1990s (D1), and for macroinvertebrate richness (F =40.6, p<0.001) and abundance (F =20.3, p<0.015) in the late 1990s (D2). During the 2000s, no significant difference was found for macroinvertebrate richness nor abundance. The richness and abundance between each season across the four periods, e.g. Spring (D1) vs Spring (D2) vs Spring (D3) vs Spring (D4), also showed significant difference at p<0.01 (Fig. 3.4).

Relationship between community composition and environmental factors

Based on the RDA analysis, the environmental variables explained a small amount of variance Code of taxa with small loadings score were removed to improve legibility; see Table 3.5 for full names of the taxa.

Spatio-temporal diversity pattern and its relation to environmental factors

Diversity and variation in important taxa

Mean α diversity (Shannon-Weiner, H') of the global communities gradually increased from the early 1990s (H'=0.6±0.02) to the early 2000s (H'=1.9±0.01), but decreased in the late 2000s (H'=1.6±0.02) (Fig. 3.6A). For each component community, α diversity of insects was the highest, followed by molluscs' and annelids', while crustaceans had the lowest diversity.

However, the trend of α diversity of each component followed the global diversity trend.

Based on the Friedman test, α diversity of the global communities among the four periods was significantly different (Friedman chi-squared F =308, p<0.001). Significant seasonal variation in α diversity (F >25.1, p<0.001) was also detected in each period of the 1990s (D1 and D2),

but not in the 2000s (Fig. 3.6B). When each season was compared across the four periods, significant difference was always observed (Fig. 3.6B). The information on α diversity of the component community is summarized on Fig. 3.6A.

From the past (D1, 1991) to the most recent (D4, 2010) period, an increasing trend in β diversity (the total variation) of global macroinvertebrate communities was observed. The mean BDTotal for the corresponding period D1, D2, D3 and D4 was 0.52±0.07, 0.50±0.03, 0.54±0.02 and 0.59±0.01 (Fig. 3.7A). Based on seasons, the total variation of the global communities was lower in spring, compared to summer in autumn, for the first two periods. The BDTotal for the first two periods in spring, summer and autumn were: 0.44, 0.54 and 0.58 (for D1), and were 0.46, 0.53 and 0.51 (for D2), respectively. For the later two periods (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), the amount of seasonal total variation in each period was relatively similar (Fig. 3.7B).

When β diversity of each component community was separately computed for each period, molluscs always had the highest variation, followed by insects and crustaceans. Annelids always had the lowest variation. Over the four periods, β diversity of molluscs and insects increased from the past to the most recent period, whereas β diversity of annelids and crustaceans fluctuated. The information on β diversity of each component community in each period is shown in Table 3.5.

Relationship between diversity and environmental factors

Multiple regression models and the AIC criterion identified five variables which most of the time showed a significant effect on α diversity. The five variables included ammonium (NH4 + ), nitrite (NO2 -), nitrate (NO3 -) and electrical conductivity (EC) (negatively associated) and dissolved oxygen (DO) (positively associated). The global LCBD indices were also mostly affected by five variables including sinuosity, NH4 + and NO2 -(negatively associated), and EC and pH (positively associated). These five variables explained a small amount of variation of global LCBD indices (Table 3.7). When the global LCBD indices were regressed against the LCBD indices of the component communities, the LCBD indices of each component always significantly influenced the variation of the global LCBD indices (Table 3.8).

Modelling and predicting

Performance variation of modelling techniques applied in the LMB

The occurrence of 199 species were predicted. The overall performance of modelling techniques used (i.e. LR, RF, SVM and ANN) was significantly different when considering the three performance measures: Kappa (Friedman chi-squared, F=12.3; p=0.006), AUC (F=11.2, p=0.01) and ER (F=350, p<0.001). The highest mean Kappa (0.19) and mean AUC (0.60) were obtained for ANN, followed by LR (mean Kappa: 0.16, mean AUC: 0.59) and RF (mean Kappa: 0.12, mean AUC: 0.55), while SVM yielded the lowest mean Kappa (0.06) and mean AUC (0.53). On the other hand, a lower mean ER (0.09) was obtained for RF and SVM, while ANN and LR had a higher mean ER (0.13 and 0.16, respectively).

Based on Kappa and AUC, the performance of the models varied for different prevalence ranges (Fig. 3.8A-B). ANN and LR performed better than RF and SVM for the prevalence range <0.1; ANN, RF and LR performed better than SVM for the prevalence range 0.1-0.2; ANN and RF performed better than LR and SVM for the prevalence range 0.2-0.3, but not significantly different according to the standard error. Based on ER, RF and SVM performed better than ANN and LR in predicting species with a prevalence range <0.1, between 0.1 and 0.2 and between 0.2 and 0.3. Based on all calculations, the model performance was not significantly different for the prevalence range ≥0.3 (Fig. 3.8A). Linear regression models showed that the prevalence always (for each modelling technique) had a positive and significant effect on the three performance measures used. When the same data were analysed using a quadratic regression model, the explained proportion of variance of Kappa and AUC, and of ER increased, and the coefficient of the quadratic term was always negative and highly significant (p<0.001, Table 3.9).

Modelling alien mollusc occurrence and their co-occurrence with native molluscs

Based on CT models built for each period, the "alien" sites were not able to be correctly predicted; only the "co-occurrence" and the "native" sites could be correctly predicted (Fig. 3.9). At sites having a low sinuosity (<1.01), the co-occurrence was dependent on chemical water quality variables (e.g. NH4 + , NO3 -, COD, pH). Sinuosity was always one of the most important variables for the models built for each period. Based on the trees, the models revealed for each period that when the sinuosity was lower than 1.01, the co-occurrence occurred where NH4 + <0.4 mg/L and NO3 -≥2.4 mg/L for D1, COD <18.9 mg/L and NO3 -in between 2.6-4.7 mg/L for D2, COD <21.8 mg/L and NO3 -in between 2.3-5.1 mg/L for D3, and where pH ≥7.2 for D4 (Fig. 3.9). Fig. 3.9 Classification trees predicting the "co-occurrence" and the "native" sites for each period. SI: Sinuosity, NH 4 + : Ammonium, NO 3 -: Nitrate, COD: Chemical Oxygen Demand, cp: pruning complex parameter.

The mean Kappa and the mean overall CCI (Overall-CCI) of the models decreased from D1 to D4 (Fig. 3.10A-B). For the four periods (D1, D2, D3, and D4), the corresponding mean Kappa was 0.34, 0.32, 0.22 and 0.16, and the corresponding mean Overall-CCI was 93%, 92%, 79% and 54%. The mean CCI of models predicting the "co-occurrence" sites (CCI-co-occurrence) was lower for the first 3 periods (28%, 25% and 23%), while it was higher for the most recent period (52%, Fig. 3.10C). On the contrary, the mean CCI of models predicting the "native" sites (CCI-Native) for the past 3 periods (99%, 98%, and 96%) was significantly higher compared to that of the most recent period (66%, Fig. 3.10D).

Optimizing the prediction of alien mollusc occurrence

The CT models optimized by incorporating the field and cloned data (i.e. a dataset obtained by independently duplicating the field data points by k different individuals) were able to hindcast and forecast the class "alien" and the other two classes. The overall performance of the optimized models decreased from D4 to D1 for the hindcasting scenario, while it decreased from period D1 to D4 for the forecasting scenario (Fig. 3.11A-B). The calibrated and validated models using the two resampling approaches also yielded a similar trend in predictive performance and in stability for each scenario (see Article 5). The CCI of the class "alien" followed the decreasing trend of the overall performance for each scenario (Fig. 3.12).

General Discussion

Overall community composition and diversity in the two river systems

The "diversity" of macroinvertebrates in the LMB were found to be greater than those in Flemish rivers. This could be explained by the fact that most invertebrates from the LMB were identified to the species level, leading to a higher diversity; while the invertebrates from Flemish rivers were identified only to family or genus level. However, these findings indeed may not reveal the true composition and diversity for the LMB because it is very large, compared to Flemish rivers, but has been little studied. Thus, it is mostly likely that the number of reported taxa from the LMB is underestimated; the further it is investigated, the higher number of taxa and diversity will be. On the other hand, Flemish rivers have been extensively investigated and monitored regularly.

Regardless of the surface areas, tropical Asia has a unique geological arrangement (Sodhi et al. 2004). The evolutionary history and together with ecological interactions can explain the high biodiversity in tropical Asia [START_REF] Agrawal | Evidence for adaptive radiation from a phylogenetic study of plant defenses[END_REF]). This region, as well as other tropical regions, is regarded as the 'cradle' and the 'museum' of species diversity because it is characterized by a high speciation and low extinction rates, compared the temperate regions [START_REF] Mckenna | Tropical forests are both evolutionary cradles and museums of leaf beetle diversity[END_REF]. Clear evidences have been shown by global diversity of many taxonomic groups including stoneflies, dragonflies, mussels, crabs and shrimps (De Grave et al. 2008;[START_REF] Fochetti | Global diversity of stoneflies (Plecoptera; Insecta) in freshwater[END_REF][START_REF] Kalkman | Global diversity of dragonflies (Odonata) in freshwater[END_REF][START_REF] Yeo | Global diversity of crabs (Crustacea: Decapoda: Brachyura) in freshwater[END_REF]Bogan 2008); a higher number of taxa (species, genera, and families) are mostly found in tropical Asia than temperate and some other tropical regions.

Based on the computation of β diversity for each component community, molluscs and insects had a higher total variation for both systems, compared to crustaceans and annelids. This could be due to the responses of different taxonomic group of each components to environmental changes and to the ecological process of each system that can support different taxonomic groups. For example, mayflies and stoneflies (i.e. genera) have been found to have a high similarity between temperate Europe and tropical Asia, compared to some other geographic regions [START_REF] Barber-James | Global diversity of mayflies (Ephemeroptera, Insecta) in freshwater[END_REF][START_REF] Moor | Global diversity of caddisflies (Trichoptera: Insecta) in freshwater[END_REF]. A high gastropod species diversity have been reported from Europe (e.g. Southern France and Spain, Southern Alps and Balkans regions) and tropical Asia (the LMB) [START_REF] Strong | Global diversity of gastropods (Gastropoda; Mollusca) in freshwater[END_REF]. This may indicate high β diversity of insects and molluscs is favoured by ecological and environmental conditions from both regions which promote their distribution and abundance, and thus lead to a higher β diversity. A more detailed discussion on the differences and on responses to environmental variables is provided separately in the following paragraphs for each system.

Spatio-temporal changes of communities and their relation to environmental factors

The LMB

Spatial organization and environmental factors

The compositional and diversity differences found for the four macroinvertebrate assemblages were related to spatial environmental heterogeneity observed from the upstream to the downstream parts of the LMB (Dobson et al. 2002;Heino et al. 2005;Wang et al. 2012;Salmah et al. 2014). This can be explained by the proportion of the richness, abundance, diversity and the represented indicator taxa of the different taxonomic groups for each assemblage. Although assemblage IIa, IIb1 and IIb2 were characterized by a low and similar macroinvertebrate richness, abundance and diversity, each component (i.e. annelids, crustaceans, molluscs and insects) contributed differently to the overall macroinvertebrate community composition for the particular assemblages. On the other hand, assemblage I had a high richness, abundance, diversity and a high number of indicator taxa. Therefore, these findings clearly suggested that the community composition of the four assemblages, particularly assemblage I, is distinctly organized and structured along the environmental gradients of the LMB.

Site-specific characteristics, connection in habitats and mesohabitats and variation in river morphology can strongly affect community structure and composition found from the upstream to the downstream parts of rivers (Thorp and Delong 1994;Zilli and Marchese 2011;Mazão and Bispo 2016). Among the measured environmental variables, physical conditions and land cover types are the key factors determining macroinvertebrate assemblages along the LMB.

The surface area of watershed, river depth and altitude are the most important variables discriminating the upstream assemblages (IIa and IIb2) from the downstream assemblage (I).

This discrimination is clearly explained by the axis 1 of the LDA accounting for 65.7 % of the total variance of the data and by a high predictive reliability for each particular assemblage (Fig.

3.3A-B

). Because of the observed physical conditions, together with other related variables (e.g. water temperature, dissolved oxygen and land cover types), the assemblages were characterized by a very different community structure, composition and diversity.

The intermediate assemblage (IIb1) in-between was positively linked to three types of land cover: agricultural land, inundated and wetland. Along axis 3 of the LDA ordination, these land cover types are the key factors best explaining the macroinvertebrate composition of assemblage IIb1 (Fig. 3.3C-D). However, the proportion of inundated and wetland constituted a small percentage of the land cover, and together with agricultural land, they only explained a limited amount of the variance of the data for axis 3 of the LDA model (11.2%). As a result, the predictive performance for this assemblage was not as high as the performance for the other assemblages (see Sor et al., 2017a). A broad range of values of physical-chemical variables (e.g. dissolved oxygen, river width and depth, the surface area of watershed and altitude) found at sites belonging to this assemblage could also explain the performance yielded. On the other hand, assemblage IIa was associated with two important factors (electrical conductivity and urban areas) according to the LDA axis 2, which explained a higher variance of the data (23.1%). Therefore, these two factors, together with other measured environmental variables, could correctly predict this assemblage with a higher reliability. The possible explanation for the association found is that many sites that belong to this assemblage are located along the tributaries (e.g. the Mun and Chi river basins in Thailand) and the main channel where cities were built. The tributaries are surrounded by intensified agriculture and the cities are exposed to a high level of anthropogenic disturbance (Dao et al. 2010;Kudthalang and Thanee 2010).

Therefore, the runoff from the surrounding agricultural areas and the discharge of urban wastewaters may cause the increase in electrical conductivity (Wetzel 2001).

Spatial diversity pattern and environmental factors

The distinct organization of community composition and α diversity from the up to the downstream leads to a high β diversity. Alpha diversity pattern follows the richness and abundance patterns of each assemblage; a high diversity was observed at the downstream sites, which are characterized by a large watershed surface area, deep and wide rivers and a high water temperature, and a low diversity found for the upstream and tributary sites, which are characterized by highland area, clear water, high DO concentration and dense forest cover.

The most important taxa, which mainly contributed to the high β diversity, found at the downstream were dominated by annelids, crustaceans and molluscs. This is because they are more abundant in the downstream part of the LMB, as in the case of two mollusc species:

Corbicula leviuscula and Limnoperna siamensis. Unsurprisingly, this result supports previous studies (Arscott et al. 2005;Collier and Lill 2008;Wang et al. 2012). On the other hand, at the upstream sites, most of the important taxa belong to insects. This could be due to that fact that clear water and high values of dissolved oxygen, which are mostly found in tributaries and upstream sites, are mainly preferred by insect taxa (Dobson et al. 2002;Collier and Lill 2008;Królak and Korycińska 2008). For instance, the two important taxa found with the highest SCBD indices, Polypedilum sp. (I033) and Ablabesmyia sp. (I017) were highly associated with sites having high values of SD, which were mostly in tributaries (see Sor et al., 2017a).

The β diversity of macroinvertebrates was estimated as the total variance (BDTotal) of the communities found at the sampling sites and computed the contributions of individual sampling sites (LCBD indices) to BDTotal. LCBD values are scaled to add up to 1 over the whole study; the mean LCBD value in the LMB was thus 1/63 = 0.016. The variation in LCBD indices (range: 0.010-0.023) appeared to be related to the environmental factors that mainly explained the community composition and α diversity. LCBD indices were positively associated with river depth, surface area of watersheds and electrical conductivity. High values of these variables are most of the characteristics of sites situated in the delta and along the main channel of the LMB (Sor et al. 2017a). This indicates that a high degree of uniqueness in composition of benthic macroinvertebrates in this tropical system mostly occurred along the mainstream of the rivers, which are highly associated with anthropic disturbance (see Article 2). Small LCBD indices were linked to the tributary sites, which are characterized by a high values of Secchi depth (rivers with clear water). Due to this fact, tributaries mainly support particular assemblages of particular macroinvertebrate taxa, e.g. Ephemeroptera, Plecoptera, Trichoptera and Diptera (Dobson et al. 2002;Wang et al. 2012).

When regressed with LCBD indices of component communities, the global LCBD indices was highly linked to the taxonomic compositions of the mollusc and insect communities (LCBDM + LCBDI indices). However, the degree of uniqueness in composition of macroinvertebrate communities (global LCBD indices) is expected to be contributed by the component communities. In the LMB, LCBDs of mollusc and insect communities, which had a higher total variation (BDMTotal = 0.78 and BDITotal = 0.74, respectively), explained most of the global LCBD variation because these two groups had higher abundances and wider distributions than the annelid and crustacean communities, which had lower total variation (BDATotal = 0.72 and BDCTotal = 0.38, respectively). De'ath (2002) and Davidson et al. (2010) explained that taxa with low richness and low occurrence lead to a less variance of the community composition, and this is similar to the present findings for the annelid and crustacean communities.

Flemish rivers

Spatio-temporal changes in composition and environmental factors

The benthic community composition in Flemish river system was found to gradually change (e.g. macroinvertebrate richness and abundance) from the early 1990s to the late 2000s. The changes can be seen from the RDA plots (Fig. 4.1A-B). The circular arrangement of the monitoring periods is due to the gradual changes in species composition, and thus some taxa highly contributed to the variation in composition in the past periods, but not in the recent periods, and vice-versa (Table 3.6, Fig. 4.1B). However, the seasonal changes in recent periods (2000s) were not significantly different, except for the past periods (1990s). This could also mean that the water quality condition across seasons in each period of 2000s is more or less the same, and that allows similar composition of macroinvertebrates to occupy the habitats.

The temporal changes over the four periods observed corresponded to the improvement of water quality condition from the past to the recent period (Leuven et al. 2009). Clearly, the centroids of community composition in the 1990s were linked to a high value of ammonium, nitrite, total phosphorus, and COD, which indicate a poor water quality. Whereas the centroid of community composition in the early 2000s was more related to a high value of DO, and thus the communities was composed of many moderate-to-high pollution sensitive insect taxa, e.g. Simuliidae (I030), Hydropsychidae (I119), Baetis (I037) (Gabriels et al. 2010). In the late 2000s, the community centroid was more related to a high value of pH and EC. This could be due to the fact that many samples from this period were mainly recorded from brackish polder watercourses and from the main harbour watercourses, e.g. river basin Ijzer, Brugse Polders (Bru-Pol), Gentse Kanalen (Gen-Kan) and Beneden-Schelde (Ben-Sch) (Fig. 4.1B) where a high level of seawater intrusion occurs and a high intensity of human-related activities takes place. Brugse Polders, Gen-Kan: Gentse Kanalen, Ben-Sch: Beneden-Schelde, Bov-Sch: Boven-Schelde, and Dij-Zen: Dijle Zenne. Taxa (red colour) that had small loading scores on the RDA axes were removed to improve visibility; see Table 3.5 for full name of the taxa.

Spatio-temporal diversity patterns and environmental factors

Alpha diversity of macroinvertebrates was undoubtedly related to the improvement of water quality in Flanders. A gradual increase in α diversity from the early 1990s to the early 2000s could be explained by the gradual increase in the overall number of taxa and abundance of all macroinvertebrates and of each component (Fig. 3.4, Table 3.4). However, the number of taxa of the global communities in the late 2000s slightly decreased or remained more or less stable, compared to the early 2000s, but the their abundance spectacularly increased. This pattern was also observed for the four component communities, of which insects and molluscs had a higher α diversity, compared to annelids and crustaceans (Fig. 3.6). This could be due to the dominance of some taxa (e.g. alien taxa), which highly benefit from the water quality improvement, and thus increases the homogenization of the communities and subsequently reduces the diversity [START_REF] Rahel | Homogenization of freshwater faunas[END_REF][START_REF] Rahel | Biogeographic barriers, connectivity and homogenization of freshwater faunas: tt's a small world after all[END_REF]. A clear evidence can be observed from the dominance of an alien mollusc Physella over the native molluscs Physa in Flemish rivers in the recent periods (see Appendix T4). The seasonal differences in α diversity were found only for the 1990s (significantly low in Spring), but not the recent periods (2000s), as found for richness and abundance explained earlier.

Similar to other studies, α diversity of the animals in Flemish rivers is affected by nitrate content (e.g. ammonium, nitrite, nitrate), electrical conductivity and dissolved oxygen (Mereta et al. 2012;Md Rawi et al. 2013;Suhaila and Che Salmah 2014). The higher concentration of nitrate content and electrical conductivity, the more it influences the living conditions of the benthic communities (Friberg et al. 2010;Boets et al. 2013), and thus negatively affects α diversity, as found in the present study (Table 3.7). A positive association between α diversity and dissolved oxygen was also observed, and this is generally the case found in previous investigations [START_REF] Baptista | Diversity and habitat preference of aquatic insects along the longitudinal gradient of the Macaé river basin, Rio de Janeiro, Brazil[END_REF]Md Rawi et al. 2013;Suhaila and Che Salmah 2014).

Beta diversity of the global communities across sampling sites in Flemish rivers was relatively high (range: 0.49-0.60). A higher β diversity was found for the recent periods, compared to the past. The same pattern was found for β diversity of each component (Table 3.5). This indicates that the community composition varies a lot from the past periods (1990s), when the water quality was very poor, to the recent periods (2000s), when the water quality was substantially improved. The water quality improvement also influenced the seasonal β diversity of global communities in each period. Similar to α diversity patterns, a lower β diversity was found in Spring, compared to Summer and Autumn, during the 1990s, but not during the 2000s.

The five most important taxa (i.e. Naididae (A11), Gammaridae (C07), Asellidae (C15), Chironomidae non thummi-plumosus (I016) and Chironomidae thummi-plumosus (I017)) were very widespread and abundant (see Appendix T4) over the four periods. This is the reason that these five taxa greatly contributed (between 5%-14%) to β diversity in each period. However, Naididae (A11), Asellidae (C15) and Chironomidae thummi-plumosus (I017) were closely linked to the centroid of site distribution in the early and late 1990s (Fig. 4.1B), suggesting that they are able to withstand the poor water quality (Feld and Hering 2007;Gabriels et al. 2010).

The other taxa that were responsible for the temporal changes in community composition include Enchytraeidae (A08), Haplotaxidae (A09) and other taxa (e.g. Bythinella, Blephariceridae and Arctocorisa), which occurred in the 1990s, but disappeared in the 2000s (Table 3.6, Appendix T4). Moreover, the different composition of communities is also caused by a wide distribution and high abundance of Dryopidae, Gerris, Calopteryx, Physella, and the present of new taxa (e.g. Menetus (an alien mollusc), Panopeidae (crustacean), Somatochlora and Procloeon (insect)) in the 2000s (Table 3.6, Appendix T4).

The β diversity among sites for each period was estimated from the LCBD indices of corresponding period. In most cases, the LCBD values were positively associated with pH and EC (as found in the LMB), while negatively associated with nitrate content (ammonium and nitrite) and sinuosity. The association found suggests that sites with a high degree of uniqueness (i.e. high LCBD values), which may indicate high or low species richness as shown in Legendre and De Cáceres (2013), mainly occur in straight rivers that have low concentration of nitrogen content, but high concentration of EC and pH. This seems to be in the case of sites situated in the main harbour watercourses and in brackish polder where rivers have a low sinuosity, but have high values of conductivity and pH (Fig. 4.1B). The physical-chemical conditions of these watercourses may reduce taxonomic composition. Only some taxa belonging to Oligochaeta, Chironomidae (Diptera) and Physidae (Molluscs), which are pollution tolerant and can occur in areas with high pH and EC values (Rodrigues Capítulo et al. 2001;Feld and Hering 2007;De Troyer et al. 2016), remain, and thus taxonomic composition may become unique (high LCBD indices), but not more diverse, by having their number of taxa reduced.

The LCBD indices of global communities were found to be highly contributed by the four component communities' LCBDs, as reflected by their significant correlation (Table 3.8). This result is somewhat different from that found for the LMB; only molluscs and insects that were highly and significantly contributed to the global LCBD values. However, the finding for Flemish rivers is more logic because, as mentioned earlier, the global LCBD indices is expected to be contributed by the LCBDs of each component community.

Model development, performance and predictions

Modelling techniques and their application in the LMB

According to Kappa and AUC, which have been shown to be a better measure compared to error rate that is highly biased (Fielding and Bell 1997;Manel et al. 1999;Sor et al. 2017b), ANN performed the best across the complete prevalence range and a lower prevalence range (i.e. <0.1, in the case of rare species), compared to other techniques (Fig. 3.8). Some authors (e.g. Mastrorillo et al., 1997;Pearson et al., 2002;Segurado et al., 2004) agree that ANN provides advantages over other techniques for predicting species occurrence, whilst others found RF to be better (Grenouillet et al. 2011;Gallardo and Aldridge 2013). This could be due to different input predictors, parameterization settings and the applicability of each technique.

Setting the complex parameter of the models has to be carefully taken into account because they may lead to an over-fitted model, which produces a less general result with a lower applicability to different situations (Babyak 2004). Model applicability is dependent on the type of data (e.g. missing values and data distribution, De'ath and Fabricius, 2000;Therneau and Atkinson, 1997) and on where the data is derived (e.g. different ecological regions, Guisan and Thuiller, 2005;Randin et al., 2006). Most of the studies evaluate and apply modelling techniques based on the data collected within particular regions, but the applicability of those techniques to a different geographic range is hardly assessed (Fielding and Haworth 1995;Özesmi and Mitsch 1997;Kleyer 2002;[START_REF] Everaert | Comparison of the abiotic preferences of macroinvertebrates in tropical river basins[END_REF][START_REF] Forio | Analysing the effects of water quality on the occurrence of freshwater macroinvertebrate taxa among tropical river basins from different continents[END_REF]). As such, the varied performance of modelling techniques may reveal a unique behaviour of each technique and its suitable applicability (e.g. ANN) for the Lower Mekong Basin and the neighbouring areas rather than for other geographic areas.

The performance of modelling techniques was highly dependent on species prevalence (Table 3.9), suggesting that each technique may have different behaviour when predicting species with different prevalence range. For all models, the best performance was found for common species with an intermediate prevalence Fig. 3.8). In between this prevalence range, there is a smaller bias for the models to select presence/absence data for training and validation sets.

A more or less equal distribution between the presence and absence, which are likely to result in many correctly predicted instances of both the true positive and true negative fractions, could be the main reason responsible for the best predictive performance (Manel et al. 2001;Allouche et al. 2006). For very common/generalist (prevalence ≥0.6) and for rare species (prevalence <0.1), the models yielded a low performance. In these cases, there is a high imbalance between the presence and absence data. Therefore, the models can correctly predict many true positive instances and a few or perhaps no true negative instances for common species, and vice versa for rare species. This explains the low performance of the models because the measures (i.e. Kappa and AUC) are designed to reflect model performance in absence and presence instances simultaneously (Cohen 1960;Cicchetti and Feinstein 1990;Zweig and Campbell 1993). Thus, a few or no instances of either the true positive or the true negative fraction results in a low performance.

Predicting alien species occurrence and their co-occurrence with native molluscs

Based on a set of chemical water quality variables, the CT models were able to correctly predict the co-occurrence alien and native molluscs, and the occurrence of native molluscs; but not the occurrence of the alien alone. Sinuosity was always one of the most important determining factors because it formed the main root in all models. Sites having a low sinuosity (<1.01), which mainly corresponds to straight rivers, may be subjected to a high number of passing ships, which is considered the main pathway of invasions (Bij de Vaate et al. 2002;Nunes et al. 2015). Moreover, straight rivers shorten travelling distances, resulting in more frequent transportations and thus allow a large and frequent amount of ballast water being released.

Consequently, due to a higher number of introductions, the survival rate of alien molluscs increases (Gollasch 2006). The hotspots of mollusc invasion (the alien and the co-occurrence sites) were mainly situated in brackish polder watercourses and at large rivers that have a shorter distance to the ports in the Rhine delta, to the coast of the North Sea and to the other large rivers (e.g. Meuse River). This observation is also supported by evidence in Boets et al. (2016) and in Grabowski et al. (2009) who found that alien fauna mostly inhabited large rivers where intensive navigation takes place. The followed key variables responsible for the co-occurrence included NH4 + , NO3 -, COD and pH. These variables are the major factors influencing the distribution and diversity of freshwater fauna (Leuven et al. 2009;Wang et al. 2012). Thus, for each period, where sampling sites had a low sinuosity (<1.01), the co-occurrence was mainly dependent on the chemical water quality status (see Article 4).

The overall performance of the models (Kappa and Overall-CCI) in predicting the cooccurrence was moderate to good. When using data from the past period, the performance was higher than when using data from the most recent period (Fig. 3.10A-B). This fact could be attributed to less complex biotic interactions and to the limited number of occurrences of alien species in the past. For example, most alien and native species had a low occurrence and abundance in the past periods (Appendix T4). In this context, alien species may invade those sites where competition is low and where a few native species occur. Some alien taxa have a strong ability to compete for food and niche with the native species. For examples, in the past, the alien Physella occurred at only few sites with a few individuals, but then widely spread with high abundance in the recent periods, leading to a reduced abundance of the native Physa (Appendix T4). This could be a result of competition for food and niche with the alien Physella. Araújo and Luoto (2007) and Meier et al. (2010) stated that the biotic interaction is important in predicting species distribution, and in general when included in the model, it increases the predictive performance. Provided that information on biotic interaction was not included in the models and that the alien molluscs already occupy a wide range of environmental conditions in the most recent data, the overall performance of the models based on the most recent period is somewhat lower. Moreover, during the continuing expansion phase of invasion (the late 2000s), the range of environmental conditions where alien species occur increased. At the same time, many native species recovered their occurrences and abundance, due to improved water quality.

These findings might also explain the higher co-occurrence and thus the lower overall performance of the models.

However, the prediction of the "co-occurrence" of alien and native molluscs was more reliable for the recent period compared to the first three periods. This could be explained by the increased co-occurrence sites observed in the most recent period (see Article 4). This is similar to the models predicting the occurrence of native molluscs. The predictive models yielded a high performance when the prediction was based on the periods (i.e. the first three periods) that have a large sample size of the "native" sites. This is quite logic as for predictive models, the more input samples provided, the better the models learn and as a result, a higher predictive performance can be obtained (Stockwell and Peterson 2002;Hernandez et al. 2006).

Optimizing the prediction of alien mollusc occurrence

Since the CT models were not able to predict the sole occurrence of alien mollusc, which could be due to a low number of instances of this particular class, the CT models were optimized via a hindcasting and forecasting approach. Based on the best CT model configuration (for both in terms of performance and interpretation) obtained, a similar sample size of each class appears to be the best input, and thus should be considered in every predictive model as has been generally suggested (e.g. Manel et al. 2001;Allouche et al. 2006). A clear evidence was revealed from the CT models that used the field data (field observations). For these models, the prevalence of each predicted class was not equally distributed and the models were unable to predict the putative occurrence of the alien molluscs. However, by using the combined data, their occurrence could be correctly predicted at past periods (D1, D2, and D3) when the CT models were calibrated based on the data of the most recent period (D4) and at the recent periods (D2, D3 and D4) when the models were calibrated based on the data of D1. These results correspond to the field past and current situation.

One of the challenges in predicting a species' occurrence is obtaining a balanced class distribution for the response variable, because a small sample size of each predicted class can be the source of instability and errors in species distribution models (McPherson et al. 2004;Allouche et al. 2006). However, results from this case study demonstrated that having the exact same class distribution (i.e. stratified split) is not always necessary during the calibration and validation process because the models can make the correct prediction and yield a similar performance when the prevalence reaches a certain threshold. However, the model is incapable of predicting a class with a very low prevalence, as was the case of the class "alien" in this study. As collecting new field data from past/recent periods is unfeasible/costly, optimization appears necessary (Hirzel and Guisan 2002). In respond to this, I optimized the prediction by considering cloned data for the models, resulting in a better prediction and higher reliability.

Cloned data have been applied and recently suggested for hierarchical models in ecology (Lele et al. 2007;Ponciano et al. 2009;Lele et al. 2010). The increasing number of clones used in the models can increase the maximum likelihood of the prior class distribution, but it does not affect the statistical accuracy, which mainly depends on the information of the field data and the model calibration and validation process (Lele et al. 2007). This suggests that CT models can be successfully optimized and improved to predict the actual occurrence of alien molluscs by incorporating cloned data into the models.

General Conclusion and Perspectives

General conclusion

The different macroinvertebrate community composition, which leads to a different amount of variation (i.e. β diversity), is mostly likely driven by environmental conditions characterizing each system. A higher diversity is found for the LMB, compared to Flemish rivers.

Identification effort could be a reason causing the differences; most benthic animals from the LMB were identified to the species level, whilst those from Flemish rivers were identified only to family or genus level. However, these results tend to be revealing only a portion of knowledge on macroinvertebrates for the LMB, because this basin is very large, but has been seldom studied. Thus, the number of reported taxa is most likely to be underestimated; the further investigations will greatly increase the systematic and ecological knowledge on macroinvertebrate from this basin. On the other hand, Flemish rivers have been extensively investigated and monitored regularly over the past decades. Nevertheless, the two systems showed some similarities; molluscs and insects had a higher total variation and highly contributed to the total variation of the global communities, compared to crustaceans and annelids. Moreover, a high degree of uniqueness in community composition (i.e. high LCBD values) frequently occurred at sites where a higher level of anthropic disturbance was observed.

In the LMB, the most important variables influencing the macroinvertebrate composition and diversity include altitude, surface area of watersheds, river width and depth, Secchi depth, DO, EC and water temperature. Together with the highly variable topography, geology, hydrology and different land cover, the composition and diversity of macroinvertebrates in the LMB is spatially and distinctly organized from the up to the downstream of the basin.

In Flanders, river morphology (e.g. sinuosity) and physical-chemical water quality variables (e.g. DO, pH, NH4 + , NO3 -, NO2 -, COD and EC) are the key variables. Water quality improvement from the past to the recent period had gradually enhanced the composition and diversity of benthic animals, leading to an increased β diversity in the recent period. The improvement of water quality also promoted macroinvertebrate composition and diversity, which significantly varied among seasons (i.e. Spring vs Summer vs Autumn) in the past 1990s, to be seasonally similar in the recent periods (2000s). However, good water quality appeared to drive community homogenization. This is clearly indicated by a decreased local (α) diversity in the 2000s.

Among modelling techniques applied in the LMB, ANN performed the best across the complete species prevalence range. ANN also yielded the highest performance when predicting the occurrence of rare species. In Flemish rivers, the CT models could only predict the cooccurrence of alien and native molluscs, and the occurrence of native molluscs; but not the sole occurrence of alien species. Via a hindcast-and forecasting approach, the CT models were successfully optimized to reliably forecast and hindcast the sole occurrence of alien molluscs, by incorporating cloned observations into the models. Whether in the past or recent periods, the results of this optimization correspond to field observations.

Implications for management and restoration

Diversity of global and component communities is biologically and ecologically important for ecosystem structure and functioning. Thus, a sudden change in the composition of any taxonomic group may result in disproportionate or unexpected responses of other taxa in the system (Naeem 1998), and consequently alter ecological processes (Covich et al. 1999). In this context, the indicator species, the species/taxa with high among-sites variance (high SCBD indices) and the uniqueness in community composition at sampling sites (LCBD indices) identified/analysed in this study have provided insights into the ecological importance and environmental degradations of sampled sites/river reaches (see Sor et al. 2017a, Article 2). With these results, together with the key correlated environmental factors and the implemented modelling techniques, habitat quality and suitability for the most vulnerable native and pollution sensitive taxa could be more conveniently monitored and analysed. This knowledge could provide an advantage to support decision-making concerning conservation, management and restoration planning of these keystone taxa and their communities at local and regional scales.

With the challenge of increasing rate of invasion, which has become a major concern for global economy and environment (Sala et al. 2000), numerous threats have put high pressure on native biodiversity. Although some alien species (e.g. Dreissena mussels and Corbicula Asian clams, the well-known exotic species for European and American river systems), also play a key role in ecological functioning, e.g. improving water quality via filtering process [START_REF] Higgins | The effect of dreissenid invasions on chlorophyll and the chlorophyll: total phosphorus ratio in north-temperate lakes[END_REF][START_REF] Sousa | Invasive bivalves in fresh waters: impacts from individuals to ecosystems and possible control strategies[END_REF], their wide spread and high abundance substantially affect other species from the local to ecosystem levels. In the case of Dreissena and Corbicula, they compete with native species, reduce plankton communities, cause a decline of dissolved oxygen and transfer organic resource from the water column to the sediments [START_REF] Caraco | Dissolved oxygen declines in the Hudson River associated with the invasion of the zebra mussel (Dreissena polymorpha)[END_REF][START_REF] Descy | Modelling the impact of benthic filter-feeders on the composition and biomass of river plankton[END_REF][START_REF] Caraco | Top down control from the bottom: Regulation of eutrophication in a large river by benthic grazing[END_REF][START_REF] Vaughn | Unionid mussels influence macroinvertebrate assemblage structure in streams[END_REF]. Moreover, once they have successfully colonized new habitats, eradication is hardly possible (Regan et al. 2006). This consequently leads to native species replacement, food web reorganization and community composition simplification (Gurevitch and Padilla 2004;Bernauer and Jansen 2006;Didham et al. 2007). Due to these facts, applicable risk assessment approaches are required to mitigate ecological and socioeconomic impacts posted by alien species, and to develop sound prevention and management options [START_REF] Panov | Assessing the risks of aquatic species invasions via European inland waterways: from concepts to environmental indicators[END_REF]. In this study, an assessment (i.e. the case in Flemish rivers) using predictive models including a hindcasting and forecasting approach provides more insights into the ecology of alien species (e.g. the environmental conditions where alien species solely exist and co-exist with native species), which can serve as a basis for invasion control and be used to support environmental management and conservation planning. However, shall be this approach or other novel assessment tools implemented, both negative and positive effects should be taken into account, because as mention earlier, alien species can also be a good ecosystem engineers or a food sources for other species (Palmer et al. 1997;[START_REF] Chowdhury | Ecosystem engineering by mussels supports biodiversity and water clarity in a heavily polluted lake in Dhaka, Bangladesh[END_REF].

Perspectives

Information on and insights into benthic macroinvertebrate ecology and its implications for management and restoration planning provided by this study are of great significance contributing to increase scientific knowledge of the two basins, particularly the LMB which has been seldom studied. However, there is much more to be investigated. For both systems, a deeper investigation of each component community should be conducted as they differently respond to environmental conditions. It is also recommended for further rigorous examination of temporal changes at some specific sites exposed to different anthropic disturbance or have different land cover characteristics. Between these sites, temporal species composition, variations, taxonomic and functional diversity, biological traits and the degree of uniqueness in community composition (i.e. LCBD indices) should be analysed rigorously. These can provide useful and fundamental knowledge on how the communities will in the future change under different environmental conditions (e.g. from a less to a substantial disturbed sites/locations) (Legendre and Salvat 2015), and thus a management or restoration planning can be foreseen.

Moreover, should there be new data collection, other important variables including nutrients, sediment loads and habitat variables (bed rock, mud, leave litters etc.) and climatic variables (for large spatial scale like the LMB) have to be taken into account because they greatly influence macroinvertebrate community composition in both the tropical and temperate systems (Nicola et al. 2010;Cai et al. 2012;[START_REF] Pearson | Dynamics of invertebrate diversity in a tropical stream[END_REF].

Regarding the predictive modelling approach, an ensemble modelling evaluated with an array of performance measures (e.g. Kappa statistics, true skill statistics, ROC curve, error rate or the correctly classified instances) is recommended when comparing the performance of different techniques or when forecast or hindcast species distribution (Araújo and New 2007;Grenouillet et al. 2011). Moreover, incorporating the missed variables including those mentioned in an earlier paragraph, biotic interactions and dispersal vectors (as in the case of alien species) into the models is also suggested. This incorporation will improve the predictive performance and reliability of the models (Araújo and Luoto 2007;Boets et al. 2014;Parravicini et al. 2015).

Lastly, testing model transferability into a new geographic region (e.g. outside the region where the models were calibrated) is vital to validate the model applicability. This is suggested because most of studies have evaluated and applied modelling techniques based on the data collected within particular regions, but have not access their applicability outside the studied areas. This case also applies to the present study; the CT models were successful optimized to hindcast and forecast the occurrence of alien species occurrence only in Flemish river systems.

Therefore, these optimized models should be in the future validated using data collected outside Flanders, e.g. river systems in Netherland or neighbouring regions.

Appendices Appendix T1

List of sampling sites, corresponding codes used in this study (Site code) and in the Mekong River Commission report (MRC site code), and the channel and the year sampled. * indicates the site which was sampled further away from its original coordinates and which was considered as a different site in our study.

Appendix T2

List of taxa recorded in the present study and their occurrences. The numbers in the parentheses tell the number of taxa reported for each family/order/clade. The families/orders/clades that are not followed by any number mean that they are represented by only one observed taxon.

Appendix T4

List of taxa and their occurrences recorded in Flemish rivers from 1990 to 2010. The numbers of the occurrence and abundance of each taxon were calculated based on the median values of each period. D1: 1991-1995, D2: 1996-2000, D3: 2001-2005, and D4: 2006-2010. Note: Mollusc taxa in bold are alien genera which were used in the predictions in the case study of Flemish rivers. In this study, we investigated patterns of spatial variation in macroinvertebrate assemblages in the Lower Mekong Basin (LMB) and examined their relationship with environmental factors. Cluster analysis was used to group macroinvertebrate samples and Linear Discriminant Analysis was performed to discriminate the major factors associated with the macroinvertebrate assemblages. Four clusters could be distinguished based on the dissimilarity between macroinvertebrate assemblages. The assemblages related to the tributaries and the upstream parts (cluster II) were characterized by a lower richness, abundance, diversity and a lower number of indicator taxa compared to the assemblage found downstream in the Mekong delta (cluster I). Aquatic insects and their indicator taxa (e.g. Caenodes sp., Dipseudopsis sp. and Gomphidae sp.), preferring a high-altitude environment with a high dissolved oxygen concentration and a high density of wood/shrub and evergreen forests, were the most predominant group in the assemblages occupying the tributaries and the upstream parts (cluster IIa). The assemblage found in the delta, consisting largely of molluscs and a moderate richness and abundance of worms, crustaceans and dipteran insects, was mainly represented by Corbicula leviuscula and C. moreletiana (molluscs), Namalycastis longicirris and Chaetogaster langi (worms), Corophium minutum and Grandidierella lignorum (crustaceans), and Cricotopus sp. and Clinotanypus sp. (dipteran insects). This assemblage was associated with a large watershed surface area, deep and wide rivers and a high water temperature. The intermediate assemblage (cluster IIb1) in-between could be discriminated based on land cover types including inundated, wetland and agricultural land, and was represented most by molluscs. Strikingly, the assemblage occupying the upstream parts (cluster IIa), which is related to intensified agriculture and a moderate conductivity, was characterized by a higher macroinvertebrate diversity compared to the mountainous and less impacted tributaries. This could mean that the natural stress is high in these systems for some taxa, leading to a lower overall taxonomic richness and abundance. Nevertheless, the number of taxa and the diversity of macroinvertebrates remained relatively high across the basin, especially in the delta assemblage. Therefore, the LMB deserves a particular attention for conservation.

Taxonomic groups

PART II: PUBLICATIONS

Introduction

Tropical regions are remarkably rich in biodiversity (Sodhi et al., 2004), with 16 out of the 25 terrestrial biodiversity hotspots of the world being located in the tropical zone [START_REF] Mittermeier | Hotspots: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions. Agrupación Sierra Madre[END_REF]. Southeast Asia has a unique geological history (Sodhi et al., 2004), and can be separated into four biodiversity hotspot regions: The Philippines, Sundaland, Wallacea and Indo-Burma [START_REF] Mittermeier | Hotspots: Earth's Biologically Richest and Most Endangered Terrestrial Ecoregions. Agrupación Sierra Madre[END_REF]. Through Indo-Burma, the Mekong River, which is the longest river in Southeast Asia, flows from its source in the Tibetan plateau to the South China Sea (Zalinge et al., 2003). This river harbours diverse biotic communities and is the breeding area of numerous endemic, threatened and endangered species. The aquatic and terrestrial communities of the Mekong form a biologically important food web that supports a high biodiversity (Sodhi et al., 2004). Moreover, it is an economically important region, since aquatic fauna such as fish (∼1200 species, [START_REF] Rainboth | Species identification field guide for fishery purposes[END_REF], molluscs, crustaceans and insects, are all highly dependent on this basin as a breeding ground (Davidson et al., 2006;Zalinge and Thuok, 1998).

Macroinvertebrates are a key component of freshwater ecosystems (Bogan, 2008;Palmer et al., 1997). In river systems, macroinvertebrate communities are differently organized and structured along environmental gradients. Many studies have shown that in the upstream regions, the invertebrate communities are mainly characterized by a high abundance of insects (e.g. Ephemeroptera, Odonata and Trichop-tera) (Arab et al., 2009;Collier and Lill, 2008;[START_REF] Jiang | Longitudinal patterns of macroinvertebrate communities in relation to environmental factors in a Tibetan-Plateau river system[END_REF], whereas the downstream communities are characterized by a high abundance of molluscs, crustaceans and worms (Arscott et al., 2005; Collier and Lill, 2008;Pérez-Quintero, 2011). This spatial variation in macroinvertebrate community composition is influenced by trophic variables, e.g. trophic level or food sources (Cai et al., 2012;Nicola et al., 2010); water chemistry, e.g. dissolved oxygen, water conductivity, pH (Al-Shami et al., 2013;[START_REF] Heino | Biodiversity of aquatic insects: spatial gradients and environmental correlates of assemblage-level measures at large scales[END_REF][START_REF] Kumar | The distribution and diversity of benthic macroinvertebrate fauna in Pondicherry mangroves, India[END_REF] and physical river conditions, e.g. river depth, river width, size of catchment area and characteristics of land cover (Allan, 2004;[START_REF] Beisel | Stream community structure in relation to spatial variation: the influence of mesohabitat characteristics[END_REF]Chadwick et al., 2006;Cortes et al., 2013).

The Lower Mekong Basin (LMB) is characterized by a long and large floodplain (Eastham et al., 2008) and is known for its high biodiversity (Sodhi et al., 2004). However, the knowledge of macroinvertebrates in the LMB is limited. Given that this river basin is being impacted by various anthropogenic disturbances such as agricultural activities, aquaculture, urbanization and mining (Köhler et al., 2012;Nhan et al., 2007;Sodhi et al., 2004), there is an urgent need to understand the patterns of spatial organization and community structure of macroinvertebrates in this basin and their relation to environmental factors. Up to date, only a few studies have been published on the basin, e.g. community structure and composition of littoral invertebrates in the Mekong delta (Wilby et al., 2006), the diversity and distribution of crustaceans and molluscs in the Indo-Burma region (Cumberlidge et al., 2011;Köhler et al., 2012), and the assessment of water quality using physicochemical variables and benthic macroinvertebrates in the Chi river in Thailand (Kudthalang and Thanee, 2010). Yet, no attempt has been made to examine the spatial patterns of macroinvertebrate assemblages and their relation to key environmental variables at a large scale.

The objectives of the present study were i) to investigate patterns of spatial variation in macroinvertebrate assemblages in the Lower Mekong Basin, ii) to analyse the variability of macroinvertebrate composition among the assemblages, and to determine their key indicator taxa (the most representative taxa) for particular assemblages, and iii) to identify the important environmental variables that are associated with the macroinvertebrate assemblages in the basin. We expected that physical conditions of habitats, compared to other measured variables, have a strong correlation to the composition and diversity of macroinvertebrate assemblages.

Materials and methods

Study area

The Mekong River Basin (MRB) is divided into the Upper Mekong Basin (UMB) and the Lower Mekong Basin (LMB). The UMB on the Tibetan plateau in China is composed of narrow, deep gorges and small, short tributaries, whereas the LMB stretches from Yunnan province in South China to the delta in Vietnam and covers approximately 70% of the total length of the MRB (Eastham et al., 2008). The LMB consists of a large floodplain and long, broad tributaries and it drains more than 76% of the Mekong basin. The climate of the LMB is dominated by a tropical monsoon rainfall system, which is characterized by a dry (November -April) and a wet (May -October) season generated by the northeast monsoon and the southwest monsoon, respectively. The most intensive rainfall falls from July to September, while the lowest precipitation is observed between January and April (Adamson et al., 2009). The annual rainfall of the LMB varies from 1000-1600 mm in the driest regions to 2000-3000 mm in the wettest regions (Hoanh et al., 2003). A higher precipitation is found in the eastern mountainous regions of Laos and in northeast Thailand (Eastham et al., 2008).

The largest floodplain water body of the LMB is the Tonle Sap Lake (TSL) in Cambodia (Adamson et al., 2009), which is the largest freshwater lake in Southeast Asia (Sarkkula et al., 2003). The TSL is connected to the Mekong through the Tonle Sap River, and thus creating an exceptional hydrological cycle. In the wet season, the TSL receives an excess water from the Mekong River and expands its surface area from 2500 km 2 to 15,000 km 2 . In the dry season when the rain ceases and water levels drop in the Mekong, a reverse flow occurs; the drained water from the TSL flows to the Mekong delta (Arias et al., 2011). The Mekong delta is characterized by a number of man-made canals, which are mostly used for domestic and agricultural activities (Kummu et al., 2008).

Data collection and processing

Benthic macroinvertebrates were sampled at 60 sampling sites along the main channel of the LMB and its tributaries by the Mekong River Commission (MRC) (Fig. 1). This sampling was carried out once a year in March during the dry season from 2004 to 2008. At each sampling site, macroinvertebrates were sampled from three locations in the benthic zone: near the left and right banks, and in the middle of the rivers. At each location, a minimum of three samples (where intersample variability is low, e.g. tributaries) and a maximum of five samples (where inter-sample variability is higher, e.g. the main channel and the delta) were collected using a Petersen grab sampler which has a sampling area of 0.025 m 2 . With the grab sampler, four sub-samples were taken and pooled to give a single sampling unit covering a total area of 0.1 m 2 . In total, between nine (3 samples × 3 locations) and fifteen (5 samples × 3 locations) pooled samples were collected at each sampling site. Each pooled sample was rinsed using a sieve (300 μm mesh size). In the field, the samples were sorted and then preserved by adding 10% formaldehyde to obtain a final concentration of about 5%. In the laboratory, they were identified to the lowest taxonomic level possible and counted using a compound microscope (40-1200 magnification) or a dissecting microscope (16-56 magnification). Macroinvertebrate abundance data per sampling unit was averaged across all samples (between 9 and 15 samples) collected from each sampling site.

At the sampling sites, geographic coordinates and altitude were determined with a GPS (Garmin GPS 12XL). River width was measured in the field using a Newcon Optik LRB 7 × 50 laser rangefinder. Other physical-chemical variables were measured at the three locations where macroinvertebrates were sampled. River depth was measured using a line metre. With a handheld water quality probe (YSI 556MP5), water temperature, dissolved oxygen, pH and water conductivity were measured at the surface (0.1-0.5 m) and at a depth of 3.5 m or at a maximum depth of the river (wherever less than 3.5 m) and then the average value was recorded for each location. Water transparency was measured with a Secchi disc by lowering it into the water and recording the depth at which it was no longer visible. The physical-chemical data of each sampling site was represented by the average value across the three sampling locations. The surface area and land cover data of watersheds drained at each sampling site were determined using a Geographic Information System (ArcGIS 10.4, ESRI). Geographic data (ArcGIS shapefiles) about the LMB (land cover types, river networks, basin boundaries and subcatchments derived from topographical maps) was provided by the MRC.

In total, 108 samples were collected from the 60 sampling sites. In 2008, 3 sampling sites were sampled further away from their original sampling coordinates, and thus we considered them as different sampling sites (see Supporting Information Appendix S1). Therefore, a total of 63 sampling sites were taken into account in the present study. Because of unequal sampling efforts (i.e. unequal and different number of samples at each site during the 5-year sampling period) and missing values of environmental variables, we used median values from the collected data to represent each site in our analysis, as suggested by McCluskey and Lalkhen (2007). More precisely, for the sites (49 sites) which were sampled only one time during the 5-year sampling period (see Appendix S1), the one sample collected from each site was used as the representative sample. The remaining 14 sites contained two samples (11 sites), three samples (2 sites) and four samples (1 site). From these 14 sites, the median values were used to represent the corresponding sites.

Data analysis

We used Ward's hierarchical clustering method to group the 63 sampling sites into different clusters. The Bray-Curtis dissimilarity distance (Legendre and Legendre, 2012) of macroinvertebrate samples between sites was calculated using the Hellinger transformation in the package "vegan" of R (Rao, 1995). Macroinvertebrate richness, abundance, Shannon-Wiener diversity were used to analyse the variability in macroinvertebrate assemblage composition among the different clusters. An "assemblage" of macroinvertebrates refers to different taxonomic groups (i.e. annelids, crustaceans, insects and molluscs) living in a given geographic area (e.g. a defined cluster) [START_REF] Fauth | Simplifying the jargon of community ecology: a conceptual approach[END_REF]. The indicator taxa of macroinvertebrate assemblages in the defined clusters were determined using the Indicator Value (IndVal, Dufrene and Legendre, 1997) with the package "labdsv" of R (Roberts, 2013). The Indicator Value of a taxon is an index ranging from 0 to 1, indicating the least to the most important taxa occurring in a group of sites. A value of 1 is obtained when every individual of the taxon is found only in the group and when it occurs in all sites of that group. A high number of taxa with significant Indicator Values may provide information on the habitat they prefer to share. Taxa having Indicator Values with a pvalue ≤ 0.01 were retained as the most important taxa representing the assemblage in each cluster (consisting of a group of sites).

Significant differences in macroinvertebrate assemblage composition (i.e. richness, abundance, diversity) and in environmental conditions among different clusters were tested using a one-way ANOVA or a Kruskal-Wallis test. One-way ANOVA was applied when residuals of the models were normal (Shapiro-Wilk test, p > 0.05, and homoscedastic (Bartlett's test, p > 0.05)); otherwise, the non-parametric test (Kruskal-Wallis) was used. To evaluate the differences between clusters, multiple pair-wise comparisons were conducted using a Kruskal-Wallis test.

Linear Discriminant Analysis (LDA) was performed, using the package "ade4" of R (Chessel, 2006), to assess which measured environmental variables best accounted for the differences in macroinvertebrate assemblage composition among clusters grouped by the hierarchical clustering. Before performing the LDA, environmental variables were tested for multivariate homogeneity of within-group covariance (Borcard et al., 2011). The contribution of each variable to the discrimination between clusters was represented by the standardized factorial coefficient, projected as an arrow on the LDA plot. Cross-validation was used to evaluate the predictive performance of the LDA model. The overall quality of the model and the reliability of the prediction for the member of each cluster were evaluated by Cohen's Kappa Statistic (Kappa). All statistical analyses were performed in R (R Core Team, 2013).

Results

Macroinvertebrate diversity

Two hundred and ninety-nine macroinvertebrate taxa accounting for a total of 21,810 individuals were identified; of which, 131 taxa (44%) were insects, 98 (33%) were molluscs, 38 (13%) were crustaceans and 32 (10%) were annelids. The most commonly found insect orders were Diptera (37 taxa), Ephemeroptera (32), Odonata ( 22) and Trichoptera (20). For molluscs, most taxa belonged to the order Unionida (18), Veneroida (15) and Caenogastropoda (50); whereas for crustaceans, most taxa belonged to the order Decapoda (18) and Amphipoda (9). Annelids were mainly represented by the order Haplotaxida (15). Shannon diversity of macroinvertebrates at each site ranged from 0.8 to 3.3.

Macroinvertebrate assemblages and indicator taxa

Based on the Bray-Curtis dissimilarity and the hierarchical analysis, the 63 sites were grouped into four clusters (Fig. 2). Cluster I consisted of 11 sites located in the Mekong delta in Vietnam; cluster IIa, of 11 upstream sites in Laos and Thailand; cluster IIb1, of 15 sites mostly in middle part in Cambodia, a few sites in Laos and Thailand, and 1 site in Vietnam; and cluster IIb2, of 26 sites mostly located in the tributaries of the LMB (Fig. 1). The number of sampling sites and environmental variables characterizing each cluster are shown in Table 1.

The four clusters, which were later in the text considered as four different assemblages, had a significant difference in macroinvertebrate richness (Kruskal-Wallis H = 27.7, p < 0.001), abundance (Kruskal-Wallis H = 23.7, p < 0.001) and Shannon diversity (Kruskal-Wallis H = 11.7, p < 0.01) (Fig. 3, Table 1). The highest macroinvertebrate richness, abundance and diversity were found for assemblage I, followed by assemblage IIa and IIb1. The lowest macroinvertebrate richness, abundance and diversity were found for assemblage IIb2 (Fig. 3A andB, Table 1). Molluscs represented the highest proportion of richness and abundance for assemblage I and IIb1. For assemblage IIa and IIb2, insects represented the highest proportion, followed by molluscs. Crustaceans and annelids made up the lowest proportion of richness and abundance for assemblage IIa, IIb1 and IIb2, but they represented a higher proportion compared to insects for assemblage I (Fig. 3C andD).

The four assemblages were represented by a different number of indicator taxa: 53, 14, 2 and 2 taxa, for assemblage I, IIa, IIb1 and IIb2, respectively. Molluscs, crustaceans and annelids represented a relatively high number of indicator taxa compared to insects for assemblage I. The assemblage IIa was mostly represented by insects, while assemblages IIb1 and IIb2 were represented by two different indicator taxa. The detailed information on indicator taxa for each assemblage is provided in Table 2.

Relationship between macroinvertebrate assemblages and environmental variables

Most of the environmental variables showed significant differences between the four macroinvertebrate assemblages. The detailed information on the environmental differences between the assemblages is provided in Table 1. The results of the LDA used to discriminate the macroinvertebrate assemblages based on the physical-chemical and land cover types are shown in Fig. 4. The factorial axes 1, 2 and 3 respectively explained 65.7%, 23.1% and 11.2% of the total variance of the data. The global performance of the LDA was high (Kappa = 0.86). The predictive reliability (Kappa) for assemblage I, IIa, IIb1 and IIb2 was 0.90, 0.84, 0.80 and 0.90, respectively. Along axis 1, assemblage I was situated opposite to assemblage IIa and IIb2. Assemblage I was positively correlated with the surface area of watershed, river depth, river width and water temperature, but negatively associated with altitude and dissolved oxygen. Whereas assemblage IIa was positively correlated with water conductivity and urban areas, and assemblage IIb2 was positively linked to high altitude, dissolved oxygen, Secchi depth, wood-/shrub-land and evergreen forests. Based on axis 1 and 3, assemblage IIb1 was positively linked to inundated, wetland and agricultural areas (Fig. 4A-D).

Discussion

Spatial variation in macroinvertebrate assemblages and habitat characteristics

Macroinvertebrate assemblages can vary substantially according to habitat characteristics, particularly with large variations in habitat types and resources (Costa and Melo, 2008;[START_REF] Ilg | Extreme flood events favour floodplain mollusc diversity[END_REF][START_REF] Vaughn | The functional role of burrowing bivalves in R. Sor et al. Limnologica[END_REF]. The combination of environmental factors such as water velocity, river depth, organic material load, watershed and substratum structure, together with the high amount of food resources and different water quality conditions found in the LMB (Chea et al., 2016), are responsible for determining macroinvertebrate assemblages (Al-Shami et al., 2013;[START_REF] Beisel | Stream community structure in relation to spatial variation: the influence of mesohabitat characteristics[END_REF][START_REF] Lamouroux | Biological traits of stream macroinvertebrate communities: effects of microhabitat, reach, and basin filters[END_REF][START_REF] Lorenz | Upstream river morphology and riparian land use overrule local restoration effects on ecological status assessment[END_REF][START_REF] Pan | A floodplain-scale lake classification based on characteristics of macroinvertebrate assemblages and corresponding environmental properties[END_REF]. Moreover, macroinvertebrate assemblages are also structured by different land cover characteristics, e.g. agricultural land, wood-/shrub-land and urban areas (Allan, 2004;Cortes et al., 2011Cortes et al., , 2013)). Therefore, these physical-chemical variables and land cover types are likely the main factors determining the spatial organization and composition of macroinvertebrate assemblages in this large floodplain river.

The spatial organization of macroinvertebrate assemblages is mainly related to the relationship between habitat conditions and life-history traits at different spatial scales or along the longitudinal gradients of rivers (Collier and Lill, 2008;[START_REF] Tonkin | Do productivity and disturbance interact to modulate macroinvertebrate diversity in streams?[END_REF]. In this regard, key physical-chemical variables and land cover types can act as environmental filters, which allow the most suited taxonomic groups of macroinvertebrates to be present in each assemblage [START_REF] Menezes | Beyond taxonomy: a review of macroinvertebrate trait-based community descriptors as tools for freshwater biomonitoring[END_REF]. We found a higher insect richness and abundance for assemblage IIb2, which is mostly related to the tributaries and some sites of the main channel at the upstream parts. This assemblage is represented only by one insect and one annelid indicator taxon, which were widely distributed in the assemblage. This result suggests that natural stress in the mountainous and tributary system is probably too high for some taxa to occur (Feld and Hering, 2007). The fact that this assemblage was dominated by insects agrees with the expectation that insect orders mainly inhabit mountainous and shaded habitat regions (Ferrington, 2008;[START_REF] Heino | Biodiversity of aquatic insects: spatial gradients and environmental correlates of assemblage-level measures at large scales[END_REF]Suhaila and Che Salmah, 2014). The presence of evergreen forests, bamboo, crops and woody debris may provide food resources and egg deposition sites for many insect taxa [START_REF] Medhurst | Headwater streams and forest management: does ecoregional context influence logging effects on benthic communities?[END_REF][START_REF] Sweeney | Effects of streamside vegetation on macroinvertebrate communities of White Clay Creek in eastern North-America[END_REF]. Clear water and a high level of dissolved oxygen found in high altitude streams characterize the tributaries, which are preferred by insect taxa (Collier and Lill, 2008;Dobson et al., 2002;Królak and Korycińska, 2008). Moreover, sampling sites in this assemblage are comparable to the Holarctic region where Diptera are more abundant (Ferrington, 2008), and are characterized by the increasing latitude where most Ephemeroptera occur [START_REF] Pearson | Gradients in regional diversity of freshwater taxa[END_REF].

Macroinvertebrate assemblage IIa was mainly related to the upstream sites situated along the main channel and some tributaries. In this assemblage, insects were also the dominant group in terms of richness, abundance and indicator taxa (7 out of 14 taxa, Table 2). However, annelid and mollusc indicator taxa were also found in this assemblage. A high water conductivity (22.8 ± 6 mS/m), resulting from anthropogenic disturbance (e.g. runoff from the surrounding agricultural areas and the discharge of urban wastewaters), could be one of the factors responsible for the community composition found, as it has been demonstrated that conductivity has a strong impact on macroinvertebrate diversity [START_REF] Kumar | The distribution and diversity of benthic macroinvertebrate fauna in Pondicherry mangroves, India[END_REF]Lods-Crozet et al., 2001). The high contribution of insects to this assemblage could be explained by a high level of dissolved oxygen present at some sites (Table 1), as indicated in previous studies (Dobson et al., 2002;Królak and Korycińska, 2008), while disturbance-tolerant taxa of Oligochaeta and Gastropoda may be well adapted at sites having a high water conductivity (Feld and Hering, 2007;Wang et al., 2012). Moreover, the main channel generally provides a higher variation in habitats and a higher nutrient and sediment load, which are also favourable for annelids and molluscs [START_REF] Haag | Role of ecological factors and reproductive strategies in structuring freshwater mussel communities[END_REF].

In-between the upstream and downstream parts of the river basin, the macroinvertebrate assemblage was dominated by molluscs, followed by insects. Characteristics of sites in the mid-reaches such as rock and deciduous forests play an important role in supplying organic matter, which is preferred by snails (grazers) and insect collectors and shredders (Thorp and Delong, 1994). Many sites belonging to this assemblage that are located around the Tonle Sap Lake, Cambodia, are characterized by inundated, agricultural land and wetland (Fig. 4C and D, Table 1). These site specific features could promote the richness and abundance of molluscs found at these sites. Moreover, it is noteworthy that the areas around the Tonle Sap Lake support a high mollusc production (Ngor et al., 2016). Therefore, it is not surprising that molluscs made up the largest proportion of the macroinvertebrate assemblage and that snail species were put forward as indicator species for this assemblage. The tributary sites that are located upstream, characterized by shaded areas and surrounded by large aquatic plants, are likely responsible for the second highest proportion of insect richness and abundance and the presence of an insect indicator taxon.

In the LMB, downstream sites were mainly associated with a large surface area of watersheds, large floodplains, wide and deep rivers. These conditions promote a high nutrient and sediment load from the upstream parts [START_REF] Blair | From bedrock to burial: the evolution of particulate organic carbon across coupled watershed-continental margin systems[END_REF], which can support many species and a high abundance of macroinvertebrates (Cai et al., 2012;Nicola et al., 2010). Moreover, the high temperature related to the downstream sites may also enhance the richness and abundance of molluscs [START_REF] Vaughn | The functional role of burrowing bivalves in R. Sor et al. Limnologica[END_REF]. High nutrient input and high temperature due to sunlight observed at the downstream delta provide optimal conditions for the phytoplankton community to reach a high abundance [START_REF] Hecky | Nutrient limitation of phytoplankton in freshwater and marine environments: a review of recent evidence on the effects of enrichment[END_REF][START_REF] Statzner | Questions and comments on the river continuum concept[END_REF]Vannote et al., 1980), which in turn enhance the richness and abundance of phytoplankton feeders such as molluscs [START_REF] Vaughn | The functional role of burrowing bivalves in R. Sor et al. Limnologica[END_REF]. However, although the downstream assemblage was dominated by molluscs, annelids and crustaceans, dipteran insects (14 out of 33 occurring taxa) also made up a relatively high proportion. As a result, each taxonomic group in this assemblage is represented by many indicator taxa. This suggests that the downstream sites, the delta, provide good habitat conditions to support a diverse fauna. The high richness and abundance of all taxonomic groups and the high diversity in this assemblage might be explained 1) by the large surface area of the watershed and the large width and depth of the river, which may provide an optimal nutrient load and different mesohabitats and microhabitats, respectively (Al-Shami et al., 2013;[START_REF] Haag | Role of ecological factors and reproductive strategies in structuring freshwater mussel communities[END_REF][START_REF] Jacobsen | Structure and diversity of stream invertebrate assemblages: the influence of temperature with altitude and latitude[END_REF]Mereta et al., 2012;Sedell et al., 1989), and 2) by a slow water flow, as reflected by the low change in altitude (6.6 ± 1.8 m) of the sampling sites, which may provide good habitat conditions for less-mobile taxa like molluscs, annelids and crustaceans [START_REF] Castella | A methodological approach to the use of terrestrial invertebrates for the assessment of alluvial wetlands[END_REF][START_REF] Haag | Role of ecological factors and reproductive strategies in structuring freshwater mussel communities[END_REF].

The compositional and diversity differences found for the four macroinvertebrate assemblages were related to spatial environmental heterogeneity observed from the upstream to the downstream parts of the LMB (Dobson et al., 2002;Heino et al., 2005;Salmah et al., 2014;Wang et al., 2012). This can be explained by the proportion of the richness, abundance, diversity and the represented indicator taxa of the different taxonomic groups for each assemblage. Although assemblage IIa, IIb1 and IIb2 were characterized by a low and similar macroinvertebrate richness, abundance and diversity, each taxonomic group (i.e. annelids, crustaceans, molluscs and insects) contributed differently to the overall macroinvertebrate community composition for the particular assemblages. On the other hand, assemblage I had a high richness, abundance, diversity and a high number of indicator taxa. Therefore, these findings clearly suggested that the community composition of the four assemblages, particularly assemblage I, is distinctly organized and structured along the environmental gradients of the LMB.

Relationship between assemblages and key environmental variables

Site-specific characteristics, as mentioned above, are known to influence local community and composition of macroinvertebrates. At an assemblage level, connection in habitats and mesohabitats and variation in river morphology strongly affect community structure and composition from the upstream to the downstream parts of rivers (Mazão and Bispo, 2016;Thorp and Delong, 1994;Zilli and Marchese, 2011). Our study demonstrated that among the measured environmental variables, physical conditions and land cover types are the key factors determining macroinvertebrate assemblages along the LMB. These results support the previous finding that macroinvertebrate community structures in the Chishui river basin and in the riverconnected lakes of the Yangtze river are strongly influenced and well predicted by physical habitat variables, e.g. water depth and altitude [START_REF] Jiang | Structure of macroinvertebrate communities in relation to environmental variables in a subtropical Asian river system[END_REF][START_REF] Pan | A floodplain-scale lake classification based on characteristics of macroinvertebrate assemblages and corresponding environmental properties[END_REF]. In other tropical/sub-tropical (Pearl and Qiangtang rivers, southern China) and temperate river basins (e.g. the Olo, Corgo, Pinhao and Tua rivers in northern Portugal), characteristics of land cover have been reported to better explain and predict benthic macroinvertebrate communities (Allan, 2004;Cortes et al., 2011Cortes et al., , 2013)).

In the LMB, the surface area of watershed, river depth and altitude are the most important variables discriminating the upstream assemblages (IIa and IIb2) from the downstream assemblage (I). This discrimination is clearly explained by the axis 1 of the LDA accounting for 65.7% of the total variance of the data and by a high predictive reliability for each particular assemblage. Because of the observed physical conditions, together with other related variables (e.g. water temperature, dissolved oxygen and land cover types), the assemblages were characterized by a very different community structure, composition and diversity.

The intermediate assemblage (IIb1) in-between was positively linked to three types of land cover: agricultural land, inundated and wetland. Along axis 3 of the LDA model, these land cover types are the key factors best explaining the macroinvertebrate composition of assemblage IIb1 (Fig. 4C andD). However, the proportion of inundated and wetland constituted a small percentage of the land cover, and together with agricultural land, they only explained a limited amount of the variance of the data for axis 3 of the LDA model (11.2%). As a result, the predictive performance for this assemblage (Kappa = 0.80) was not as high as the performance for the other assemblages. A broad range of values of physical-chemical variables (e.g. dissolved oxygen, river width and depth, the surface area of watershed and altitude) found at sites belonging to this assemblage could also explain the performance yielded. On the other hand, assemblage IIa was associated with two important factors (water conductivity and urban areas) according to the LDA axis 2, which explained a higher variance of the data (23.1%). Therefore, these two factors, together with other measured environmental variables, could correctly predict this assemblage with a higher reliability (Kappa = 0.85). The possible explanation for the association found is that many sites that belong to this assemblage are located along the tributaries (e.g. the Mun and Chi river basins in Thailand) and the main channel where cities were built. The tributaries are surrounded by intensified agriculture and the cities are exposed to a high level of anthropogenic disturbance (Dao et al., 2010;Kudthalang and Thanee, 2010). Therefore, the runoff from the surrounding agricultural areas and the discharge of urban wastewaters may cause the increase in water conductivity (Wetzel, 2001).

Recommendations for management of the LMB

Macroinvertebrate diversity is of great importance for ecosystem structure and functioning as it interacts with both biotic and abiotic factors. Crustaceans and insects are generally responsible for regulating decomposition, shredding detritus and bioturbation. Molluscs also contribute to bioturbation, sediment formation and filtering of water, while most annelids regulate decomposition and autotrophs although some of them also promote bioturbation and sediment formation (Palmer et al., 1997). The different processes enhanced by these taxa promote ecosystem functioning in the LMB, particularly in its lower reaches. Assemblage I, occupying the downstream part of the LMB, was characterized by a high macroinvertebrate diversity and was represented by a number of indicator taxa from each group of macroinvertebrates, which could be expected to increase ecosystem stability. These taxa are key components of the ecosystem and of the food web supporting a rich biodiversity of fish, which are a main source of proteins for local people. Thus, a change in the composition of one of these taxa may result in disproportionate or unexpected responses of other taxa in the LMB (Naeem, 1998), and consequently alter ecological processes (Covich et al., 1999) such as organic matter processing (Palmer et al., 1997).

For management purposes, indicator taxa identified in this study representing the ecological importance of each site and each community are of great significance. By using the results of this study such as the indicator taxa and the correlated environmental factors for each assemblage, habitat quality and suitability for the most vulnerable native taxa could be more conveniently monitored, which could provide an advantage for decision-making concerning conservation and management of these keystone taxa. Moreover, this study can serve as a baseline for future research in a biodiversity hotspot region that has not been investigated intensively.

Conclusion

The macroinvertebrate assemblages in the Lower Mekong Basin (LMB) were characterized by an increasing richness, abundance and diversity from the upstream to the downstream sites. Sites located in the upper tributaries (cluster IIb2) and upper main channel (cluster IIa) were characterized by a high altitude, high levels of dissolved oxygen, a high water conductivity, large fractions of wood-/shrub-land and evergreen forests. Macroinvertebrate assemblages in these two clusters were dominated by insects. In the Mekong delta (cluster I), the macroinvertebrate assemblage largely consisted of a great number of molluscs and an average number of annelids, crustaceans and dipteran insects. This assemblage could be discriminated based on a large surface area of watershed, the river depth and width and a high water temperature. The assemblage (cluster IIb1) found in-between the tributaries and the upper main channel and the delta region was associated with sites characterized by inundated, wetland and agricultural land and was represented most by molluscs. Overall, our study found that the number of macroinvertebrate taxa and the diversity remain relatively high across the basin, especially in the delta region which can be considered a hotspot zone for biodiversity. Therefore, the LMB deserves a particular attention for conservation.

Introduction

The variation in community composition among sites, or beta (β) diversity (Legendre and De Cáceres, 2013;[START_REF] Whittaker | Vegetation of the Siskiyou Mountains, Oregon and California[END_REF], is of primary interest to community ecology. Beta diversity is an important component of biodiversity as it links local (α) to regional (γ) diversity, and it varies as a function of the spatial scales and gradients of the study areas [START_REF] Anderson | Navigating the multiple meanings of β diversity: a roadmap for the practicing ecologist[END_REF]Legendre and Legendre, 2012;[START_REF] Whittaker | Evolution and measurement of species diversity[END_REF][START_REF] Whittaker | Vegetation of the Siskiyou Mountains, Oregon and California[END_REF]. Therefore, understanding the variation in species composition among sites, i.e. β diversity, enables community ecologists to disclose evolutionary and ecological processes at work in a community of interest [START_REF] Valdujo | Environmental correlates of anuran beta diversity in the Brazilian Cerrado[END_REF], by analyzing and testing such processes in a way that indicates how they affect and maintain biodiversity in the ecosystem (Legendre and De Cáceres, 2013).

Biotic interactions in communities, e.g. intra-specific competition, have been reported to influence β diversity [START_REF] Matthiessen | Dispersal frequency affects local biomass production by controlling local diversity[END_REF][START_REF] Valdujo | Environmental correlates of anuran beta diversity in the Brazilian Cerrado[END_REF]. In identical ecological patches, interactions among species within their own taxonomic group or between different taxonomic groups can lead to different patterns of β diversity [START_REF] Hillebrand | Regional and local impact on species diversity -from pattern to processes[END_REF][START_REF] Tonkin | Variable elements of metacommunity structure across an aquatic-terrestrial ecotone[END_REF], and thus affect ecosystem functioning. Environmental gradients, habitat heterogeneity [START_REF] López-González | Ecological biogeography of Mexican bats: the relative contributions of habitat heterogeneity, beta diversity, and environmental gradients to species richness and composition patterns[END_REF], and natural and humanderived disturbances [START_REF] Lamy | Understanding the spatiotemporal response of coral reef fish communities to natural disturbances: insights from beta-diversity decomposition[END_REF]Legendre and Salvat, 2015) have been shown to also influence β diversity. For aquatic macroinvertebrates, β diversity is mainly related to drainage basins and within-stream environmental factors, while it has been reported not to be significantly related to habitat degradation, eutrophication, longitude and altitude (Friberg et al., 2010;Md Rawi et al., 2013). However, longitude and altitude have been found to be substitute variables (proxies) for major drivers patterning β diversity of macroinvertebrates at broad geographical scales (J. Wang et al., 2012). The environmental variables related to the geographical proxies may play important roles in structuring the broad-scale pattern of β diversity in a given region.

Several papers have reported patterns of β diversity in tropical ecosystems for plants and vertebrate animals (e.g. [START_REF] Legendre | Partitioning beta diversity in a subtropical broad-leaved forest of China[END_REF][START_REF] López-González | Ecological biogeography of Mexican bats: the relative contributions of habitat heterogeneity, beta diversity, and environmental gradients to species richness and composition patterns[END_REF][START_REF] Mena | Species turnover on elevational gradients in small rodents[END_REF][START_REF] Wearn | Grain-dependent responses of mammalian diversity to land use and the implications for conservation setaside[END_REF]. β diversity of macroinvertebrates has also recently been analyzed by several authors (e.g. Costa and Melo, 2008;[START_REF] Leigh | Hydrological connectivity drives patterns of macroinvertebrate biodiversity in floodplain rivers of the Australian wet/dry tropics[END_REF][START_REF] Ligeiro | Spatial scale and the diversity of macroinvertebrates in a Neotropical catchment[END_REF], but only a few studies have taken place in South-East Asia (e.g. Al-Shami et al., 2013, Salmah et al., 2014). As the ecosystems in that region are highly endangered and heavily impacted by human disturbances (Salmah et al., 2014;Sodhi et al., 2004;[START_REF] Strayer | Freshwater biodiversity conservation : recent progress and future challenges[END_REF], assessing the patterns of macroinvertebrate β diversity and their relationships to geo-environmental factors and to related biotic communities is urgently needed.

The Mekong River Basin is divided into Upper and Lower Mekong Basins (LMB). The LMB, covering an area of about 609,000 km 2 (77% of the whole basin) (Zalinge et al., 2003), includes portions of four densely populated countries: Thailand, Laos, Cambodia and Vietnam. This basin harbors diversified communities of fish and invertebrates, forming biologically important food webs that support high biodiversity (Sodhi et al., 2004). Many aquatic taxonomic groups such as fishes, mollusks, crustaceans and insects are highly dependent on this basin as a breeding ground (Davidson et al., 2006;Zalinge and Thuok, 1998). In spite of high suspected biodiversity in the LMB, the β diversity and community patterns of its aquatic taxonomic groups, particularly the macroinvertebrates, have seldom been studied. The biomonitoring surveys conducted by the Mekong River Commission (MRC) represent the only major work conducted on aquatic macroinvertebrates in the LMB. In this study, we used this bio-monitoring data to explore the β diversity pattern of aquatic macroinvertebrates. Analyzes of this dataset, collected from sites sampled over 5 successive years (2004)(2005)(2006)(2007)(2008), should significantly contribute to increase our scientific knowledge of the LMB.

Beta diversity can be computed in different ways [START_REF] Koleff | Measuring beta diversity for presence-absence data[END_REF][START_REF] Whittaker | Vegetation of the Siskiyou Mountains, Oregon and California[END_REF]. A classical approach is to compute β diversity as β = γ/ᾱ, where γ is the total number of species in a given region and ᾱ is the average number of taxa for a sample set within the region [START_REF] Whittaker | Vegetation of the Siskiyou Mountains, Oregon and California[END_REF]. This classical measurement is still preferred by many authors [START_REF] Higgins | Patterns of functional and taxonomic organization of stream fishes: inferences based on α, β, and γ diversities[END_REF][START_REF] Jost | Partitioning diversity into independent alpha and beta components[END_REF][START_REF] Sor | Diversity of Monogononta rotifer species among standing waterbodies in northern Cambodia[END_REF] although new approaches have been developed [START_REF] Anderson | Navigating the multiple meanings of β diversity: a roadmap for the practicing ecologist[END_REF][START_REF] Legendre | Analyzing beta diversity: partitioning the spatial variation of community composition data[END_REF]Legendre and De Cáceres, 2013).

In this study, we used the total variance of the macroinvertebrate communities among the study sites of the LMB as a measure of beta diversity (BDTotal) and partitioned it into "Local Contributions to Beta Diversity" (LCBD) and "Species Contributions to Beta Diversity" (SCBD) (Legendre and De Cáceres, 2013). We identified the important taxa contributing most to total β diversity, i.e. those with high among-site variance, as well as the geo-environmental factors that were associated with the macroinvertebrate communities throughout the sites. In addition, we investigated the influence of the LCBD indices of the component communities (i.e. annelids, crustaceans, mollusks and insects) on the LCBD indices of the global macroinvertebrate community composition (including all component communities). Our questions of interest are the following: 1) Is there a moderate or a large amount of variation in macroinvertebrate community composition among the sites in the LMB? 2) What are the taxa that contribute most to the total β diversity? We expect the important taxa of annelids, crustaceans and mollusks, measured as richness and abundance, to be associated with sites located downstream, whereas the important taxa of insects should be associated with sites located farther upstream, as has been shown by Arscott et al. (2005) and Królak and Korycińska (2008). 3) Are there sites that have exceptionally unique taxonomic compositions? We hypothesize that some sampling locations exhibit significant uniqueness in taxonomic composition. 4) What are the geo-environmental conditions that characterize the sites with significant LCBD indices? We expect the LCBD indices to increase with river width and pH, following the β diversity patterns found in tropical streams in Malaysia (Al-Shami et al., 2013), and decrease with latitude and altitude, following the β diversity patterns observed in major geographical diversity gradients (J. Wang et al., 2012). 5) What are the component communities that mainly influence the LCBD indices of the global macroinvertebrate communities?

Materials and methods

Macroinvertebrate and geo-environmental variables collection

From 2004 to 2008, the Mekong River Commission (MRC) conducted biomonitoring surveys and sampled macroinvertebrates at 60 sites along the LMB once a year in March during the dry season (Fig. 1). To harmonize the data being collected, the sampling locations were selected from different habitats such as those in or close to villages or towns, at rivers with substantial shipping, next to crop fields and meadows with livestock, upstream or downstream of dams or weirs, and at more pristine areas surrounded by forest with only few houses. At each sampling site, benthic macroinvertebrates and geo-environmental variables were collected at the same time. For the detailed information on the collection process, we refer to Sor et al. (2017).

In 2008, 3 sampling sites were sampled farther away from their original sampling coordinates, and thus they were regarded as new sampling sites (see Appendix T1 in Part I: Synthesis). Therefore, we considered a total of 63 sampling sites in the present study.

Data processing and statistical methods

For the 63 sampling sites, 108 samples of biological and geo-environmental variables were available. Due to unequal sampling efforts, a small number of sites were sampled only once, twice or thrice during the 5-year sampling period. Since this is the first survey of macroinvertebrates ever conducted in the LMB and the sampling protocol insured that the collected samples were comparable among sites, these data are important to obtain a first assessment of beta diversity. Therefore, we used median values from data collected on macroinvertebrate and geo-environmental variables to represent each site in our analyzes, as suggested for small sample size by McCluskey and Lalkhen (2007). The community composition data was partitioned into a global macroinvertebrate community data table (including all component communities), and component community data tables (for annelid, crustacean, mollusk and insect communities).

The community composition data were Hellinger-transformed at the beginning of the analyzes (Legendre and Gallagher, 2001;Legendre and Legendre, 2012). For Hellinger-transformed data, the total variance, or total β diversity (BDTotal), of a community composition data table is an index between 0 and 1, and it can be partitioned into local contribution (LCBD) and species contribution (SCBD) indices. An LCBD value is an index showing the degree of uniqueness in taxonomic composition in each site, computed as the relative contribution of a site to BDTotal, so that the LCBD indices sum to 1, whereas an SCBD index shows the relative degree of variation of a taxon across all sites. The BDTotal, LCBD and SCBD indices were computed using the function "beta.div" available in the adespatial package in R (Dray et al., 2016). The Hellinger transformation was used because the corresponding Hellinger distance is one of the dissimilarity functions admissible for beta diversity analyzes (Legendre and De Cáceres, 2013;Legendre and Gallagher, 2001); it does not give high weights to the rare species. To identify significant uniqueness in taxonomic composition of the sampling sites, the LCBD indices were tested for significance against a significance level α = 0.05. The p-values were corrected for multiple testing using the Holm correction to reduce the experimentwise type I error rate of multiple tests. In addition to LCBD, Hellinger-transformed data also allow researchers to compute SCBD indices; this is not allowed by most other admissible dissimilarity functions (Legendre and De Cáceres, 2013;Legendre and Gallagher, 2001). In the following paragraphs, BDTotal, LCBD and SCBD designate the indices of the global macroinvertebrate communities, whereas BDATotal, BDCTotal, BDMTotal, BDITotal, and LCBDA, LCBDC, LCBDM and LCBDI designate the BDTotal and LCBD indices for annelid, crustacean, mollusk and insect communities, respectively. SCBD indices that were higher than the mean of SCBD values identified the taxa that were the most important contributors to BDTotal. Before associating these important taxa with the geo-environmental factors, we normalized these variables using the indications provided by function "boxcoxfit" in the geoR package in R [START_REF] Ribeiro | Package "geoR[END_REF]. The Box-Cox transformation was performed on the geo-environmental variables because this transformation attempts to normalize the variables, thus meeting the assumptions of linear models and residuals' normal distributions [START_REF] Ahola | Effects of access to extra cage constructions including a swimming opportunity on the development of stereotypic behaviour in singly housed juvenile farmed mink (Neovison vison)[END_REF]. Then, we conducted a Redundancy Analysis (RDA, Legendre and Legendre 2012) on the Hellinger-transformed abundance data. To identify which component community was more related to which part (downstream or upstream) of the LMB, we computed Pearson correlations between the richness (number of taxa) and abundance (number of individuals) of the important taxa pertaining to each component community, on the one hand, and to the geographical factors on the other hand.

We independently ran simple and multiple regression analyzes to determine which, among the geographical and environmental variables, mainly accounted for the variation of the LCBD indices. To identify the strength of the regression models, we computed stepwise selection with the Akaike Information Criterion (AIC). The models having the lowest AIC and highest adjusted R 2 were considered to have the strongest influence on the LCBD indices. To investigate the influence of component communities on the LCBD indices, we computed LCBDA, LCBDC, LCBDM and LCBDI, and regressed the LCBD indices (for the global communities) on the LCBD indices of the four component communities. We computed four types of linear regression models: 1) Simple regression models, e.g. LCBD ~ LCBDA; 2) 2component multiple regression, e.g. LCBD ~ LCBDA + LCBDC; 3) 3-component multiple regression, e.g. LCBD ~ LCBDA + LCBDC + LCBDM and 4) all-component multiple regression, LCBD ~ LCBDA + LCBDC + LCBDM + LCBDI. Model selection, based on the AIC, was conducted to obtain a descriptive assessment of the components that contribute most to the variation of the global LCBD indices. Finally, variation partitioning was applied to quantify the variance that main component communities and the significant geo-environmental best explain the variation of the global LCBD indices. All statistical analyses were performed in R (R Core Team, 2013).

Results

General macroinvertebrate composition and environmental variables

In total, 21,810 individuals representing 299 taxa and 90 families were identified in the dataset (see Appendix T2 in Part I: Synthesis). Taxonomic richness was highest at the Mekong delta sites (Fig. 1). Among the taxa, 32 belonged to annelids (2,672 individuals), 38 to crustaceans [START_REF] De Troyer | Water quality assessment of streams and wetlands in a fast growing East African City[END_REF]054), 98 to mollusks (10,603) and 131 to insects (6,481). The most common families of annelids were Naididae (47% of occurrence) and Nereididae (16%); of crustaceans were Palaemonidae (26%) and Corophiidae (16%) and of mollusks were Unionidae (18%), Corbiculidae (14%), Viviparidae (12%) and Stenothyridae (9%). Insect communities were characterized by Diptera (28%), Ephemeroptera (24%), Odonata (17%) and Trichoptera (15%). Three taxa were most widely distributed; two belonged to insects: Ablabesmyia sp. (73% occurrence) and Polypedilum sp. (70%) and one was a mollusk, Corbicula tenuis (67%). In addition to being widely distributed, these 3 taxa were among the top 10 most abundant. Of the total individuals, Ablabesmyia sp. accounted 2.9%, Polypedilum sp. for 3.8%, whereas the 3 most abundant species, Corbicula leviuscula, Limnoperna siamensis and Corbicula tenuis, accounted for 8.4%, 6.1% and 5.8%, respectively. The data on taxonomic richness and abundance, and the environmental variables are summarized in Table 1.

Beta diversity and important taxa -habitat relationship

The total β diversity of macroinvertebrates in the LMB was BDTotal = 0.80, and there was a total of 60 taxa contributing most to the BDTotal. The value of BDTotal is very high, considering that the maximum that can be obtained for Hellinger-transformed data is 1, when all sites have entirely different species compositions. This great variation was also observed for each component community: annelids (BDATotal = 0.72), mollusks BDMTotal = 0.78) and insects (BDITotal = 0.74), excepted for crustaceans (BDCTotal = 0.38). Over all 299 taxa, 60 important taxa had SCBD indices larger than the mean SCBD (0.003), 29 of which belong to insects, 18 to mollusks, 7 to annelids and 6 to crustaceans (see Table 3 Based on the correlation analyzes, the richness and abundance of the important taxa of annelids, crustaceans and mollusks were significantly and negatively correlated with latitude and altitude; the richness and abundance of the important taxa of insects were significantly and positively correlated with latitude and altitude (Table 2). As the Mekong River generally runs from north to south, decreasing latitude and decreasing altitude are both associated with going downstream. The first two axes of the RDA model (Fig. 2a,b) accounted for 15.7% of the total variance of the community of important taxa. Axis 1 (9.4%) of the plot showed the river gradients opposing river depth and river width (left of the plot), which were strongly associated with communities in the Mekong delta, to altitude and latitude, which were associated with communities in the upstream sites. Along axis 2 (6.3%), the communities from some tributaries and main channel sites located upstream in the LMB (Thailand and Laos) were found in the positive part of the axis and associated with high values of dissolved oxygen, while many other tributary sites were found in the lower part of the axis and related to high values of Secchi depth.

Uniqueness in taxonomic composition and its association with geo-environmental factors

Five sampling sites exhibited significant global LCBD indices at the p = 0.05 level after Holm correction for multiple testing. These sites were CM3 (LCBD = 0.023, p = 0.006), TM2 and TM4 (LCBD = 0.021, p = 0.006), and TM1 and VB4 (LCBD = 0.020, p = 0.012) (Fig. 1). CM3 also exhibited significant LCBD indices for all component communities (LCBDA = 0.026, p = 0.032; LCBDC = 0.038, p = 0.006; LCBDM = 0.023, p = 0.006; LCBDI = 0.027, p = 0.006). LCBD indices of the global communities and of component communities, which indicate the uniqueness in taxonomic composition at the sites, are provided in Appendix T3 in Part I: Synthesis. LCBD values are scaled to add up to 1 over the whole study; the mean LCBD value in this study was thus 1/63 = 0.016.

Results of the simple and multiple regressions between the global LCBD indices and geoenvironmental variables are shown in Table 3. Only water conductivity, river depth and Secchi depth were significantly associated with the global LCBD indices, and these variables remained significant after the stepwise selection of the model (Table 3). Water conductivity and river depth were positively associated and accounted for 14% and 12%, and Secchi depth was negatively associated and accounted for 6% of the variation of LCBD indices (Table 3).

Influence of component communities and environmental conditions on global LCBD indices

Among the simple regression models, LCBDI indices were the strongest determinants of the variation of the global LCBD indices. In multiple regressions (2-, 3-, and 4-components), the combination of LCBDI and LCBDM indices best explained the variation of the global LCBD indices (Table 4). The detailed results of the regression models are shown in Table 4. Based on the variation partitioning results (Fig. 3), the variation of the global LCBD indices was well explained by that of the LCBDI and LCBDM (adjusted R 2 = 0.84), whereas the three significant environmental variables explained a smaller fraction (adjusted R 2 = 0.26) of the variation of the global LCBD indices. For the most part, the environmental variable fraction represented explanation shared with LCBDI and LCBDM (adjusted R 2 = 0.24).

Discussion

Macroinvertebrate communities of the LMB

Macroinvertebrate communities in the LMB has been scarcely studied, in particular at a broad large spatial scale. To our knowledge, most of the existing studies were conducted at local sub-watershed scales of the LMB (e.g. [START_REF] Clavier | Spatial and temporal variation of benthic macroinvertebrates in the Nam Gnom Basin receiving discharged waters from the Nam Theun 2 Reservoir (Lao PDR)[END_REF][START_REF] Cuong | Arthropod diversity and abundance in wild rice, Oryza rufipogon, in the Mekong Delta, Vietnam[END_REF][START_REF] Getwongsa | Impacts of agricultural land use on stream benthic macroinvertebrates in tributaries of the Mekong River, northeast Thailand[END_REF][START_REF] Quang | Biodiversity of littoral macroinvertebrates in the Mekong River[END_REF], except for a macroinvertebrate pilot study by [START_REF] Pathoumthong | Macroinvertebrate pilot study for ecological health monitoring in the Lower Mekong Basin[END_REF], which was conducted over thirteen sampling sites of the LMB and reported 218 macroinvertebrate taxa. The number of macroinvertebrate species (299) identified in the present study is the largest ever reported from the basin. Each component (i.e. annelids, crustaceans, mollusks and insects) comprised a higher number of taxa than previous reports of field studies (Table 5). Moreover, we found that 36 macroinvertebrate families were represented by only one species (see Appendix T2 in Part I: Synthesis). This indicates that the LMB could support numerous rare or endemic taxa as found for mollusks that at least 111 species are endemic to the LMB (Köhler et al., 2012).

Table 5 The number of recorded taxa of different taxonomic groups of macroinvertebrates recorded in the present study and in previous reports. The number of arthropod taxa in the parenthesis is the combination of crustacean and insect taxa.

Variation of important taxa and their relationship to habitat characteristics

As stated in our second research question, the annelids, crustaceans and mollusks that have high SCBD indices are more abundant in the downstream part of the LMB, as reflected by the strong negative correlation with latitude and altitude (Table 2). Unsurprisingly, this result also supports previous studies (Arscott et al., 2005;Collier and Lill, 2008;B. Wang et al., 2012). The important taxa included Corbicula leviuscula (code B34 in Fig. 2b), C. lamarckiana (B32), Limnoperna siamensis (B03) and Branchiura sowerbyi (A16), which mostly and abundantly occurred in the downstream part. However, an insect (i.e. Cricotopus sp., I023) was also an important taxon found in the delta. This could be due to the fact that, many species in this genus are capable of withstanding low oxygen concentrations, are resistant to heavy metals, able to withstand high salt concentrations or pollution and can feed on rice [START_REF] Boesel | A review of the genus Cricotopus in Ohio, with a key to adults of species of the northeastern United States (Diptera, Chironomidae)[END_REF][START_REF] Sinclair | Discrimination of Cricotopus species (Diptera: Chironomidae) by DNA barcoding[END_REF], which are all characteristics observed in the delta.

Most of the important taxa characterizing the upstream sites belong to insects. In tropical as well as temperate regions, clear water and high values of dissolved oxygen are mostly found in tributaries and upstream sites, which are mainly preferred by insect taxa (Collier and Lill, 2008;Dobson et al., 2002;Królak and Korycińska, 2008). Of the 3 taxa found with the highest SCBD indices (Polypedilum sp.; Ablabesmyia sp. and Cryptochironomus sp.), Polypedilum sp. (code I033 in Fig. 2b) and Ablabesmyia sp. (I017) were highly associated with sites having high values of SD (Fig. 2b), which were mostly observed in tributaries. Cryptochironomus sp. (I024) and other taxa such as Bezzia sp. (I013) and Anagenesia sp. (I051) were more associated with high values of dissolved oxygen, which occurred at three of the sites with significantly unique taxonomic composition (TM1, TM2, TM4). Surprisingly, Corbicula sp. (B36) and Oligochaeta sp. (A11) were also more associated with these sites. These two taxa may have important taxonomic and ecological value because they were restricted to the main channel shared by Thailand and Laos and its nearby sites.

Beta diversity and uniqueness in community composition

The β diversity of macroinvertebrates in tropical river systems, particularly in South-East Asia, has not been extensively studied (Boyero et al., 2009;[START_REF] Dudgeon | Tropical Stream Ecology[END_REF]. Furthermore, the published papers (Al-Shami et al., 2013;Salmah et al., 2014) did not estimate the β diversity of macroinvertebrates as the total variance (BDTotal) of the communities found at the sampling sites and computed the contributions of individual sampling sites (LCBD indices) to total β diversity. This measure (i.e. BDTotal) quantified "the variation in macroinvertebrate composition among studied sites in the LMB", to which is referred as β diversity by ecologists [START_REF] Anderson | Navigating the multiple meanings of β diversity: a roadmap for the practicing ecologist[END_REF][START_REF] Legendre | Analyzing beta diversity: partitioning the spatial variation of community composition data[END_REF]Legendre and De Cáceres, 2013;[START_REF] Whittaker | Evolution and measurement of species diversity[END_REF][START_REF] Whittaker | Vegetation of the Siskiyou Mountains, Oregon and California[END_REF]. The BDTotal computed here is an independent derived quantity that can certainly measure community differentiation of studied taxa, and thus more suitable to analyse beta diversity of macroinvertebrates in the LMB, when compared to the classical approach (i.e. the additive or multiplicative) which is dependent on alpha and gamma diversity. The great variation of macroinvertebrate composition (BDTotal = 0.80) found may reveal complex evolutionary and ecological processes operating at a site-to-global spatial scale of the LMB.

The contributions of sampling sites (LCBDs) to BDTotal can indicate the ecological uniqueness of each sampling site in terms of community composition and provide valuable information on the level of habitat degradation of sampling sites. These ecological indications can be used to support ecological assessments, restoration and conservation planning of the LMB. For example, we found that sites with large LCBD indices, which are the most different from the centroid of the distribution of the sites in a PCA ordination and hence the most interesting to examine in detail, mostly occurred along the main channel of the LMB. In particular, the 5 sites that had significant uniqueness in species composition (after Holm correction) occurred along the main channel and not in tributaries (Fig. 1a). The discussion on what triggered these sites to have higher degrees of uniqueness in species composition than others is provided in the following paragraphs.

Environmental factors responsible for uniqueness in community composition

Three environmental factors were found to be positively (water conductivity and river depth) or negatively (Secchi depth) associated with the degree of site uniqueness in taxonomic composition (Table 3). These 3 factors collectively explained 26% (adjusted R 2 = 0.26) of the variance in degrees of site uniqueness in taxonomic composition (global LCBD indices). However, our results did not find the types of relationships between β diversity and geoenvironmental factors found in previous studies conducted over smaller areas (Al-Shami et al., 2013;[START_REF] Wang | Patterns of elevational beta diversity in micro-and macroorganisms[END_REF]. This could be due to the dominant effect of anthropic pressure, which is spread along the LMB (Dao et al., 2010;Kudthalang and Thanee, 2010).

Previous papers have shown that conductivity had a positive influence on macroinvertebrate diversity (Lods-Crozet et al., 2001;[START_REF] Rizo-Patrón | Macroinvertebrate communities as bioindicators of water quality in conventional and organic irrigated rice fields in Guanacaste, Costa Rica[END_REF]. However, we found that most of the sites with significant uniqueness in taxonomic composition had low taxonomic richness. High values of conductivity were mostly measured in the main channel sites (e.g. sites TM2, TM4 and the nearby sites) where they receive runoffs and discharge of urban wastewaters from intensified agriculture from surrounding river basins and cities (Dao et al., 2010;Kudthalang and Thanee, 2010;Sor et al., 2017), and consequently lead to high conductivity (Wetzel, 2001). When sources of pollution (i.e. high concentration of inorganic dissolved solids) enter the rivers, only pollution-or disturbance-tolerant taxa (e.g. Oligochaeta, Chironomidae (Diptera) and Gastropoda) can resist (Feld and Hering, 2007;B. Wang et al., 2012). Pollution is the connection between high conductivity and low taxonomic richness.

Deep rivers (i.e. with high values of river depth) that have low Secchi depth can be considered proxies for anthropic activities [START_REF] Baird | Mekong River fish conservation zones in southern Laos: assessing effectiveness using local ecological knowledge[END_REF]Dao et al., 2010), which is why they appear to influence LCBD indices. For example, we found small values of LCBD indices at most of the tributary sites where clear water and low pollution are observed, whereas large values of LCBD indices (e.g. > mean LCBD value) were found at sites with high river depth (e.g. most sites in the delta) and at other sites along the main channel of the upper part of the basin (Fig. 2a, see Appendix T3 in Part I: Synthesis) where high levels of anthropic disturbance were observed. A clear evidence of the association between anthropic activities and high LCBD indices is found at the sites with significant LCBD indices (e.g. CM3, TM1, TM2 and TM4), all of which receive a moderate to high pressure of human impacts. Site CM3 is surrounded by houses, animal wastes and rubbish disposal. Site TM1 seems to be in a very high pressure area since it is opened to many anthropic activities such as animal and human waste disposal, artificial bank creation, local markets, dense population (~10,000 inhabitants), constructions, fishing and boat traffic. Sites TM2 and TM4 are also exposed to waste disposal and fishing, floating houses, tourism (TM2) and agriculture (TM4) (Dao et al., 2010;Kudthalang and Thanee, 2010;Sor et al., 2017). As a result, these sites supported low numbers of taxa (TM1, 9 taxa; TM2, 12; TM4, 21; CM3, 12), which indicate that ecological restoration is needed for these sites and their surroundings.

On the other hand, site VB4 has a significant LCBD value with a moderate number of taxa (34 taxa); the site with highest richness in our study was VB1 (74 species), located close to VB4. VB4 is located at the border between Cambodia and Vietnam and comprises a set of natural land cover (e.g. wood-, shrub-, grass-, inundated and wetland) on the west side of the river. Although the other side of the river has some anthropic activities including houses, fishing and small-scale business, VB4 still had a unique and rich taxonomic composition (Fig. 1a,b). Thus, VB4 and the surrounding sites/areas, particularly on the west side with natural land covers, may have high conservation value. Sites with high LCBD values may have high or low species richness, as shown in Legendre and De Cáceres (2013) and found in the research reported here.

Influence of component communities and environmental conditions on global LCBD indices

The variation partitioning indicated the most striking relationship (adjusted R 2 = 0.84) between the uniqueness in taxonomic composition of the macroinvertebrate communities (global LCBD indices) and the combination of uniqueness in taxonomic compositions of the mollusk and insect communities (LCBDM + LCBDI indices). Note that the global LCBD indices are not simply the sum of the component community LCBD indices; LCBD indices are computed separately for the global study and each component group as the squared distances of the sites to the multivariate ordination centroid. However, the degree of uniqueness in taxonomic composition of macroinvertebrate communities (global LCBD indices) is expected to be contributed by the component communities. In the LMB, LCBDs of mollusk and insect communities, which had a higher total variation (BDMTotal = 0.78 and BDITotal = 0.74, respectively), explained most of the global LCBD variation because these two groups had higher abundances and wider distributions than the annelid and crustacean communities, which had lower total variation (BDATotal = 0.72 and BDCTotal = 0.38, respectively). De'ath (2002) and Davidson et al. (2010) mentioned that taxa with low richness and low occurrence explained less variance of the community composition, and this is similar to our findings for the annelid and crustacean communities.

Co-occurrence among different component communities can directly or indirectly constrain the spatial distributions and the taxonomic abundance of the component communities [START_REF] Miller | Direct and indirect species interactions in an early old-field plant community[END_REF][START_REF] Wootton | Putting the species together: testing the independence of interactions among organisms[END_REF]. Predators (e.g. Odonate taxa and some of the Diptera) may prey upon the taxa of other taxonomic communities, and thus affect the taxonomic occurrence and abundance of the global macroinvertebrate communities. [START_REF] Golfieri | Odonates as indicators of the ecological integrity of the river corridor: development and application of the Odonate River Index (ORI) in northern Italy[END_REF] reported that the abundance of the Odonates is closely linked to the abundance of their prey in the ecosystems. However, the Odonates preferring high water quality in the upstream sites may not directly influence the annelids or crustaceans, most of them being associated with habitats with lower water quality (annelids) or brackish water (crustaceans) in the downstream sites.

For the communities that had a wide distribution in the LMB (e.g. insects and mollusks), their co-occurrence may be the result of niche expansion or competition, and thus they may have indirect interactions by competing for or facilitating resource availability. For example, an increasing topographical complexity of the streambeds, which can alter the near-bed flow, might enhance feeding success of mussel and suspension-feeding caddisfly communities [START_REF] Cardinale | Species diversity enhances ecosystem functioning through interspecific facilitation[END_REF][START_REF] Vaughn | Community and foodweb ecology of freshwater mussels[END_REF].

Several studies suggested that, at large spatial scale, the environment is more important than biotic interactions in governing species composition and distribution (Luoto et al., 2006;[START_REF] Pearson | Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful?[END_REF], while other suggested both (Araújo and Luoto, 2007). Our finding showed that the combination of component communities had a great influence (adjusted R 2 = 0.84) on the degree of uniqueness in macroinvertebrate community composition (global LCBD indices). However, this is restricted to important factors such as precipitation, land use cover, nutrients and sediment loads, which have been reported to better explain the community composition and distribution of macroinvertebrates (Cai et al., 2012;Nicola et al., 2010), but are not available for our study. Nonetheless, our findings further suggest that the combination of biotic and abiotic conditions explain jointly an appreciable amount of the global LCBD variation (adjusted R 2 = 0.24) and seem functionally important together in governing the uniqueness in community composition, and thus influencing beta diversity and composition of macroinvertebrate communities in the LMB.

Conclusion and remarks

The present study revealed the highest number of macroinvertebrate species ever reported from the LMB. The large diversity of different components (annelids, crustaceans, insects and mollusks) led to a great amount of variation, or beta diversity, in overall community composition among studied sites. The important taxa of annelids, crustaceans and mollusks were mostly found in the downstream sites, particularly in the delta, whereas the important taxa of insects were more related to the upstream sites. Most of sites located along the main channels had a high degree of uniqueness in macroinvertebrate taxonomic composition (i.e. high LCBD indices), of which the sites with significant LCBD indices had an exceptionally low richness, which is most likely due to anthropic impacts. An exception was found for one site located in the delta that had a significant LCBD value and moderate macroinvertebrate richness. This is perhaps because of the natural land covers observed on the west side of the river. Mollusk and insect communities, and three environmental variables (water conductivity, river depth and water transparency) were found to be mainly responsible for the variation in LCBD indices.

Our results provide valuable ecological information for selecting locations for conserving different taxonomic groups of macroinvertebrates at broad and small spatial scales. For example, site CM3 and the three other sites (TM1, TM2 and TM4) with significant LCBD indices and low richness are of particular interest for restoration planning, as these locations are experiencing severe degradation of local environments. Site VB4 and the surrounding sites/areas on the west side of the river deserve attention for protection since VB1 had very high richness and VB4 had a significant LCBD index and high richness. The combination of LCBD indices and species richness of the four component communities can thus be used for restoration and conservation planning.

Introduction

Various modelling techniques have been widely implemented in different ecological systems, e.g. terrestrial, freshwater lentic and lotic, and marine ecosystems (Guo et al., 2015;Lek et al., 1996;Lencioni et al., 2007;Park et al., 2003;[START_REF] Schröder | Statistical classification of terrestrial and marine ecosystems for environmental planning[END_REF]. The techniques applied are generally used to investigate or to explain the relationship between the occurrence or abundance of studied species and environmental variables or to predict the relationships being measured (Boets et al., 2013;Goethals et al., 2007). The use of modelling techniques to combine both explaining and predict-ing such relationships is also commonly applied (Call et al., 2016;Roura-Pascual et al., 2009).

The performance of data-driven predictive models is affected by several factors including species prevalence [START_REF] Brotons | Presence-absence versus presence-only modelling methods for predicting bird habitat suitability[END_REF]Hernandez et al., 2006;[START_REF] Stokland | Species distribution modelling-effect of design and sample size of pseudo-absence observations[END_REF]. In most cases, models predicting species which have unequal occupied and unoccupied samples/sites result in a low performance. With a species having a high prevalence, models tend to become better at predicting the presence of that species, and vice-versa for less occurring species [START_REF] Mcpherson | Effects of species' ecology on the accuracy of distribution models[END_REF]. Both cases consequently lead to a low model performance when considering the correct prediction of both the presence and absence of a species. Moreover, it has been demonstrated that species prevalence affects the performance of models in a nonlinear way. For example, Guo et al. (2015) and Manel et al. (2001) reported the nonlinear effect of species prevalence on the performance of models predicting the occurrence of fish and macroinvertebrates. A similar finding has also been revealed for http://dx.doi.org/10.1016/j.ecolmodel.2017.03.006 0304-3800/© 2017 Elsevier B.V. All rights reserved. models predicting the distribution of plants and birds (Allouche et al., 2006;McPherson et al., 2004).

Applications of predictive models have provided knowledge and understanding of the ecology and behaviour of studied taxa, which could support decision making, management and conservation planning. For instance, Chen et al. (2015) used different predictive models as an assessment approach to explain and predict the success of invasive species in China. In addition to the increased use of predictive models, an ensemble modelling framework is recommended when aiming to identify important factors influencing model performance (Araújo and New, 2007). With the ensemble modelling approach, some modelling techniques such as Random Forest and Artificial Neural Networks are found to yield a better predictive performance (Grenouillet et al., 2011;Guo et al., 2015;Segurado and Araujo, 2004). However, although there have been studies assessing the performance of predicting models from an ensemble modelling framework, many have not considered analysing the performance based on a complete prevalence range nor comparing the performance based on different prevalence ranges.

The Lower Mekong Basin (LMB) which is known for its high biodiversity (Sodhi et al., 2004) is a breeding ground of numerous endemic, threatened and endangered species of fish, molluscs and crustaceans (Davidson et al., 2006;[START_REF] Zalinge | It's big, unique and important: fisheries in the Lower Mekong Basin: as seen from a Cambodian perspective[END_REF]. Therefore, it is useful to get more insight into this region based on predictive models which are applicable for different taxonomic groups inhabiting this particular area. To date, the data covering a large spatial scale of the LMB is only available for fish and macroinvertebrates, which were collected by the Mekong River Commission (MRC). The fish data were collected only from the main channel [START_REF] Poulsen | Fish Migration and the Maintenance of Biodiversity in the Mekong River Basin[END_REF], while macroinvertebrates were collected from both the tributaries and the main channel (Dao et al., 2010). In this study, we used the macroinvertebrate data, sampled over 5 successive years (2004)(2005)(2006)(2007)(2008), to build predictive models, which can provide insights on a wide range of keystone species occupying the LMB as well as the neighbouring regions.

The objectives of the present study are to utilize different modelling techniques to 1) predict the occurrence of macroinvertebrate species in the LMB and analyse how the species prevalence (i.e. the occurrence rate of each species in the total samples) affects the behaviour of modelling techniques' performance, and 2) compare the performance of the applied techniques based on the complete prevalence range (i.e. 0.0-1.0), and based on different prevalence ranges (i.e. at a 0.1 interval).

Methods

Data collection and processing

Benthic macroinvertebrates were sampled at 63 sampling sites along the main channel of the LMB and its tributaries by the MRC. This sampling was carried out once a year in March during the dry season from 2004 to 2008. To obtain as much information as possible on macroinvertebrates inhabiting the main river and the tributaries, the MRC collected samples at three locations from the benthic zone of each sampling site: near the left and right banks, and in the middle of the rivers. At each location, a minimum of three samples (where inter-sample variability is low, e.g. tributaries) to a maximum of five samples (where inter-sample variability is higher, e.g. the main channel and the delta) were collected using a Petersen grab sampler. With the grab which has a sampling area of 0.025 m 2 , four sub-samples were taken and pooled to give a single sample covering a total area of 0.1 m 2 . In total, between nine (3 samples × 3 locations) and fifteen (5 samples × 3 locations) pooled samples were collected at each sampling site. Each pooled sample was rinsed using a sieve (0.3 mm mesh size). In the field, samples were sorted and then preserved by adding 10% formaldehyde to obtain a final concentration of about 5%. In the laboratory, the samples were identified to the lowest level possible and counted using a compound microscope (40-1200 magnification) or a dissecting microscope (16-56 magnification). The abundance data of macroinvertebrates per sample (a total area of 0.1 m 2 ) was averaged across all samples (between 9 and 15 samples) collected from each sampling site.

At the sampling site, geographical coordinates and altitude were determined with a GPS (Garmin GPS 12XL). All physical-chemical variables were measured at the three locations where macroinvertebrates were sampled. River width was measured in the field using a Newcon Optik LRB 7 × 50 laser rangefinder, and the river depth was measured using a line metre. Water temperature, dissolved oxygen, pH and water conductivity were measured using a handheld water quality probe (YSI 556MP5). To get a more reliable determination of each variable, the measurement reading was taken at the surface (0.1-0.5 m) and at a depth of 3.5 m or at a maximum depth of the river (wherever less than 3.5 m) and then the average value was recorded for each location. Water transparency was measured with a Secchi disc by lowering it into the water and recording the depth at which it was no longer visible (Dao et al., 2010). The recorded data of each physical-chemical variable was based on the averaged value across the three sampling locations of each site. The distance from the sea was measured by drawing a line from the sea to each locality using GIS-software (ArcGIS version 10.0).

A total of 108 samples were collected from the 63 sampling sites (Fig. 1). Because of unequal sampling efforts (i.e. unequal and different number of samples at each site during the 5-year sampling period) and missing values of environmental variables, we used median values from the collected data to represent each site in our analyses, as suggested by McCluskey and Lalkhen (2007). Therefore, 63 samples remained for the analyses. In total, 299 taxa were obtained from the dataset, of which 131 taxa were insects, 98 were molluscs, 38 were crustaceans and 32 were annelids. The most commonly identified insects belonged to Diptera (37 taxa), Ephemeroptera (32), Odonata (22) and Trichoptera (20). For molluscs, Caenogastropoda (50 taxa), Unionida (18) and Veneroida (15) were represented the most. Most crustaceans belonged to Palaemonidae (10 taxa) and Corophiidae (6 taxa), while most annelids belonged to Naididae (15 taxa) and Nereididae (5 taxa). The detailed information of taxonomic resolution is provided in the Supplementary data Appendix A.

The abundance data of macroinvertebrates from the 63 sites were converted to presence-absence data to analyse how species prevalence (presence/absence) affects the performance of predictive models. Species prevalence was defined as the occurrence rate of each species in the total samples. The species prevalence of a species is an index ranging from 0 to 1, indicating the lowest to highest occurrence rate of that species over all samples. The obtained prevalence values from all macroinvertebrate species formed a complete prevalence range for the present study. For a later analysis, the complete prevalence range was grouped into different ranges based on an interval of 0.1. In other words, the species having a prevalence value between 0.0 and 0.1 were aggregated in a group, and the species having a prevalence between 0.1 and 0.2 were aggregated in another group, and so forth (see Appendix A in Supplementary material).

Predictions

In our predictive models, we used the presence/absence of each species as the response variable. The measured environmental variables used as the input predictors were: altitude, river width, river Article 3 depth, distance from the sea, water temperature, dissolved oxygen, water conductivity, pH and Secchi depth. The summary of all environmental variables is provided in Table 1. Before building the models, the data of these predictor variables were normalized using the zero minimum (Fig. 2). This normalization approach is appro-priate for variables that highly vary (different measuring scales). The normalized values were based on the variables' standard deviation, and the range of all variables remains constant at the same time [START_REF] Jarocka | The impact of the methods of the data normalization on the result of linear ordering[END_REF]. For every prediction, species that occurred only in one instance (a prevalence of ∼0.02) were removed from the data. We removed these species because we used Leave-One-Out cross-validation (LOO), due to our small sample size, to validate the models. With LOO, it is not feasible to split data into training and validation sets for species that have only one occurrence instance.

Model selection and validation

The occurrence of macroinvertebrate species was predicted using four modelling techniques which are commonly used in ecology, namely Logistic Regression (LR), Random Forest (RF), Support Vector Machine (SVM) and Artificial Neural Network (ANN). LR is an extension of the general linear model or so-called generalized linear model (GLM, McCullagh and Nelder, 1989). In LR, a binary response (presence/absence) and a "logit" link function are used. The "logit" link function systematically defines the relationship between the mean value of each response variable and the systematic part of the model, which is a function of the predictor variables. Values between 0 and 1 are the outputs for LR. A threshold value of 0.5 was used to indicate if a species was present or not. If the output value was >0.5, it indicated the occurrence (the presence) of the macroinvertebrate species and vice versa.

RF is a classification method that grows an ensemble of 500 trees. Each tree casts a unit vote for the most popular class according to the input variables. The growing process of RF is by bootstrap aggregating (or bagging), where a tree is randomly grown from the dataset, and this process can give substantial gains in the accuracy of predicting models [START_REF] Breiman | Bagging predictors[END_REF], and it requires no pruning. The random split selection of the bagging process provides nodes of the tree where the split is selected using a subset of the input predictors that are chosen for each tree.

SVM is a supervised machine learning approach that is based on a linear classifier using a maximum hyperplane to separate two classes [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF]. For nonlinear classification, SVM uses kernel functions implicitly to map training data into a higher-dimensional feature space and to compute separating hyperplanes in a way to maximize the margin between classes. The maximum separation hyperplane is defined by a set of support vectors, which are a function of the training data that lie on/the closest to the separating margin. For each class, there is always at least one support vector present. The radial basis kernel function was applied for all SVM models.

ANN is a non-linear statistical data modelling tool used for prediction and classification. ANN is based on the basic function of biological neural networks and on a set of linked computed neurons. One of the most popular approaches of neural networks is a multilayer feedforward neural network (Goethals et al., 2007;Lek et al., 1996;[START_REF] Lek | Artificial neural networks as a tool in ecological modelling, an introduction[END_REF]. This is a backpropagation network that trains the data using a backpropagation algorithm. This network comprises three layers: an input layer, one or more hidden layers and an output layer. Neurons from one layer are in one direction linked to all neurons in the subsequent layers. In this way, connection weights among the neurons are calculated during the training phase. For detailed information on the above mentioned techniques as well as for practical examples of the used techniques, we refer to Van Echelpoel et al. (2015).

To validate the models built using the four modelling techniques, Leave-One-Out cross-validation (LOO) was used. LOO is one of the cross-validation methods using each observation as the validation when the model is trained with the remaining N-1 observations (where N is the total number of observations). It is a time-consuming approach as at a time, each observation is temporarily removed from the dataset, and this process is repeated until every single observation is rotationally left out. LOO has an advantage when N is small (as is the case in this study) because this validation method can provide a nearly unbiased estimation of the accuracy [START_REF] Cawley | Efficient leave-one-out cross-validation of kernel fisher discriminant classifiers[END_REF][START_REF] Efron | Estimating the error rate of a prediction rule: improvement on cross-validation[END_REF].

All modelling procedures were performed in the R language program (R Core Team, 2013). LR was performed using the glm() function of the stats package (R Core Team, 2013), RF using the ran-domForest() function of the randomForest package [START_REF] Breiman | Random forests[END_REF], ANN using the nnet() function of the nnet package [START_REF] Ripley | Package Nnet[END_REF] and SVM using the svm() function of the e1071 package [START_REF] Dimitriadou | Package e1071[END_REF].

Model performance measures and statistical analyses

To evaluate the four modelling techniques, three types of performance measures were used: Cohen's Kappa Statistic (Kappa), area under receiver operating characteristic curve (AUC) and error rate (ER) of the prediction. AUC was calculated as the quantitative principal performance metric. The ER (i.e. misclassification rate) was estimated measuring the fraction of all instances that were not correctly predicted. Each performance measure of the four modelling techniques was regressed against the species prevalence using two types of models: a linear and a quadratic model.

To examine overall performance differences among the four modelling techniques, we performed a Friedman test (nonnormal distribution). Mean performance values and standard errors were compared to identify which technique better performed based on the complete prevalence and based on different prevalence ranges. P-values ≤0.05 were considered to indicate significant differences. The Friedman test was performed using the friedman.test() function of the stats package in R (R Core Team, 2013).

Results

Effect of species prevalence

The species prevalence ranged from 0.03 to 0.73. Linear regression models showed that the prevalence always (for each modelling technique) had a positive and significant effect on the three performance measures used (Table 2). The proportion of variance of Kappa and AUC explained by the linear models ranged from 0.02 to 0.15, while the proportion of variance of error rate (ER) explained ranged from 0.59 to 0.63. When the same data were analysed using a quadratic regression model, the explained proportion of variance of Kappa and AUC increased from 0.06 to 0.31, and the explained proportion of variance of ER increased from 0.66 to 0.84. The coefficient of the quadratic term was always negative and highly significant (p < 0.001, Table 2).

Based on Kappa and AUC, the four modelling techniques had a low performance when predicting species with a prevalence <0.1 and ≥0.6. A high performance was observed when predicting species with an intermediate prevalence between 0.4 and 0.6. Based on ER, the prediction error of the models increased from the lowest to the highest prevalence (Fig. 3).

Model performance

The overall performance among the four modelling techniques was significantly different when considering the three performance measures: Kappa (Friedman chi-squared, F = 12.3; p = 0.006), AUC (F = 11.2, p = 0.01) and ER (F = 350, p < 0.001). The highest mean Kappa (0.19) and mean AUC (0.60) were obtained for ANN, followed by LR (mean Kappa: 0.16, mean AUC: 0.59) and RF (mean Kappa: 0.12, mean AUC: 0.55), while SVM yielded the lowest mean Kappa (0.06) and mean AUC (0.53). On the other hand, a lower mean ER (0.09) was obtained for RF and SVM, while ANN and LR had a higher mean ER (0.13 and 0.16, respectively).

Based on Kappa and AUC, the performance of the models varied for different prevalence ranges (Fig. 4). ANN and LR performed better than RF and SVM for the prevalence range <0.1; ANN, RF and LR performed better than SVM for the prevalence range 0.1-0.2; ANN and RF performed better than LR and SVM for the prevalence range 0.2-0.3, but not significantly different according to the standard error. Based on ER, RF and SVM performed better than ANN and LR in predicting species with a prevalence range <0.1, between 0.1 and 0.2 and between 0.2 and 0.3. Based on all calculations, the model performance was not significantly different for the prevalence range ≥0.3 (Fig. 4).

Discussion

The effect of species prevalence

We found a highly significant quadratic effect of species prevalence on the performance of the four modelling techniques for all three measures, which clarifies the nonlinear relationship found

Table 2

Results of linear regression models (y = ˛ + ˇ1x) and quadratic regression models (y = ˛ + ˇ1x + ˇ2x 2 ) for the effects of the prevalence of macroinvertebrate species on the three performance measures (Kappa, AUC and ER) in previous studies. Species prevalence has been reported to influence model performance in a nonlinear way. This finding has been demonstrated from the prediction of freshwater aquatic species, e.g. fish and macroinvertebrates (Guo et al., 2015;Manel et al., 2001) and of terrestrial plant and bird species (Allouche et al., 2006;McPherson et al., 2004). Moreover, our results also indicated that the performance of predictive models responds in a positive linear way as long as the prevalence remains below 0.2 for Kappa and AUC or below 0.3 for ER (Fig. 3). At a higher prevalence than the intermediate range (e.g. 0.4-0.6), the performance of the models, based on Kappa and AUC, showed a negative relationship. Although we did not have many samples with a high prevalence to support this assumption, evidence that a higher prevalence than the intermediate negatively affects the performance was found from previous studies (e.g. [START_REF] Guisan | Predicting reptile distributions at the mesoscale: relation to climate and topography[END_REF][START_REF] Mcpherson | Effects of species' ecology on the accuracy of distribution models[END_REF]Stockwell and Peterson, 2002). However, the negative relation between the higher prevalence and the ER was not observed (Fig. 3). Nevertheless, when considering the whole prevalence range for all the three performance measures, our findings suggest rather a highly significant quadratic than a simple linear effect, as reflected by a higher proportion of variance explained by the quadratic regression models (Table 2).

Quadratic effects of species prevalence have been shown to influence Kappa and ER, but not AUC (Allouche et al., 2006;Manel et al., 2001). Our findings also indicate a quadratic effect of species prevalence on AUC. This notable difference may be due to the input predictor variables, and the type of response variables (e.g. fauna vs. flora). Geomorphological variables (e.g. river depth, river width and distance from the sea) and water chemistry (e.g. dissolved oxy-gen, water conductivity and pH) have shown to strongly influence the occurrence and distribution of macroinvertebrates (Chadwick et al., 2006;Rizo-Patrón et al., 2013), whereas climatic variables and vegetation density are known to have important effects on plant and bird distributions, respectively [START_REF] Amissah | Rainfall and temperature affect tree species distribution in Ghana[END_REF][START_REF] Toledo | Distribution patterns of tropical woody species in response to climatic and edaphic gradients[END_REF]. For example, Allouche et al. (2006) found no significant quadratic effect of prevalence on AUC when modelling the distribution of woody plants using three climatic variables (i.e. a mean value of annual rainfall, of daily temperature and of minimum temperature). On the other hand, McPherson et al. (2004) found to some extent a quadratic dependency of AUC on the prevalence when using land surface and air temperature, vapour pressure deficit and normalized difference in vegetation index to predict the distribution of birds. This indicates that different predictor variables used to predict the same type of response variable or vice-versa are more likely to yield a different model performance, and thus affect the relationship between species prevalence and model performance when they are regressed against each other.

For all models, the best performance (i.e. Kappa and AUC) was found for species with an intermediate prevalence (e.g. 0.4-0.6, Figs. 3 and4). In between this prevalence range, there is a smaller bias for the models to select presence/absence data for training and validation sets. A more or less equal distribution between the presence and absence, which are likely to result in many correctly predicted instances of both the true positive and true negative fractions, could be the main reason responsible for the best predictive performance (Allouche et al., 2006;Manel et al., 2001).

For very common (prevalence ≥0.6) and for rare species (prevalence <0.1, see Appendix A in Supplementary material), the models yielded a low performance (i.e. Kappa and AUC). In these cases, there is a high imbalance between the presence and absence data. Therefore, the models can correctly predict many true positive instances and a few or perhaps no true negative instances for common species, and vice versa for rare species. This explains the low performance of the models because Kappa and AUC are designed to reflect model performance in absence and presence instances simultaneously (Cicchetti and Feinstein, 1990;Cohen, 1960 Zweig andCampbell, 1993). Thus, a few or no instances of either the true positive or the true negative fraction results in a low performance. On the other hand, opposite results were found for rare species when the performance was based on the error rate (ER); the predicting models had a highly reliable performance (a low prediction error) (Figs. 3 and4). Based on previous research, ER is considered to be a misleading measure (Fielding and Bell, 1997;Manel et al., 1999Manel et al., , 2001)). One main reason is that ER only takes into account misclassified instances to estimate the error rate. Indeed in our results, the ER does not reflect the explanation that rare species are well predicted. This suggests that ER is giving an ambiguous performance of models for predicting the presence of species with a lower occurrence or of rare species.

Comparison of modelling techniques' performance

Overall, based on Kappa and AUC, we found that ANN performed the best across the complete prevalence range, whereas SVM had the lowest reliability, and RF and LR had an intermediate performance. While some authors (e.g. Mastrorillo et al., 1997;Pearson et al., 2002;Segurado and Araujo, 2004) agree that ANN provides advantages over other techniques for predicting species occurrence, others found RF to be better (Gallardo and Aldridge, 2013;Grenouillet et al., 2011). Studies comparing model performance usually suggest different modelling techniques based on the performance measures. For example, Guo et al. (2015) compared the performance of nine modelling techniques and found no significant difference, but suggested RF as a better technique, especially in terms of interpretation. Whereas Segurado and Araujo (2004) compared nine modelling techniques' performance and found a significant difference and suggested ANN to perform better. However, both studies may not be fully comparable due to different criteria in selecting different environmental predictors.

On the other hand, our results, based on the overall mean ER, indicated that SVM and RF performed the best, as reflected by the lowest mean error rate (0.09 ± 0.01), whereas LR had the highest error rate (0.19 ± 0.01). The ER was highly dependent on the species prevalence, as implied by a highly explained proportion of variance of the ER (Table 2). Therefore, our findings, along with the recommendation from previous studies (e.g. Fielding and Bell, 1997;Manel et al., 2001Manel et al., , 1999)), suggest that ER is highly biased, compared to Kappa and AUC, to evaluate model performance.

The best model in terms of performance may not be the best model in terms of applicability. The performance is more likely related to theoretical knowledge of how each model works and to the parameterization settings (Araújo and Guisan, 2006;Elith and Graham, 2009). Eliminating less important variables or increasing the number of tree nodes when fitting a model may result in a higher model performance. However, these settings have to be carefully taken into account because they may lead to an overfitted model, which produces a less general result with a lower applicability to different situations (Babyak, 2004). Model applicability is dependent on the type of data (e.g. missing values and data distribution, De'ath and Fabricius, 2000;Therneau and Atkinson, 1997) and on where the data is derived (e.g. different ecological regions, Guisan and Thuiller, 2005;Randin et al., 2006). Most of the studies evaluate and apply modelling techniques based on the data collected within particular regions, but the applicability of those techniques to a different geographical range is hardly assessed (Fielding and Haworth, 1995;Kleyer, 2002;Özesmi and Mitsch, 1997). As such, our findings may reveal a unique behaviour of each modelling technique and its applicability for the Lower Mekong Basin and the neighbouring areas rather than for other geographical areas. This is because the performance of the same technique applied in different geographical regions can be variable (Randin et al., 2006). Accordingly, it is important to balance the model performance and its applicability when comparing modelling techniques.

Furthermore, each technique may have its own distinct characteristics regarding adjustment to the response variables [START_REF] Guisan | Predictive habitat distribution models in ecology[END_REF]. Based on Kappa and AUC, our study found that ANN and LR were more suitable to predict the occurrence of rare species (e.g. prevalence range <0.1); although the performance was not very high, it was still higher than the performance of RF and SVM. Rare species mostly occur in a particular geographical region and prefer a specific set of environmental conditions [START_REF] Prendergast | Rare species, the coincidence of diversity hotspots and cons[END_REF]. For example, Ephemeroptera (15 taxa which made up 31% of all insect species and 11% of all species that have a prevalence range <0.1) and Trichoptera (5 taxa, see Appendix A), two groups of sensitive species, mostly occur in mountainous areas with clean or unpolluted environments (Md Rawi et al., 2013;Suhaila and Che Salmah, 2014). Some species of molluscs (e.g. Gyraulus and Thiara) have been recorded from restricted habitats [START_REF] Choubisa | Freshwater snails (Mollusca: Gastropoda) as bio-indicators for diverse ecological aquatic habitats[END_REF]. Many annelid species in the family Naididae and crustacean species in the family Palaemonidae have been mainly found from a more polluted environment and from estuaries or brackish water, respectively (De Grave et al., 2008;[START_REF] Martins | Tubificidae (Annelida: Oligochaeta) as an indicator of water quality in an urban stream in Southeast Brazil[END_REF]. Perhaps, these site-specific environmental conditions could be well predicted by models using ANN and LR. The performance of the same modelling techniques (ANN and LR) and of RF were comparable when predicting species with a prevalence range between 0.1 and 0.2 (Fig. 4). For species having a wider occupancy (e.g. dipteran insects and a mollusc species Corbicula tenuis), the performance among the modelling techniques was not significantly different. This could be due to the fact that those widespread species (e.g. species with an intermediate prevalence or higher) are able to respond to a wide range of environmental conditions, which could be predicted by each type of technique used. Therefore, our findings suggest that species prevalence should be carefully taken into account when assessing model performance and predicting species occurrences. Selecting modelling techniques to predict a given species should also depend on the characteristics of the data and the purpose of the prediction.

Conclusion and remarks

Overall, we found a highly significant quadratic effect of species prevalence on the performance of the four modelling techniques. Compared to error rate (ER), the dependency of Kappa and AUC on the species prevalence was rather low, making them a better measure of model performance. A maximum performance was obtained when the species prevalence range was situated between 0.4 and 0.6. ANN generally provided a better overall performance than other modelling techniques and yielded a higher reliability when predicting species with a low occurrence. Our findings could offer useful knowledge regarding the understanding and possible proposition of modelling techniques for other species of interest that inhabit the Lower Mekong Basin and the neighbouring areas.

This study provides clear insights into the application of different modelling techniques when predicting species occurrence and how their performances vary for different species prevalence ranges. Each modelling technique has its strengths and weaknesses. Thus the selection of an appropriate technique should depend on data availability and the purpose of the study, and should balance between model performance and applicability. Modelers may consider a technique (e.g. ANN) that is seen to be generally robust across all species and across species with a small distribution range or a low occurrence (e.g. prevalence range <0.1). Environmental predictors and species prevalence (as shown in this paper), should be carefully taken into account for studies attempting to assess the distribution of a given species.

Introduction

Invasive species have become a major concern for the global economy and environment (Sala et al. 2000). A large proportion of the economy has been spent on the management of agriculture, grassland and various natural ecosystems to mitigate the effects of alien invasive species [START_REF] Williams | The economic cost of invasive non-native species on Great Britain[END_REF][START_REF] Hulme | Weed risk assessment: a way forward or a waste of time?[END_REF]). Moreover, their spread threatens native species of the same taxonomic groups and surrounding biotic communities via e.g. species replacement, food web reorganization and community composition simplification (Gurevitch and Padilla 2004;Bernauer and Jansen 2006;Didham et al. 2007). Once invasive species have successfully colonized new habitats and co-exist with native species, eradication is rarely possible (Regan et al. 2006). Consequently, the rate of replacement of native species by invasive species increases, which can thus lead to an overall decrease of native species [START_REF] Olden | Ecological and evolutionary consequences of biotic homogenization[END_REF].

Invasion success depends on traits of the invaders and the suitability of invaded environments [START_REF] Kolar | Progress in invasion biology: predicting invaders[END_REF]. Some taxa, e.g. Corbicula spp., are highly successful invaders due to their rapid Article 4

spreading ability and their capability to withstand a wide range of environmental conditions (Werner and Rothhaupt 2007;Pigneur et al. 2014). Habitat modifications, resulting in changed physical and chemical conditions, often promote the local abundance and regional distribution of alien species (Didham et al. 2007). Increased trade (shipping) and improved chemical water quality may also promote the number and abundance of alien species [START_REF] Iksr | Das Makrozoobenthos des Rheins[END_REF]Boets et al. 2016). Therefore, identifying the environmental conditions in which alien species solely exist or co-exist with native species and determining those locations that could be invaded in the future will provide essential knowledge to support environmental management and conservation planning.

River systems in Europe have been exposed to the introduction of alien macroinvertebrate species. In the river Rhine, for example, alien species contribute 11.3% of the total macroinvertebrate species richness (Leuven et al. 2009). Among the macroinvertebrate invaders, molluscs constitute a large proportion (Leuven et al. 2009;Nunes et al. 2015). However, for most river systems in Europe, e.g. river systems in Flanders, the environmental conditions in which only alien species occur or the conditions preferred by both alien and native species (co-occurrence) are poorly studied. Recently, an inventory and habitat suitability model of alien macrocrustaceans in Flanders was conducted (Boets et al. 2013;Boets et al. 2016). Moreover, Boets et al. (2016) reported that alien mollusc species, e.g. the New Zealand mud snail (Potamopyrgus antipodarum J.E.Gray, 1843) and the acute bladder snail (Physella acuta (Draparnaud, 1805)), are highly abundant in the river systems of Flanders. As such, there is an urgent need to gain insight into the environmental conditions preferred by alien molluscs and to determine the conditions that allow alien molluscs to co-occur with native molluscs, as a basis for invasion control (e.g. locations and type of actions which deserve priority).

The aim of our study is to 1) provide an analysis of the spatio-temporal occurrence of alien molluscs and their co-occurrence with native molluscs in the river systems of Flanders over the past two decades (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), and 2) identify key determining physical-chemical variables associated with the sole occurrence of alien molluscs and their cooccurrence with native molluscs. We expected that a combination of some key variables of water quality could provide suitable conditions promoting alien molluscs to occur and to co-occur with native molluscs.

Material and methods

Data collection and treatment

The Flemish Environment Agency (VMM) has collected biological and environmental data in Flanders since 1989. The samples have been collected at more than 2500 sites spread over different water bodies. In this monitoring program, a standard handnet was used to collect macroinvertebrates following the method described by Gabriels et al. (2010). At sampling sites where the kick sampling method was not possible, artificial substrates were used. Seven alien and 27 native mollusc genera were identified. Electrical conductivity (EC), pH and dissolved oxygen (DO) were measured in the field with a hand-held probe (Cond 315i, oxi 330, wtw, Germany and 826 pH mobile, Metrohm, Switzerland). All additional chemical variables, i.e. ammonium (NH 4 + ), chemical oxygen demand (COD), biological oxygen demand (BOD), total phosphorus (Pt), nitrate (NO 3 -), nitrite (NO 2 -), Kjeldahl nitrogen (KjN), and orthophosphate (oPO4), were retrieved from the monitoring dataset compiled by the VMM and which is accessible online (www.vmm.be). Nutrient analysis was performed spectrophotometrically in accordance with ISO 17025. GIS software (version 9.3.1) applied to the Flemish Hydrographic Atlas was used to determine the slope and the sinuosity of a watercourse using the difference in height in between two points 1000 m apart, and on a stretch of 100 m, respectively. For further detailed information on the determination of physicalchemical variables, we refer to Boets et al. (2016).

Data from 1991 to 2010 were used for the analyses. The data were divided into 4 periods, each encompassing 5 years of sampling effort (i.e. 1991 to 1995, 1996 to 2000, 2001 to 2005 and 2006 to 2010). This division provided more samples with which to model the preferred environmental conditions, and can provide useful information on changes in the occurrence of alien molluscs and their cooccurrence with native molluscs for each period. Due to limited frequency of occurrence for most alien mollusc genera (Table 1, 2), we decided to merge all alien genera to form one categorical variable. This provided us with a higher number of instances for our predictive models and thus a better and more robust development of the model. Moreover, we were not aiming to make predictions for individual taxa but rather to reveal common environmental conditions that most of the alien molluscs prefer. In the same way, all native genera were also merged to form one categorical variable. Environmental preferences of each genus of alien and Article 4

Table 1. The occurrence instances of each alien mollusc genus and of all genera that are merged together, compared to all collected samples (7695 samples) within the studied period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). The sum of occurrences of different alien genera is 2522 instances, and the overlapping occurrences of alien genera are 494 instances. 1991-1995 1996-2000 2001-2005 2006-2010 1991-1995 1996-2000 2001-2005 native molluscs are provided in the Supplementary material Appendix 1. Each sampling site for each period was then categorized as: a "native" site (i.e. a site having only native molluscs present), an "alien" site (i.e. a site having only alien molluscs present) and a "co-occurrence" site (i.e. a site having both alien and native molluscs present), and this status was used as the response variable in analyses. Physical and chemical water quality variables that had missing values for more than 5% of the total samples were removed from the analyses. Each period thus consisted of one response categorical variable (native/alien/co-occurrence) and 13 predictor variables. To visualize the occurrence of alien molluscs and their co-occurrence with native molluscs, we produced an occurrence map using GISsoftware (ArcGIS version 9.3.1). The summary of the physical and chemical water quality variables and of the response variable is shown in Table 3.

Genus

Article 4

Modelling

A Classification Tree (CT) model was used to predict the occurrence of native and alien molluscs and their co-occurrence, and to identify the determining physical-chemical variables. The CT was chosen among other machine learning approaches according to its performance for both predictive power and the importance of input variables (Chen et al. 2015). Moreover, this decision tree model is relatively simple to implement, easy to interpret, and it tolerates missing values during both the training and testing phases (Therneau and Atkinson 1997;De'ath and Fabricius 2000).

A CT model is based on growing and pruning to select an optimal tree. In the growing phase, the decision trees were fitted using a recursive partitioning algorithm. In each growing phase, the root of a tree (the initial node) is built from the most informative input variables. From the root, the data is split into left and right branches based on the splitting rules defined by the values of selected input variables. The growing process continues up to the terminal node until all the data in that node are of the same class or until some other stopping criterion is reached. The terminal nodes are called leaves and are labelled with the corresponding class (Quinlan 1986). In the pruning phase, the tree was pruned by setting a complex parameter at cp=0.05. Where the tree had only a root, we decreased the cp to a level (e.g. cp = 0.04, 0.03, 0.02, 0.01) that at least two terminal nodes were produced. In practice, the first few splits mostly provide a very informative division of the data (Therneau and Atkinson 1997). These criteria were set in order to make the trees easily interpretable and comparable in terms of the number of variables and complexity.

For each period, a three-fold cross-validation was used to train and validate the models. To build reliable models and to avoid misidentifying the key variables determining each class of the response variable, we made 3 replications of the three-fold cross validation. For each 3-fold cross-validation, the data was shuffled and randomly split into three subsets; two subsets were used for training and one subset for validation. For the second and third replication, we re-shuffled the data and randomly split it into new training and validation sets following the same procedures. From each training and validation set, a model was built and in this way, a performance value and the importance of each variable (as a percentage) of nine different models (3 models of each three-fold cross validation × 3 replications) were calculated. A mean performance value, obtained from the nine models, was used as a final criterion for model evaluation. Cohen's Kappa Statistic (K) and Correctly Classified Instances (CCI) were used to evaluate the model performance. The higher the value of K (ranging from 0 to 1) and of CCI (ranging from 0 to 100), the better the model predicts the response variable. The importance of each variable determining the preferred environmental conditions of each class (native/alien/co-occurrence) was averaged across the nine models. To identify which variables significantly determine the preferred conditions of each class, the importance of each variable was compared based on the standard error. The same procedures and criteria were applied for modelling the data for each period across the whole studied period.

To develop and validate the models, and to construct the classification trees (for visualization), the package "rpart" in R [START_REF] Breiman | Classification and regression trees[END_REF] was used. As we had many models which produced many trees for each period, we chose to construct the tree based on all data points of each period to be the representative one. All statistical analyses and calculations (K and CCI) were performed in R (R Core Team 2013).

Results

Occurrence of alien molluscs and their co-occurrence with natives

Overall, the occurrence of alien molluscs spatially and temporally increased in fluvial systems in Flanders (Figure 1). The "alien" sites accounted for 0.6%, 0.7%, 2.3% and 7.2% of the total samples for the period 1991-1995, 1996-2000, 2001-2005 and 2006-2010, respectively. The alien taxa which showed a notable increase in occurrences include Corbicula, Dreissena, Ferrisia and Physella. In the last period, a new alien genus (Menetus) was also recorded.

Article 4

Detailed information on the occurrences and abundance of each alien mollusc genus is provided in Table 2. There was also a spectacular increase in the "cooccurrence" sites; they respectively accounted for 7.1% (66/935 samples), 7.2% (146/2020 samples), 21.1% (567/2689 samples) and 47.3% (875/1850 samples) for the corresponding periods (Table 3).

Key determining variables

For each period, the CT models were only able to reliably predict the "co-occurrence" and the "native" sites. The representative trees of the CT models for each period are shown in Figure 2. Sinuosity was always one of the most important variables for the models of each period and this, together with chemical water quality variables (e.g. NH 4

+ , NO 3 -, COD, pH, Figure 2) and EC, were the key factors determining the predictive models (Figure 2, 3). When sinuosity was lower than 1.01, co-occurrence between alien and native molluscs was evident where NH 4 + <0.4 mg/L and NO 3 -≥2.4 mg/L for the period 1991-1995; COD <18.9 mg/L and NO 3 -in between 2.6-4.7 mg/L for 1996-2000; COD <21.8 mg/L and NO 3 -in between 2.3-5.1 mg/L for 2001-2005; and where pH ≥7.2 for 2006-2010 (Figure 2). Sites where only alien molluscs occurred could not be reliably predicted based on the available physical-chemical data (Figure 4).

Model performance

The mean Cohen's Kappa Statistic (K) and the mean overall Correctly Classified Instances (Overall-CCI) of the models decreased from the past to the most recent period (Figure 5A-B). For the four periods (1991-1995, 1996-2000, 2001-2005 and 2006-2010), the corresponding mean K was 0.34, 0.32, 0.22 and 0.16, and the corresponding mean Overall-CCI was 93%, 92%, 79% and 54%. The mean CCI of models predicting the "co-occurrence" sites (CCI-co-occurrence) was lower for the first 3 periods (28%, 25% and 23%), while it was higher for the most recent period (52%, Figure 5C). On the contrary, the mean CCI of models predicting the "native" sites (CCI-Native) for the past 3 periods (99%, 98%, and 96%) was substantially higher compared to that of the most recent period (66%, Figure 5D). The models predicting the "alien" sites for the 4 periods did not yield any reliable prediction. Only one model that was based on the data from the latest period correctly predicted one instance of "alien" sites. The confusion matrices obtained from the models and which were used to calculate the model performance measures (K and CCI) are provided in Appendix 2.

Discussion

Occurrence of alien molluscs and their cooccurrence with the natives

Alien molluscs have spread spectacularly in several European river systems during the last few decades, e.g. the rivers Rhine and Meuse (Bernauer and Jansen 2006;[START_REF] Collas | Effects of desiccation on native and non-native molluscs in rivers[END_REF]Pigneur et al. 2014). This phenomenon is similarly observed in our study in the rivers in Flanders. The remarkable increase in the occurrence of alien molluscs over the past two decades may reveal their outbreak or invasion success across the fluvial systems in Flanders, which resulted in an increased incidence of co-occurrence of alien and native molluscs.

Key determining variables

CT models indicated that co-occurrence is mainly determined by sinuosity and by a set of chemical water quality variables (i.e. NH 4 + , NO 3 -, COD, pH and EC). Sinuosity was always one of the most important factors determining co-occurrence as it formed the main root in all models. Sites having a low sinuosity (<1.01), which corresponds to mainly straight rivers, may be subjected to a high number of passing ships, which is considered one of the main pathways of invasions (Bij de Vaate et al. 2002;Nunes et al. 2015). Straight rivers shorten travelling distances, resulting in more frequent transportation, thus allowing a large and frequent amount of ballast water to be released. Consequently, with a higher number of introductions, the survival rate of alien molluscs increases (Gollasch 2006). We found that the hotspots of mollusc invasion (the alien and the co-occurrence sites) were mainly situated in brackish polder watercourses and in large rivers that have a shorter distance to the ports in the Rhine delta, to the coast of the North Sea and to other large rivers (e.g. Meuse River). This observation is supported by evidence in Boets et al. (2016) and in Grabowski et al. (2009) who found that alien fauna mostly inhabited large rivers where intensive navigation takes place.

Besides river morphology, water quality is one of the major factors influencing the distribution and diversity of freshwater fauna (Leuven et al. 2009;Wang et al. 2012). Key variables used to evaluate water quality are NH 4 + , NO 3 + , COD, pH, DO, Chloride, and total phosphorus (US-EPA 1986;[START_REF] Seq-Eau | Système d'évaluation de la qualité de l'eau des cours d'eau[END_REF]WWF 2007;Chea et al. 2016). Good water quality supports a high diversity of invertebrates (Leuven et al. 2009). In our study, where sampling sites had a low sinuosity (<1.01), cooccurrence was mainly dependent on chemical water quality status.

River systems in Europe as well as in Flanders have suffered from severe water quality degradation in previous decades. During these periods, some native species were reported to disappear (Bernauer and Jansen 2006) and only species that were able to withstand this water quality degradation remained. In the early 1990s, the invasion of alien molluscs seemed to be at an initial stage as reflected by the presence of only five alien taxa with limited spatial occurrence and low abundance (Figure 1, Table 2). Although water quality was degraded during these periods, we found that a low NH 4 + concentration (<0.4 mg/L) and a nitrate concentration higher than 2.4 mg/L were preferred by both the alien and remaining native molluscs. This could be because higher NH 4 + concentrations negatively affect their living conditions (Friberg et al. 2010). Moreover, an enrichment of nutrient content (e.g. NO 3 ->2.4 mg/L) could be advantageous for alien molluscs at the beginning of the invasion. Nutrient enrichment provides resources that can be used by alien molluscs and thus may enhance their proliferation [START_REF] Hall | Exotic snails dominate nitrogen and carbon cycling in a highly productive system[END_REF]Strayer 2010). However, a very high nutrient concentration (e.g. NO 3 ->25 mg/L) is also an indication of a high level of water pollution (SEQ-Eau 2003), which can negatively affect both alien and native species (Boets et al. 2013). Nonetheless, a low NH 4 + concentration and an optimum nutrient content (i.e. mean NO 3 -of 3.3 mg/L) facilitated alien and native molluscs co-existence during the early 1990s.

The late 1990s and early 2000s appeared to be an expansion phase for some alien molluscs (e.g. Ferrisia, Physella, Dreissena and Corbicula) since we found a large increase in their occurrences and abundance (Table 2). This increase may result in a substantial effect on water quality in river systems. A high abundance of filter-feeders (e.g. Dreissena and Corbicula) and Physella, which also feeds on phytoplankton, can lead to increased nutrient concentration. This is because the filter feeders consume phytoplankton, thus limiting the abundance of nutrient utilizing phytoplankton communities, allowing nutrient inputs (e.g. nitrate and phosphate) from surrounding areas to increase in river systems [START_REF] Lavrentyev | Effects of the zebra mussel on nitrogen dynamics and the microbial community at the sediment-water interface[END_REF]Pigneur et al. 2014). However, this might not have led to a severe impact on water quality because during these periods there was a successful rehabilitation program for improving water quality, restoring riverine ecosystems and improving habitat connectivity in Europe (Leuven et al. 2009). A decreasing trend of nitrate and other water quality variables (e.g. COD) was observed in the river systems in Flanders (UN 2004). Therefore, it can be inferred that improved water quality and the rehabilitation programs not only helped to recover the diversity of native species but also promoted the occurrence and abundance of alien species. Indeed, Leuven et al. (2009) indicated that when hydromorphological conditions remain unchanged, improvement in water quality promotes alien species.

In the late 2000s, a high percentage of surface water bodies (43%) in Europe were considered to be at a "good status". The number of waterbodies increased to 53% in 2015, and Flanders was one of the regions that well implemented the policy of the Water Framework Directive (EU 2015). The improvement in water quality was associated with a great increase in the spatial occurrence and abundance of alien molluscs. The two filter feeders (i.e. Corbicula and Dreissena) always expanded their range and frequency in the late 2000s even though the number of sampling sites and total sample size were lower than the early 2000s. Moreover, a new alien genus, Menetus, also emerged in the late 2000s. This suggests that the late 2000s can be considered as the expansion phase of the existing alien molluscs and the beginning of the expansion phase of the recently introduced alien genus, Menetus. These results are therefore unlikely to be an effect of the sampling strategy, but rather reflect the suitable physicalchemical conditions for alien molluscs to spread and proliferate. However, we found that the pH value, which was the second most important variable in determining the co-occurrence of alien and native molluscs (Figure 2, 3) during this period, was relatively high (pH: 6.5-8.5). This is probably related to certain specific environmental conditions linked to geographic regions. For example, a high pH value could be mainly recorded from brackish polder watercourses and from the main harbour watercourses where a high level of seawater intrusion occurs and a high intensity of human-related activities takes place. The alkaline watercourses may have a higher concentration of calcium and magnesium compared to inland watercourses, thus may be preferred by some alien molluscs (e.g. Ferrisia and Physella) since they require a large amount of calcium and magnesium to form their shells [START_REF] Brodersen | The effect of calcium concentration on the crushing resistance, weight and size of Biomphalaria sudanica (Gastropoda: Planorbidae)[END_REF]Gallardo and Aldridge 2013). Moreover, molluscs in the family Physidae, including the alien Physella and the native Physa, are pollution tolerant and can occur in areas with high pH values (Rodrigues Capítulo et al. 2001;De Troyer et al. 2016). Boets et al. (2013) also found an increase of alien species abundance with increasing pH.

Across the models of the four periods, other than NH 4 + , NO 3 -, COD and pH, EC was also selected as one of the key determining variables. Many studies have shown that high conductivity mostly favours alien species. For example, alien amphipods in the Vistula and Oder rivers of Poland, alien macrocrustaceans in fluvial systems in Flanders and alien gastropods in isolated ponds in Poland all benefited from high conductivity (Grabowski et al. 2009;Boets et al. 2013;Gallardo and Aldridge 2013;[START_REF] Spyra | Identifying factors linked to the occurrence of alien gastropods in isolated woodland water bodies[END_REF]. It is likely that alien mollusc species are well adapted to withstand high values of EC, while only a few native species (e.g. Lymnaea stagnalis (Linnaeus, 1758)), which occur in a wide range of environmental conditions [START_REF] Brown | Small-scale to large-scale and back: larval trematodes in Lymnaea stagnalis and Planorbarius corneus in Central Europe[END_REF], are able to co-exist in these areas. Our findings suggest that increased spatial co-occurrence results from the introduction or migration of alien species to the connected environments where native species are present or to new environments where Article 4 conditions fit them best, e.g. sites having a high conductivity. Colonizing new niches where most of the native species are not able to thrive is one of the main strategies found among invaders [START_REF] Verbrugge | Sensitivity of native and non-native mollusc species to changing river water temperature and salinity[END_REF]. This might further imply that, after success in co-existing with the native species, the aliens may overrun some native species and further spread to new areas. Clear evidence can be seen from the drastically increased occurrences and abundance of most alien molluscs from the past to the present, from the decrease in occurrences and abundance of the native Physa, Aplexa and Theodoxus, and from the disappearance of the native Bythinella, Margaritifera, Marstoniopsis and Myxas (Table 2).

Although not directly taken into account in our models, previous studies have demonstrated that dispersal vectors are important in making predictions on future locations that may be invaded by alien species. Indeed, recent research on the dispersal of alien macrocrustaceans in Flanders (Boets et al. 2013) showed that increased shipping and the connection between waterways promote the dispersal of alien species. Moreover, habitat conditions (bank structures and substrates) and hydrological variables (e.g. distance to ports/coast, flow regime and connectivity) can also influence the occurrence of alien species [START_REF] Josens | Native and exotic Amphipoda and other Peracarida in the River Meuse: new assemblages emerge from a fast changing fauna[END_REF][START_REF] Messiaen | Alien macrocrustaceans in freshwater ecosystems in the eastern part of Flanders (Belgium)[END_REF]). However, further research suggested that although habitat and hydrological variables can improve model reliability when predicting the spreading rate of alien species, these variables are often not the limiting factor when making predictions on the scale of Flanders (Boets et al. 2014). Nevertheless, these variables should be taken into account for future research that aims to analyse and predict particular preferred conditions of alien molluscs at a larger scale.

Model performance

Predictive models are widely applied to assess the environments or areas that alien species have invaded or would invade [START_REF] Pitt | Predicting Argentine ant spread over the heterogeneous landscape using a spatially explicit stochastic model[END_REF]Boets et al. 2013;Chen et al. 2015). In many cases, the performance of these predictive models ranges from fair to moderate [START_REF] Gabriels | Analysis of macrobenthic communities in Flanders, Belgium, using a stepwise input variable selection procedure with artificial neural networks[END_REF]Boets et al. 2013). In our study, the overall performance of the models was moderate to good. When using data from the past period, the performance was higher than when using data from the most recent period (Figure 5A-B). This could be attributed to less complex biotic interactions and to the limited frequency of occurrences of alien species in the past. In this context, alien species may invade those sites where competition is low and where few native species occur. As alien species start to spread, competition with native species increases, and thus some native species which had a small environmental range may disappear (i.e. Bythinella, Margaritifera, Marstoniopsis, Myxas) or decrease their occurrence and abundance (e.g. Aplexa, Theodoxus). This is epitomized by alien and native molluscs in the family Physidae. In the past, the alien Physella occurred at only few sites with a few individuals, while the native Physa abundantly and widely occurred. A contrasting relationship between the two taxa was observed for the last period, in part due to a declining trend in the abundance and spatial occurrence of the native Physa. This could be a result of competition for food and niche with the alien Physella. Biotic interactions are important in predicting species distribution (Araújo and Luoto 2007;Meier et al. 2010), and when included will generally increase the predictive performance of a model. The exclusion of biotic interactions in our study may explain the overall lower performance of the models based on the most recent period. Moreover, during the late 2000s the range of environmental conditions where alien species occur increased (Figure 4) while many native species recovered their range and density (Table 2), due to improved water quality. This higher co-occurrence might also explain the lower performance of the models.

Similarly, higher co-occurrence may have influenced the higher reliability in the prediction of the "co-occurrence" of alien and native molluscs for the recent period compared to the first three periods. Likewise, models predicting the "native" sites yielded a high performance when the prediction was based on the periods (i.e. the first three periods) that have a large sample size of the "native" sites. This is quite logical as for predictive models the more input samples provided the better the models learn, and as a result, a higher predictive performance can be obtained (Stockwell and Peterson 2002;Hernandez et al. 2006). On the other hand, the models were not able to predict the "alien" sites (Appendix 2), due to a low number of instances of this particular class. Although the number of samples of the "alien" sites increased in the most recent period, it was still not sufficient for the models to learn and make a correct prediction. Small sample sizes, together with the opportunistic and generalistic characteristics of the alien species [START_REF] Nehring | Four arguments why so many alien species settle into estuaries, with special reference to the German river Elbe[END_REF], are therefore considered the main reasons for the models to yield a very low performance. Additional observations or a particular optimization approach is thus recommended to better predict the "alien" sites and evaluate the predictive power of the models.

Article 4

Conclusion

From the past to the most recent situation, there is an increasing trend in the spatial co-occurrence of alien and native molluscs in Flanders. Co-occurrence was predicted to mainly occur in rivers having low sinuosity and good chemical water quality. In addition, our most recent data indicated that alien molluscs have reached a relatively high number of sites where natives were not present, indicating either that alien molluscs have invaded more new sites or replaced native species at sites where they previously occurred. Given that our models were not able to make reliable predictions for environmental conditions preferred by alien molluscs, additional predictors and observations or perhaps a particular optimization approach is needed to predict the habitat conditions where alien molluscs are able to dominate the community. These results provide important information regarding the past and current coexistence of alien and native molluscs in Flanders. Our findings may be used to support management and conservation planning.

INTRODUCTION

Freshwater ecosystems are exposed to a high degree of isolation and endemism, and to a range of environmental and human pressures, which make these systems vulnerable to biological invasions (Dudgeon et al., 2006;[START_REF] Richter | Threats to imperilled freshwater fauna[END_REF]. The distribution of alien invasive species in freshwater ecosystems is increasing and is becoming one of the leading concerns for global economy and biodiversity (Sala et al., 2000). Major ecological effects of invasions cover replacement of native species, changes in nutrient cycling, reorganization of food webs and simplification of community composition (Bernauer and Jansen, 2006;Didham et al., 2007;Gurevitch and Padilla, 2004), all of which in turn alter the natural mechanisms of ecosystem processes.

Insights for practical assessment and management of invasive species have been revealed from previous studies which are based on predictive models (Boets et al., 2013;Chen et al., 2015;[START_REF] Elith | Species distribution models: ecological explanation and prediction across space and time[END_REF][START_REF] Hill | A predicted niche shift corresponds with increased thermal resistance in an invasive mite, Halotydeus destructor[END_REF]. For example, modelling approaches have been used to predict distributional areas of invasive/alien plants, herbs and macroinvertebrates (Boets et al., 2013;Chen et al., 2015;Roura-Pascual et al., 2009;Thuiller et al., 2005;Vicente et al., 2011). Approaches employed include generalized linear models, generalized additive models, classification and regression trees, random forest and artificial neural networks, all of which are data-driven or knowledge based, and are considered important tools for predicting species distributions (Guisan and Thuiller, 2005).

While many available predictive approaches are widely implemented, a wide variation of their performance has been frequently reported (Elith et al., 2006;Guisan et al., 2007;Segurado and Araujo, 2004). Some models even yield contrasting predictions of habitat suitability (e.g. Guisan et al. 2007;Evangelista et al. 2008;Roura-Pascual et al. 2009). As such, when calibrating and validating predictive models, data characteristics such as sample size, species prevalence or environmental predictors have to be carefully taken into account because they can influence the models' performance (Dormann et al., 2008;Luoto et al., 2006). Furthermore, predictive models are sensitive to parameterization and selection criteria during the modelling process (Araújo and Guisan, 2006;Elith et al., 2006), and thus can result in an uncertainty of current or past/future projections of species distributions (Buisson et al., 2010;Svenning et al., 2008).

In this context, a model that can provide better insights is of prime value for assessment and management planning. Provided that the eradications of successful invaders is hardly possible (Regan et al., 2006), assessment and control of invasive/alien species are thus an urgent need to maintain biodiversity. A long-term investigation on the invasion process/success requires a rigorous analysis, which is mostly based on statistical models and their interpretation [START_REF] Blossey | Before, during and after: the need for long-term monitoring in invasive plant species management[END_REF]. As most of the predictive models show varied performance when predicting species occurrence (Buisson et al., 2010;Svenning et al., 2008), improvement of the reliability of a predictive model via testing it through "hindcasting" -predicting what happened in the past -is subsequently recommended [START_REF] Pearman | Prediction of plant species distributions across six millennia[END_REF][START_REF] Sanders | Hindcasting helps scientists improve forecasts for life on Earth[END_REF]. Only when a model can provide a clear understanding of what is likely to take place in the ecosystem, managers and decision makers can construct policies to mitigate unwanted impacts. Therefore, an optimization of a particular predictive model, which can be robust and yields a relatively stable performance, will provide a better insight for further research conducted to confidently predict species distributions.

In this study, we optimized the reliability of a predictive algorithm (i.e. Classification Tree) in hindcasting and forecasting alien mollusc occurrence in river systems in Flanders, Belgium. We built 144 models based on a dataset of 20 years (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) comprising 7695 field data points, and 144 models based on a combined dataset, which is referred to the combination of the field data and cloned data (i.e. a dataset obtained by independently duplicating the filed data by k different individuals). In the calibration and validation process, we used random and stratified resampling approaches to split the data into different folds and set different levels of the models' complex parameters. We then evaluated how the two types of the dataset (field and combined data) and the two resampling approaches, under same and different parameterization settings, affect the model optimization. We expected that the models which were calibrated and validated using the stratified-split data are more robust and yield a greater reliability, compared to those being calibrated and validated using the random-split data, as has been shown in Hirzel & Guisan (2002).

MATERIAL AND METHODS Data collection and processing

The Flemish Environment Agency (VVM) has been collecting biological and environmental data in Flanders since 1989. The samples were collected at more than 2500 sites spread over different water bodies. For detailed information on how the data were collected, we refer to Sor et al. (2016). Data from 1991 to 2010 was used in the present study. This dataset was then divided into 4 periods; period D1: data from 1991 to 1995, D2: 1996 to 2000, D3: 2001 to 2005 and D4: 2006 to 2010. The response variable of each site was grouped as a categorical variable: a "native" site (i.e. a site having only the native taxa present), an "alien" site (i.e. a site having only the alien taxa present) and an "overlap" site (i.e. a site having both the native and alien mollusc taxa present). Since the number of occurrence instances of most alien mollusc genera is limited (Table 1), we decided to group them into one categorical variable. In this way, we had a higher number of instances for our predictive models and thus a better and more robust model that could be developed. We removed physical and chemical water quality variables containing missing values in more than 5% of all samples. Finally, each period consisted of one categorical response variable (native/alien/overlap) and 11 environmental predictors, which are summarized in Table 2.

Modelling algorithm

We used a Classification Tree (CT) model to predict the occurrence of the response variable. The CT is known as an easy-to-read and easy-to-interpret method (Guisan and Thuiller, 2005). Due to this simplicity and interpretability, and its ability to tolerate missing values during the calibration and validation process (De'ath and Fabricius, 2000;Therneau and Atkinson, 1997), CT is mainly implemented to model complex biological interactions and is suggested for management applications (Guisan and Thuiller, 2005;[START_REF] Waite | Comparison of watershed disturbance predictive models for stream benthic macroinvertebrates for three distinct ecoregions in western US[END_REF]. For a detailed description of the theoretical growing and pruning process of the CT models, we refer to Quinlan (1986) and Therneau & Atkinson (1997).

Hindcasting, forecasting and validation

Field data

For the hindcasting, we only used data period D4 to calibrate the models. Then the calibrated models were used to hindcast (validate) the response variable based on the environmental predictors of the data from period D4, D3, D2 and D1. During the pruning phase, the complex parameter of the calibrated models was set at cp=0.05. Where the tree had only a root, we decreased the cp to a level (e.g. cp between 0.04 -0.01) that at least two terminal nodes were produced. The first few splits usually provide a very informative division of the data (Therneau and Atkinson, 1997). These criteria were set in order to make the trees easily interpretable and comparable in terms of the number of variables and complexity.

For the forecasting, we only used data period D1 to calibrate the models. Then we used the calibrated models to forecast (validate) the response variable based on the environmental predictors of the data from period D1, D2, D3 and D4. The setting of the complex parameter of the calibrated models was similar to the settings as described in the "hindcasting" section.

The calibration and validation process was based on a three-fold cross-validation (3CV). The 3 fold CV was obtained by splitting the whole dataset into three folds. Two folds of the 3 fold CV were in turn used for the calibration, whereas one fold was used for the validation. In other words, we used two-thirds to calibrate and one-third to validate the models, and repeated the procedure three times so that at the end each fold was used exactly once for the validation (hindcasting and forecasting). For the data splitting, we used a random and a stratified resampling approach.

The random resampling was done by shuffling and randomly splitting the whole dataset into three folds. From this approach, the sample size of each class was not equally stratified across the three folds. To allow a reliable error estimation of the hindcasting and forecasting models for each period, we made three replicates of the 3 fold CV procedure by reshuffling and randomly splitting the whole dataset again into 3 new folds.

For the stratified resampling, we first allocated the three classes of the response variable and the environmental predictors of the whole dataset into three separated data subsets (DS, i.e. DS1, DS2, DS3). Each DS corresponds to the data of each class of the response variable. Then, we randomly divided each DS into 3 smaller data subsets (i.e. DS1a, DS1b, DS1c; DS2a, DS2b, DS2c; DS3a, DS3b, DS3c). Thus, each of the smaller data subsets of the same class had exactly the same number of samples (n, e.g. nDS1a = nDS1b = nDS1c). For any DS that could not be equally divided by 3, each of the smaller data subsets may have 1 sample more or less than the others (e.g. nDS1a ± 1). Finally, we recombined each smaller data subset to obtain new three equal-stratified data subsets, i.e. Fold1 (nDS1a + nDS2a + nDS3a), Fold2 (nDS1b + nDS2b + nDS3b) and Fold3 (nDS1c + nDS2c + nDS3c), which were later used to calibrate and validate the models. The schematic diagram showing this resampling procedure is depicted in Fig. 1. To be consistent, we also made three replicates of the 3 fold CV for this resampling approach, by reshuffling the original data and then following the above-described procedure. The calibration and validation of the models was performed separately for each resampling approach.

Combined data

In most cases, the "native" class had the largest sample size, followed by the "overlap" and the "alien" class. Due to the unequal sample size, we developed a cloned dataset by independently duplicating the original filed data of the less frequently occurring classes by k different individuals. Then we combined the cloned data with the original field data to obtain the same sample size for each class (Fig. 2). One may argue to use an equal-stratified dataset of field observations by randomly selecting a number of observation of the more occurring classes equally to the less occurring class. If following this, however, we will lose incredible valuable-information of the field observations and that leaves a small sample size for the model development, which consequently leads to a less reliable and less robust models. On the contrary, using the cloned data allows us to incorporate all field observations and provides us an equal data distribution and sufficient samples for each period. Incorporating the cloned data into the models will increase the maximum likelihood of the prior class distribution, making the models more robust (Lele et al., 2007). To run the hindcasting, forecasting and validations for each period using the combined data, we followed the criteria and procedures used in the "field data" section.

Model performance and analysis

For each scenario (hindcasting and forecasting), nine models (3 models of 3 fold CV × 3 replicates) were built for each resampling approach in each period (see Supporting Information Table S1). In this way, a performance value of the nine different models was calculated. From the field data, we built 36 models ([3 models of 3 fold CV hindcasting × 3 replicates × 2 resampling types] + [3 models of 3 fold CV forecasting × 3 replicates × 2 sampling types]) for each period. In the same way, 36 models were built based on the combined data for each period. Therefore, we totally built 144 models over the four periods for each data type (field and combined data).

Cohen's Kappa Statistic (K) and Correctly Classified Instances (CCI) were used to evaluate the model performance. The higher the value of K (ranging from 0 to 1) and of CCI (ranging from 0 to 100), the better the model predicts the response variable. To obtain the best model configuration in both scenarios, we analysed the model performance based on 1) the data types, 2) the resampling approaches and 3) the data types and resampling approaches. All the analyses were carried out using the language program R ( R Core Team, 2013).

RESULTS

Field data

For the hindcasting scenario, the yielded predictive performance (i.e. Kappa and CCI) of the CT models were Kappa = 0.16, 0.16, 0.11 and 0.07, and CCI = 55%, 58%, 60% and 50% for the data period D4, D3, D2 and D1 (Fig. 3). The yielded Kappa and CCI for the forecasting scenario were 0.35, 0.07, 0.07 and 0.03, and 93%, 91%, 79% and 47%, respectively. The performance of the models, which were calibrated and validated using the two resampling approaches, showed a relatively similar trend and stability (see Supporting information Fig. S1). However, the CT models were not able to hindcast and forecast the class "alien", although they could hindcast and forecast the other two classes (see Table S2).

Combined data

The Kappa (0.24, 0.16, 0.13, 0.02) and overall CCI (48%, 44%, 42%, 24%) of the models decreased from period D4, D3, D2 to D1 for the hindcasting scenario. For the forecasting scenario, the model performance decreased from period D1 to D4; the Kappa was 0.58, 0.17, -0.01 and -0.05, and overall CCI was 72%, 44%, 33% and 30% for the corresponding periods (Fig. 4a-b). The calibrated and validated models using the two resampling approaches also yielded a similar trend in predictive performance and in stability for each scenario (see Fig. S1). Interestingly, the models were able to hindcast and forecast the class "alien" and the other two classes. The CCI of the class "alien" followed the decreasing trend of the overall performance for each scenario; it was 76%, 62%, 44% and 9%, and 100%, 44%, 20% and 14% for the corresponding periods used in the hind-and forecasting scenario, respectively (Fig. 4c).

DISCUSSION

In this study, we optimized the hind-and forecasting of alien mollusc occurrence to obtain the best CT model configuration for both in terms of performance and interpretation. Variation in class distribution and resampling design were investigated and strongly indicated that a similar sample size of each class should be considered as has been generally suggested (e.g. Manel et al. 2001;Allouche et al. 2006). Based on field data (field observations), for which the prevalence of each predicted class was not equally distributed, the CT models were unable to predict the putative occurrence of the alien molluscs. These results indicate that the field observations could not be predicted. However, by using the combined data, their occurrence could be correctly predicted at periods other than when the CT models were calibrated, which corresponds to the field past and current situation.

Improving the model reliability by testing it to hindcast the occurrence of a given species in the past could be decisive for predicting the future occurrence of that species. Although predictive models have been substantially implemented in ecology in the last two decades [START_REF] Guisan | Predicting species distributions for conservation decisions[END_REF]Guisan and Thuiller, 2005), only a few have tested each model's performance via hindcasting (e.g. Svenning et al. 2008;[START_REF] Espíndola | Predicting present and future intra-specific genetic structure through niche hindcasting across 24 millennia[END_REF][START_REF] Maire | Hindcasting modelling for restoration and conservation planning: Application to stream fish assemblages[END_REF][START_REF] Pelletier | Historical species distribution models predict species limits in western plethodon salamanders[END_REF]. Most of the existing studies have used currently known species occurrences to hindcast or forecast the past or future putative species occurrence (Boets et al., 2013;Buisson et al., 2010;[START_REF] Maire | Hindcasting modelling for restoration and conservation planning: Application to stream fish assemblages[END_REF]Svenning et al., 2008), which sometimes results in a high uncertainty of the models used. In the present study, we used both the hindcastingand forecasting-based models to predict the occurrence alien molluscs. As a result CT models could make a correct prediction of the past and recent occurrence.

One of the challenges in predicting a species' occurrence is obtaining a balanced class distribution for the response variable, because a small sample size of each predicted class can be the source of instability and errors in species distribution models (Allouche et al., 2006;McPherson et al., 2004). However, our results demonstrated that having the exact same class distribution (i.e. stratified split) is not always necessary during the calibration and validation process because the models can make the correct prediction and yield a similar performance when the prevalence reaches a certain threshold (see Fig. S1). However, the model is incapable of predicting a class with a very low prevalence, as was the case of the class "alien" in this study (see Table S2). As collecting new field data from past/recent periods is unfeasible/costly, optimization appears necessary (Hirzel and Guisan, 2002). In respond to this, we optimized the prediction by considering cloned data for the models, resulting in a better prediction. Cloned data have been applied and recently suggested for hierarchical models in ecology (Lele et al., 2010(Lele et al., , 2007;;Ponciano et al., 2009). The increasing number of clones used in the models can increase the maximum likelihood of the prior class distribution, but it does not affect the statistical accuracy, which mainly depends on the information of the field data and the model calibration and validation process (Lele et al., 2007).

Parameterization and selection criteria during the modelling process are known to be sensitive to model performance (Araújo and Guisan, 2006;Elith and Graham, 2009). In this study, the complex parameter (cp) ranging from 0.05 to 0.01 was used to optimize the tree selection of the calibrated models. In most cases, this setting resulted in a tree having 2 to 4 nodes, with an increased coefficient of determination (R 2 ) of the model cross-validation (i.e. reducing crossvalidation error) (see Fig. S2). The trees having a few nodes and yielding an increased R 2 are known to provide a good reliability and representation of data division [START_REF] Therneau | An introduction to recursive partitioning using the rpart routines[END_REF]Therneau and Atkinson, 1997). This suggests that the complexity parameter and selection criteria set in our study may not lead to over-fitting of the hind-and forecasted models.

In the best model configuration (using the combined data), environmental variables, which are considered as important factors promoting the occurrence and abundance of alien molluscs, were incorporated in the calibrated and validated models. They include water quality variables such as nitrate, ammonium, conductivity, pH, and chemical oxygen demand (Boets et al., 2013;Grabowski et al., 2009;Strayer, 2010;[START_REF] Vermonden | Environmental factors determining invasibility of urban waters for exotic macroinvertebrates[END_REF]. The water quality conditions in Flanders, which have been substantially improved over the past decades (MIRA, 2012), could be the factors leading to a gradual decrease in global performance of both the hind-and forecasting models in our study. This result might also reveal the natural characteristics of predictive models that the uncertainty increases when longer predictions of the future are made (Buisson et al., 2010;[START_REF] Pearman | Prediction of plant species distributions across six millennia[END_REF]. Furthermore, alien species are mostly considered generalists, being able to withstand a wide range of environmental conditions (Werner & Rothhaupt 2007;Pigneur et al. 2014). This characteristic, along with their currently increased migration patterns and the changes in water quality conditions, could influence the correct prediction of the class "alien" when predicting the past or to the present (Fig. 1). Nevertheless, although the CCI decreased from recent to past periods (hindcasting scenario) and vice-versa (forecasting scenario), the CT models could still correctly predict the occurrence of alien molluscs with a low error rate (Fig. 4), demonstrating the robustness of the CT models.

While we have incorporated all important environmental factors to optimize the CT models, available data of dispersal vectors might provide a better model output. Although these vectors are not always a limiting factors, previous studies have shown its contribution in predicting the distribution of alien species (Boets et al., 2014(Boets et al., , 2013)). Regardless the dispersal vectors, however, if the alien species maintain their climatic niche, the CT models will be able to predict their future occurrence on the basis of the current environment conditions in which they have already established (Parravicini et al., 2015). Otherwise, climatic variables should be considered in the models, because such variables, e.g. temperature, have been illustrated to associate with the establishment success of alien species (Leuven et al., 2009;[START_REF] Werner | xerror=0.93↓ Replication Replication Replication 5 Field data_Hindcasting_Stratified Models calibrated on D4 and validated by D4,D3,D2,D1 cp=0[END_REF]. Consideration of incorporating these variables and applying the hind--and forecasting-based approach outside Flanders will also help to test the robustness of CT models and to improve our understanding of the behavior and ecological characteristics of alien molluscs, which are imperative for assessment and management purposes.

CONCLUSION

Since the rate at which studies employing predictive models increases and that their performance uncertainty is always noticed, evaluating the reliability of a particular model via hindcasting is a proposition for further studies aiming to predict future species occurrence.

Here we report the robustness of CT models in correctly predicting alien mollusc occurrence by incorporating cloned data into the models. With the same parameterization, the performance of the CT models is unlikely to be affected by resampling approaches (i.e. random vs stratified) used to split the data for calibrating and validating the models.

Our optimization success provides insight for the application of predictive models. Although model reliability decreased when predicting further into the past or into the future, the models still provide a correct prediction with a low error rate. The CT model is therefore suitable for assessment and management applications. Incorporating dispersal vectors and climatic variables into the models and testing model transferability into a new area may also help to increase the robustness of the CT models.

List of tables and figures with captions

Table 1. The occurrence instances of each alien mollusc genus and of all genera that are merged together, compared to all collected samples (7695 samples) within the studied period (1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). Table S2. Confusion matrices showing the instances of the observed and predicted classes obtained from the CT models of the hindcasting and forecasting scenarios for each period. The instances of each class are averaged from the 9 models which were calibrated and validated using random resampling approach. The three classes are overlap (O), Alien (A) and Native (N). Fig. 2 The CCI of the predicted classes of models using the field data.

Fig. 3 The overall performance of models using the combined data.

Fig. 4 The CCI of the predicted classes of models using the combined data.

  កូ នសូ មអរគ ុណេ កពុ ក អន កម ្ដ យជអេនក ស្រមប់ ករច ិ ញច ឹ មប ី បច ់ែថរក ន ិ ងទំ ន ុកប្រមុ ង ្រគប់ ែបបយ៉ ង។ អរគ ុណបង្របុ ស ន ិ ងបង្រសី ែដលេម ើ លែថ ន ិ ងករពរប្អ ូ ន ំងព ី តូ ចជមួ យ ពុ កែម៉ ។ េសចក្ដ ី ្រស ញ់ ករផគ ត់ផគ ង ់ ន ិ ងករេល ើ កទ ឹ កច ិ ត្តរបស់ អន កទ ំងអស់ គន គម នអ្វ ី កត់ៃថ្ល បន េឡ ើ យ។ ក្ដ ី សុ ខរបស់ អន កទ ំងអស់ គន ជក្ដ ី សុ ខរបស់ ខ ញ ុ ំ ែដរ។

Fig. 1 . 1

 11 Fig. 1.1 Freshwater food web in the Tonle Sap Lake and River ecosystem, Cambodia. The size of each circle represents the biomass (tonne/km 2 ) of each functional group (e.g. shrimps, crabs, molluscs etc. (Chea et al. 2016).

Fig. 2 . 1

 21 Fig. 2.1 The Lower Mekong Basin (LMB, A) and macroinvertebrate sampling sites (shaded dots, B). Subsamples and replicates were taken at each sampling site as illustrated in C.

Fig. 2 . 2

 22 Fig. 2.2 Map of Flanders indicating: (A) the most important watercourses and geographic locations, the polder area (grey) and the three main harbours indicated by rectangles (Boets et al. 2016), (B) different river basins (van Griensven and Vandenberghe 2006) and (C) monitoring sites between 1991-2010, which were used in the present study.

  build reliable models and to avoid misidentifying the key variables determining each class of the response variable, 3 replications of the three-fold cross validation was made. For each 3fold cross-validation, the data was shuffled and randomly split into three subsets; two subsets were used for training and one subset for validation. For the second and third replication, I reshuffled and randomly split it into new training and validation sets following the same procedures. From each training and validation set, a model was built and in this way, a performance value and the importance of each variable (in percentage) of nine different models (3 models of each three-fold cross validation × 3 replications) were calculated. A mean performance value, obtained from the nine models, was used as a final criterion for model evaluation.Cohen's Kappa Statistic (Kappa) and Correctly Classified Instances (CCI) were used to evaluate the model performance. The importance of each variable determining the preferred environmental conditions of each class (native/alien/co-occurrence) was averaged

Fig. 2 .

 2 Fig. 2.3 A) Bar charts showing the number of the field and cloned data, which were combined together and used in the models. B) Schematic diagram illustrating the splitting procedure of stratified resampling approach.

  3.1). Cluster I was situated in the Mekong delta in Vietnam; cluster IIa, along the upstream sites in Laos and Thailand; cluster IIb1, mostly in middle part in Cambodia and a few sites in Laos, Thailand and Vietnam; and cluster IIb2, mostly located in the tributaries of the LMB (Fig.3.1A).

Fig. 3 . 1

 31 Fig. 3.1 Sampling sites and the four clusters, representing four macroinvertebrate assemblages, based on the cluster analysis (A), and the dendrogram showing sites belonging to the four clusters based on the dissimilarity and Ward's hierarchical clustering method (B).

Fig. 3 . 2

 32 Fig. 3.2 Box and whisker plots of richness (A) and abundance (B) of macroinvertebrate assemblage in each cluster and its proportion of mean richness (C) and abundance (D) consisting of different components of macroinvertebrates.

  results of the LDA model used to discriminate the macroinvertebrate assemblages based on the physical-chemical variables and land cover types are shown in Fig. 3.3. Along axis 1, assemblage I was situated opposite to assemblage IIa and IIb2. Assemblage I was positively correlated with the surface area of watershed, river depth, river width and water temperature, but negatively associated with altitude and dissolved oxygen. Whereas assemblage IIa was positively correlated with electrical conductivity and urban areas, and assemblage IIb2 was positively linked to altitude, DO, Secchi depth, wood-/shrub-land and evergreen forests. Based on axis 1 and 3, assemblage IIb1 was positively linked to inundated, wetland and agricultural areas (Fig. 3.3).

Fig. 3 . 3

 33 Fig. 3.3 Results from the LDA discriminating the four clusters (I, II, IIb1, IIb2), representing four macroinvertebrate assemblages, using Axes 1, 2 and 3 that explained the indicated percentage of the total variance in the data (A, C), and correlations of the environmental factors to the corresponding axes (B, D). ALT: altitude, RW: river width, RD: river depth, SD: Secchi depth, WT: water temperature, DO: dissolved oxygen, EC: electrical conductivity, SAW: the surface area of watersheds, Agr: agricultural land, Bmb: bamboos, Crp: crops, Dec: deciduous forests, Evg: evergreen forests, Grs: grassland, Ind: inundated, Mix_evg.dec: mixed evergreen and deciduous forests, Plt: plantations, Reg: regrowth, Roc: rocks, Urb: urban areas, Wat: water, Wet: wetland, Wod: wood-& shrub-land.

Fig. 3 . 4

 34 Fig. 3.4 Box and whisker plots of macroinvertebrate richness (MR, A) and log+1 abundance (MA, B) between seasons in each period. Significant difference in richness and abundance for each season across the four periods is indicated as "a", "b", "c". * p<0.05, *** p<0.001.

  samples are more related to high values of sinuosity, and in the early (D3) and late (D4) 2000s, samples were more related to high values of DO and pH, respectively (Fig. 3.5).

Fig. 3 .

 3 Fig. 3.5 RDA ordination plots showing the association of macroinvertebrates with environmental variables. A) Plot showing sampling years and river basins; Name of river basins: Ijzer, Leie, Maas, Nete, Demer, Dender, Bru-Pol: Brugse Polders, Gen-Kan: Gentse Kanalen, Ben-Sch: Beneden-Schelde, Bov-Sch: Boven-Schelde, and Dij-Zen: Dijle Zenne. B) Plot showing macroinvertebrate taxa and environmental factors for the same analysis. Code of taxa with small loadings score were removed to improve legibility; see Table3.5 for full names of the taxa.

Fig. 3 . 6

 36 Fig. 3.6 Bar and standard error plots showing α diversity of the global and component communities in each period (A), and box and whisker plot showing α diversity (SH) of global communities between seasons (B). Significant difference in α diversity for each season across the four periods is shown as "a", "b", "c" (B). *** p<0.001.

Fig. 3 . 7

 37 Fig. 3.7 Bar and standard error plots showing β diversity of global communities between periods (A) and between seasons of each period (B).

Fig. 3 . 8

 38 Fig. 3.8 Performance of predictive modelling techniques. A) Performance based on different prevalence ranges, B) Behaviour of the performance based on the complete prevalence range.

Fig. 3 .

 3 Fig. 3.10 Mean and standard error bars indicating the overall model performance based on Kappa (A) and CCI (B), and the predictive power of the predicted class "co-occurrence" (C) and the class "native" (D). CV: cross-validation, Reps: Replications.

Fig. 3 .

 3 Fig. 3.11 Mean and standard error bars indicating the overall performance of CT models based on Kappa and CCI, when using the combined data and random resampling approach.

Fig. 3 .

 3 Fig. 3.12 Mean and standard error bars indicating the performance of CT models based on the correct prediction (CCI) of the predicted class "alien".

Fig. 4 . 1

 41 Fig. 4.1 RDA plots showing the association of (A) the centroid of macroinvertebrate communities (points) with environmental variables in each period and of (B) the relationship between different taxa, environmental variables and river basins: Ijzer, Leie, Maas, Nete, Demer, Dender, Bru-Pol:Brugse Polders, Gen-Kan: Gentse Kanalen, Ben-Sch: Beneden-Schelde, Bov-Sch: Boven-Schelde, and Dij-Zen: Dijle Zenne. Taxa (red colour) that had small loading scores on the RDA axes were removed to improve visibility; see Table3.5 for full name of the taxa.

  Contents lists available at ScienceDirect Limnologica journal homepage: www.elsevier.com/locate/limno Spatial organization of macroinvertebrate assemblages in the Lower Mekong Basin Ratha Sor a,b,c, ⁎ , Pieter Boets b,d , Ratha Chea a , Peter L.M. Goethals b , Sovan Lek a a Laboratoire Evolution & Diversité Biologique, UMR 5174, Université Paul Sabatier -Toulouse III, 118 route de Narbonne, 31062 Toulouse cédex 4, France b Laboratory of Environmental Toxicology and Aquatic Ecology, Ghent University, Campus Coupure building F, Coupure links 653, B9000 Ghent, Belgium c Department of Biology, Faculty of Science, Royal University of Phnom Penh,Russian Boulevard, 12000, Phnom Penh, Cambodia d Provincial Centre of Environmental Research, Godshuizenlaan 95, 9000 Ghent, Belgium A R T I C L E I N F O

Fig. 1 .

 1 Fig. 1. Sampling sites and the four clusters, representing four macroinvertebrate assemblages, based on the cluster analysis.

Fig. 2 .

 2 Fig. 2. Dendrogram showing sites belonging to the four clusters, representing four macroinvertebrate assemblages, based on the Bray-Curtis dissimilarity and Ward's hierarchical clustering method.

Fig. 3 .

 3 Fig. 3. Box and whisker plots of richness (A) and abundance (B) of macroinvertebrate assemblage in each cluster and its proportion of mean richness (C) and abundance (D) consisting of different taxonomic groups of macroinvertebrates.

Fig. 4 .

 4 Fig.4. Results from the LDA discriminating the four clusters (I, II, IIb1, IIb2), representing four macroinvertebrate assemblages, using Axes 1, 2 and 3 that explained the indicated percentage of the total variance in the data (A, C), and correlations of the environmental factors to the corresponding axes (B, D). ALT: altitude, RW: river width, RD: river depth, SD: Secchi depth, WT: water temperature, DO: dissolved oxygen, WC: water conductivity, SAW: the surface area of watersheds, Agr: agricultural land, Bmb: bamboos, Crp: crops, Dec: deciduous forests, Evg: evergreen forests, Grs: grassland, Ind: inundated, Mix_evg. dec: mixed evergreen and deciduous forests, Plt: plantations, Reg: regrowth, Roc: rocks, Urb: urban areas, Wat: water, Wet: wetland, Wod: wood-& shrub-land.

Fig. 1 .

 1 Fig. 1. Map of the sampling sites in the LMB. (a) LCBD indices with significant p-values uncorrected (brown dots) and corrected for multiple testing by applying Holm correction (shaded dots with star); open circles: nonsignificant LCBD indices. (b) Richness (number of taxa) of the sampling sites. Three red dots indicated with red arrows: lowest richness, 6; large shaded circle: highest richness, 74. The richness for the five sites having significant LCBD indices is shown in the parentheses. The sizes of the circles are proportional to LCBD (a) or richness (b) values.

Fig. 2 .

 2 Fig. 2. RDA ordination plots showing the association of taxa having SCBD higher than the mean value (i.e. important taxa) with the geo-environmental factors. (a) Plot showing the sites and geo-environmental factors. VB4, CM3, TM1, TM2 and TM3 in bold are the sites with significant LCBD indices. LAT: latitude, ALT: altitude, WT: water temperature, DO: dissolved oxygen, WC: water conductivity, RW: river width, RD: river depth, SD: Secchi depth. (b) Plot showing taxa with high SCBD indices and geo-environmental factors for the same analyzes. The first letter of each taxon name represents the component of macroinvertebrate communities (A: annelid, C: crustacean, B and G: mollusk, and I: insect). Geo-environmental variables and taxa that have small loading score were removed to improve legibility.

Fig. 3 .

 3 Fig. 3. Variation partitioning results showing the fractions of the global LCBD variation explained by insect and mollusk communities (LCBDI + LCBDM), by the three environmental variables (WC: water conductivity, RD: river depth, SD: Secchi depth), and by the two sets of explanatory factors.

Fig. 1 .

 1 Fig. 1. The Lower Mekong Basin (LMB) and macroinvertebrate sampling locations (shaded dots).

Fig. 2 .

 2 Fig. 2. Box and whisker plots of normalized environmental variables. Rectangles delineate the first and third quartiles, dark bars are the medians, the lower and upper bars are the minima and maxima, and the circles are outliers. ALT: altitude, DFS: distance from the sea, RW: river width, RD: river depth, SD: Secchi depth, WT: water temperature, DO: dissolved oxygen, WC: water conductivity.

Fig. 3 .

 3 Fig. 3. Effects of species prevalence on the performance measures of the predictive models. Each symbol is the mean performance value of the corresponding modelling techniques. The lines are LOWESS smoothers. Kappa: Cohen's Kappa Statistic, AUC: area under the curve, LR: logistic regression, ANN: artificial neural network, RF: random forest, SVM: support vector machine.

Fig. 4 .

 4 Fig. 4. Mean values and standard error bars showing the performance measures of each modelling technique based on different species prevalence ranges. Kappa: Cohen's Kappa Statistic, AUC: area under the curve, LR: logistic regression, ANN: artificial neural network, RF: random forest, SVM: support vector machine, n: the number of predicted species corresponding to each prevalence range.

Figure 1 .

 1 Figure 1. Sampling locations indicating the occurrence of native and alien molluscs and their co-occurrence for each period.

Figure 2 . 3 -:

 23 Figure 2. Classification trees predicting the "co-occurrence" and the "native" sites for each period. SI: Sinuosity, NH 4 + : Ammonium, NO 3 -: Nitrate, COD: Chemical Oxygen Demand, cp: pruning complex parameter.

Figure 3 .

 3 Figure 3. Mean and standard error bar showing the importance of the variables contributing most to the predictive models for each period. COD: Chemical Oxygen Demand, EC: Electrical Conductivity, NH 4 + : Ammonium, NO 3 -: Nitrate, SI: Sinuosity. Variables that contributed less to the predictive models are not shown.

Figure 4 .-Figure 5 .

 45 Figure 4. Box and whisker plots of physical-chemical variables in which each occurrence type occurred. Rectangles show first and third quartiles, dark bars are the medians, the lower and upper bars are the minimum and maximum values, and the circles are outliers. EC: Electrical Conductivity, DO: Dissolved Oxygen, COD: Chemical Oxygen Demand, NH 4 + : Ammonium, NO 3 -: Nitrate, NO 2 -: Nitrite, Pt: Total Phosphorus, oPO4: Orthophosphate.

Fig. 1

 1 Fig. 1 Schematic diagram illustrating the splitting procedure of stratified resampling approach.

Fig. 2

 2 Fig. 2 Bar charts showing the number of the field and cloned data, which were combined together and used in the models.

Fig. 3

 3 Fig. 3 Mean and standard error bar indicating the overall performance of CT models based on Kappa (a) and CCI (b), when using the field data and random resampling approach.

Fig. 4

 4 Fig. 4 Mean and standard error bar indicating the performance of CT models based on Kappa and CCI, when using combined data and random resampling approach. (a) and (b) are the overall performance, and (c) illustrates the correct prediction of the predicted class "alien".

Table 2 . 1

 21 Mean value (and standard deviation) for environmental predictors, and occurrence instances of each class of the response variable. NH 4

	+ : Ammonium, COD: Chemical Oxygen Demand, Pt: Total Phosphorus, EC: Electrical Conductivity, NO 3 -: Nitrate, NO 2 -: Nitrite, oPO 4 : Orthophosphate, DO:
	Dissolved Oxygen. The number of sampled sites for each period is shown in square brackets.
				Period
			1991-1995	1996-2000	2001-2005	2006-2010
	Variable	Unit	[509]	[991]	[1524]	[1250]
	NH 4 +	mg/L	2.6 (4.8)	1.9 (3.4)	1.6 (2.8)	2.3 (5.1)
	COD	mg/L	55 (43)	39 (42)	34 (34)	36 (41)
	Pt	mg/L	1.0 (1.4)	0.9 (1.4)	0.9 (2.3)	0.8 (1.0)
	EC	µS/cm	1320 (2618)	987 (1149)	998 (1509)	921 (949)
	NO 3 -NO 2 -	mg/L mg/L	3.3 (4.1) 0.2 (0.3)	4.0 (4.7) 0.2 (0.3)	3.6 (3.7) 0.2 (0.2)	3.0 (3) 0.2 (0.2)
	oPO 4	mg/L	0.6 (0.9)	0.5 (1.1)	0.4 (0.7)	0.5 (0.8)
	pH		7.5 (0.6)	7.6 (0.5)	7.7 (0.4)	7.6 (0.4)
	DO	mg/L	7.5 (3.5)	6.8 (3)	6.9 (3.1)	6.6 (2.9)
	Sinuosity		1.1 (0.1)	1.1 (0.1)	1.1 (0.1)	1.1 (0.1)
	Slope	m/1000m	1.0 (1.4)	1.6 (2.8)	2.0 (3.6)	1.7 (2.8)
	Response class
	Co-occurrence	66	146	567	875
	Alien		6	15	63	133
	Native		863	1860	2259	842
	Total instances	935	2021	2889	1850

Table 3 . 1

 31 Mean value (and standard deviation) for environmental variables, richness, abundance and Shannon diversity of macroinvertebrate assemblage in each cluster.

	Cluster [n]

Table 3 . 2

 32 List of indicator taxa (and their indicator values, IndVal) of macroinvertebrate assemblage in each cluster.

	Cluster I			Cluster I (continued)
	Annelid	IndVal p-value	Insect	IndVal p-value
	Aeolosoma bengalense	0.52	0.010	Arigomphus sp.	0.67	0.005
	Aulodrilus prothecatus	0.67	0.005	Cricotopus sp.	1	0.005
	Chaetogaster langi	0.85	0.005	Clinotanypus sp.	0.52	0.005
	Chaetogaster limnaei limnaei	0.6	0.005	Nectopsyche sp.	0.67	0.005
	Dero pectinata	0.67	0.005	Sigara sp.	0.6	0.005
	Dero sp.	0.74	0.005
	Dero sp.1	1	0.005	Cluster IIa
	Dero sp.2	0.95	0.005	Annelid	IndVal p-value
	Lumbriculidae sp.	0.6	0.005	Oligochaeta sp.	0.99	0.005
	Namalycastis longicirris	0.9	0.005	Polychaeta sp.1	0.6	0.005
	Orbinia johnsoni	0.52	0.010
	Polydora sp.	0.67	0.005	Mollusc
				Corbicula sp.	0.88	0.005
	Crustacean			Hubendickia sp.	0.6	0.010
	Corophium minutum	0. 8	0.005	Kareliania sp.	0.52	0.010
	Corophium sp.	0.67	0.005	Scaphula sp.	0.52	0.010
	Cyathura carinata	0.74	0.005	Stenothyra sp.	0.6	0.005
	Cyathura truncata	0.57	0.005
	Decapoda sp.	0.91	0.005	Insect
	Eohaustorius sp.	0.6	0.005	Anagenesia sp.	0.67	0.005
	Eohaustorius tandeensis	0.67	0.005	Caenoculis sp.	0.52	0.010
	Gammarus sp.	0.6	0.005	Caenodes sp.	0.74	0.005
	Grandidierella lignorum	0.78	0.005	Choropterpes sp.	0.51	0.005
	Grandidierella vietnamica	1	0.005	Dipseudopsis sp.	0.69	0.005
	Hyale hawaiensis	0.67	0.005	Heterocloeon sp.	0.52	0.010
	Hyale sp.	0.85	0.005	Micronecta sp.	0.6	0.005
	Kamaka sp.	0.6	0.005
	Macrobrachium equidens	0.6	0.005	Cluster IIb1
	Melita sp.	0.82	0.005	Mollusc	IndVal p-value
	Monocorophium sp.	0.91	0.005	Filopaludina filopaludina filosa	0.45	0.025
	Palaemon curvirostris	0.6	0.005
				Insect
	Mollusc			Pentagenia sp.	0.62	0.010
	Afropisidium clarkeanum	0.73	0.005
	Angulyagra polyzonata	0.6	0.005	Cluster IIb2
	Angulyagra sp.	0.9	0.005	Annelid	IndVal p-value
	Bithynia siamensis	0.67	0.005	Naididae sp.	0.76	0.005
	Corbicula baudoni	0.87	0.005
	Corbicula bocourti	0.74	0.005	Insect
	Corbicula leviuscula	0.97	0.005	Gomphidae sp.	0.56	0.010
	Corbicula moreletiana	0.86	0.005
	Corbicula sp.	0.95	0.005
	Gastropoda sp.	0.74	0.005
	Hyriopsis bialatus	0.64	0.005
	Limnoperna siamensis	0.99	0.005
	Limnoperna sp.	0.95	0.005
	Lymnaea viridis	0.94	0.005
	Mekongia swainsoni swainsoni	0.67	0.005
	Sinomytilus harmandi	0.9	0.005
	Stenothyra annandalei	0.6	0.005
	Stenothyra glabrata	0.85	0.005
	Trochotaia trochoides	0.52	0.005

Table 3 . 3

 3 List of taxa with high SCBD indices (above the overall mean). The first letter of each taxon

	.3). The SCBD values are

3.2 Macroinvertebrate communities and diversity patterns in Flemish rivers 3.2.1 Overall macroinvertebrate communities

	Ephemeroptera (17). Molluscs were mainly represented by the clade Hygrophila (17 taxa),
	Caenogastropoda (
	During the 20-year monitoring (1991-2010) across fluvial systems in Flanders, 207 taxa
	(6,192,056 inds) belonging to 145 families, 25 orders and 4 clades were identified. Only 123
	taxa were identified to genus level; the reaming 84 taxa were identified to only at the family
	level. Among all recorded taxa, 131 taxa (63%) and 1,752,885 inds (28%) belonged to insects,
	34 taxa (16%) and 628,468 inds (10%) to molluscs, 21 taxa (10%) and 1,589,809 inds (26%)
	to crustaceans, and 21 taxa (10%) and 2,220,894 inds (36%) to annelids. The dominant insect
	orders were Diptera (25 taxa), Odonata (25), Hemiptera (22), Trichoptera (18) and

Table 3

 3

	.4 Mean and standard error (SE) of richness and abundance of global and each component
	community for each period.
	Communities	D1	Richness (SE) D2 D3	D4	D1	Abundance (SE) D2 D3	D4
	Global	13.6 (0.12) 14.4 (0.15) 16.4 (0.15) 15.7 (0.19)	380 (32) 279 (20) 456 (24) 1009 (70)
	Annelids	3.1 (0.04) 3.3 (0.03) 3.6 (0.03) 3.4 (0.04)	204 (28) 125 (16) 108 (10) 238 (17)
	Crustaceans	1.6 (0.02) 1.5 (0.02) 1.7 (0.02) 1.7 (0.03)	55 (9)	46 (6) 138 (14) 311 (64)
	Insects	6.2 (0.11) 6.7 (0.09)	7.8 (0.1)	6.8 (0.1)	83 (8)	79 (5) 150 (11) 335 (15)
	Molluscs	2.8 (0.07) 2.9 (0.06) 3.4 (0.05) 3.8 (0.07)	38 (8)	29 (2)	60 (6)	125 (12)

Table 3 .5

 3 Beta diversity of component communities between seasons in each period.

	Component communities

Table 3 .6

 3 List of important taxa that the key contributors to the β diversity of global communities. The most important taxa that always had a high contribution were in bold. sp: spring, sm: summer, at: autumn.

	Period

145 0.130 0.110 0.142 0.098 0.093 0.093 0.077 0.070 0.093 0.101 0.103

	A12 Lumbriculidae	0.017 0.006 0.005	-	0.005 0.007	0.007 0.005 0.008	-	-	-
	A13 Glossiphonia	0.012 0.011 0.011	0.015 0.015 0.016	0.016 0.020 0.018	0.011 0.016 0.011
	A15 Helobdella	0.023 0.029 0.025	0.025 0.032 0.036	0.024 0.031 0.037	0.018 0.032 0.027
	C04 Limnadiidae	-	-	-	-	-	-	-	-	-	0.024 -	-
	C05 Corophiidae	0.012 0.005 0.007	0.012 0.007 -	0.009 0.005 -	0.019 0.019 0.007
	C07 Gammaridae	0.

067 0.051 0.078 0.075 0.071 0.090 0.085 0.068 0.092 0.101 0.087 0.080

	C09 Palaemonidae	-	0.012 0.023	-	-	0.005	-	-	0.014	-	-	-
	C15 Asellidae	0.

070 0.076 0.072 0.105 0.075 0.071 0.053 0.076 0.081 0.070 0.096 0.083

	C18 Mysidae	-	0.021 0.034	-	0.014 0.017	-	0.008 0.018	-	0.010 -
	C21 Ostracoda	-	-	0.005	-	-	0.005	0.013 0.012 0.025	0.010 -	0.005
	I001 Dryopidae	-	-	-	-	-	-	-	-	-	0.006 0.006 0.005
	I002 Dytiscidae	0.020 0.013 0.014	0.007 0.013 0.011	0.012 0.010 0.008	-	-	-
	I005 Haliplidae	-	-	0.005	-	0.005 0.005	-	0.007 0.006	-	-	-
	I007 Hydrophilidae	-	-	-	-	0.006 -	0.006 0.010 -	-	-	-
	I014 Ceratopogonidae	0.005 -	-	-	-	-	0.018 -	-	0.011 -	-
	I016

Chironomidae. Non thummi-plumosus 0.076 0.075 0.065 0.084 0.065 0.064 0.064 0.060 0.050 0.076 0.074 0.070 I017 Chironomidae. Thummi-plumosus 0.076 0.105 0.079 0.095 0.082 0.090 0.080 0.094 0.066 0.141 0.120 0.128

	I018 Culicidae	-	0.008 0.009	-	0.006 0.008	-	0.009 0.009	-	0.014 0.014
	I025 Psychodidae	-	-	-	0.005 0.006 -	0.011 0.011 0.007	-	-	0.008
	I030 Simuliidae	0.008 0.015 0.010	0.019 0.030 0.031	0.033 0.030 0.033	0.043 0.042 0.028
	I035 Limoniidae	0.006 0.006 -	0.005 0.010 0.006	0.008 0.007 -	-	-	-
	I037 Baetis	0.021 0.017 0.009	0.027 0.041 0.030	0.039 0.034 0.023	0.022 0.036 0.020
	I039 Cloeon	0.005 0.020 0.032	0.007 0.023 0.022	0.010 0.017 0.021	0.013 0.008 0.019
	I042 Caenis	0.011 -	-	0.008 -	-	0.008 0.007 -	0.013 0.005 -
	I053 Potamopyrgus	0.010 0.018 0.014	0.019 0.017 0.019	0.019 0.018 0.031	0.022 0.023 0.030
	I061 Micronecta	-	-	-	0.007 0.010 -	0.007 0.006 -	0.008 0.011 -
	I063 Sigara	0.012 0.041 0.046	0.014 0.035 0.040	0.011 0.032 0.028	0.007 0.016 0.021
	I064 Gerris	-	-	-	-	-	-	-	0.006 0.007	-	-	0.008
	I081 Calopteryx	-	-	-	-	-	-	-	-	-	-	-	0.007
	I087 Ischnura	0.006 0.009 0.025	0.013 0.008 0.021	0.008 0.008 0.016	-	0.005 0.010
	I106 Nemoura	0.010 -	-	0.005 -	-	0.013 -	-	0.015 -	-
	I119 Hydropsychidae	0.006 -	-	-	0.008 0.007	0.012 0.009 0.013	0.018 0.008 0.016
	I120 Hydroptilidae	-	-	-	-	-	-	0.005 0.007 -	0.007 -	-
	I122 Leptoceridae	0.010 -	-	0.006 -	-	0.007 0.006 0.006	-	-	-
	I123 Limnephilidae	0.005 -	-	0.006 -	-	0.017 -	-	-	-	-
	M06 Dreissena	0.011 0.005 -	0.011 -	-	0.008 0.007 0.005	0.013 0.024 0.009
	M08 Pisidium	0.049 0.027 0.031	0.038 0.043 0.028	0.050 0.038 0.034	0.040 0.031 0.038
	M09 Sphaerium	0.010 0.014 0.012	0.010 0.010 0.007	0.012 0.014 0.009	0.009 0.008 0.010
	M12 Bithynia	0.015 0.018 0.013	0.023 0.017 0.019	0.019 0.022 0.020	0.015 0.023 0.019
	M13 Pseudamnicola	-	-	-	-	-	-	-	-	-	-	0.011 0.012
	M17 Lymnaea	0.045 0.027 0.032	0.021 0.032 0.026	0.031 0.031 0.020	0.015 0.014 0.020
	M20 Physa	0.015 0.032 0.044	0.014 0.034 0.048	0.010 0.025 0.032	-	-	-
	M21 Physella	-	-	-	-	-	-	-	0.008 0.012	0.013 0.023 0.046
	M23 Anisus	0.008 0.005 0.006	0.008 0.008 0.005	0.009 0.007 0.006	-	-	-
	M27 Gyraulus	-	0.008 0.008	0.006 0.008 0.009	0.006 0.008 0.010	-	0.008 0.012
	M31 Planorbis	-	-	0.007	-	0.007 -	-	-	-	-	-	-
	M33 Valvata	0.031 0.019 0.019	0.015 0.020 0.021	0.026 0.021 0.026	0.020 0.022 0.020

Table 3 .7

 3 Coefficient correlation and adjusted R 2 of multiple regression models between α diversity, global LCBD indices and environmental variables.

	Period Season		Most frequently significant variables
	Alpha diversity	NH4 +	NO 2 -	EC	NO 3 -	DO	Adj.R 2
	D1	Spring	-0.33*	-	-0.50***	-0.30**	-	0.27
		Summer	-0.46***	-	-0.16**	-	-	0.25
		Autumn	-0.15*	-	-	-	0.22*	0.18
	D2	Spring	-0.15*	-0.77**	-0.29***	-	0.23**	0.29
		Summer	-0.24***	-0.51***	-0.20***	-0.06*	0.22***	0.32
		Autumn	-0.09*	-	-0.20***	-	0.18**	0.19
	D3	Spring	-0.19**	-	-0.17**	-	0.18*	0.17
		Summer	-0.20***	-0.59***	-0.24***	-0.10**	0.18**	0.24
		Autumn	-0.21***	-0.65***	-0.22***	-0.08*	-	0.30
	D4	Spring	-0.39***	-	-0.14*	-	0.35**	0.31
		Summer	-0.21***	-	-0.20***	-0.10*	-	0.14
		Autumn	-0.21***	-	-0.19***	-0.09*	-	0.22

Global LCBD indices NH4 + NO 2 - EC pH Sinuosity Adj.R 2 D1

		Spring	-	-	-	0.0104***	-	0.10
		Summer	-	-	0.0002*	0.0016*	-0.001*	0.07
		Autumn	-	-	-	0.002**	-	0.09
	D2	Spring	-0.0002*	-	-	0.0026*	-0.0016**	0.05
		Summer	-0.0001***	-0.0003***	0.0001***	0.0013***	-0.0002*	0.12
		Autumn	-0.0001**	-	0.0001**	0.0012**	-0.0006**	0.11
	D3	Spring	-0.0002*	-	-	-	-	0.02
		Summer	-	-0.0003***	-	0.0007**	-0.0003***	0.06
		Autumn	-	-	0.0001***		-0.0003*	0.09
	D4	Spring	-	-	-	-	-	0.06
		Summer	-0.0001*	-0.001***	-	-	-0.0006***	0.11
		Autumn	-0.0001**	-0.0006***	-	0.0008*	-	0.07
	* p<0.05, ** p<0.01, *** p<0.001

Table 3 .8

 3 Coefficient correlation of multiple regression models between the global LCBD and LCBDs of component communities in different seasons of each period.

				LCBD indices of component communities
	Period	Season	LCBD A	LCBD C	LCBD I	LCBD M	Adj.R 2
	D1	Spring	0.24***	0.16***	0.33***	0.11*	0.64
		Summer	0.21***	0.09***	0.40***	0.02	0.66
		Autumn	0.15***	0.13***	0.32***	0.02	0.60
	D2	Spring	0.20***	0.11***	0.38***	0.10***	0.62
		Summer	0.18***	0.11***	0.40***	0.08***	0.72
		Autumn	0.16***	0.14***	0.39***	0.11***	0.67
	D3	Spring	0.12***	0.12***	0.48***	0.09***	0.66
		Summer	0.14***	0.11***	0.37***	0.06***	0.62
		Autumn	0.14***	0.17***	0.34***	0.05**	0.65
	D4	Spring	0.08***	0.20***	0.39***	0.07**	0.71
		Summer	0.16***	0.16***	0.36***	0.04*	0.72
		Autumn	0.13***	0.14***	0.39***	0.08***	0.64
	* p<0.05, ** p<0.01, *** p<0.001

Table 3 .9

 3 Results of linear regression models (y=a+b 1 x) and quadratic regression models (y=a+b 1 x+b 2 x 2 ) for the effects of the prevalence of macroinvertebrate species on the three performance measures (Kappa, AUC and ER) of the modelling techniques. Kappa: Cohen's Kappa Statistic, AUC: area under the curve, ER: error rate, LR: logistic regression, RF: random forest, ANN: artificial neural network, SVM: support vector machine. Asterisks indicate significance levels of regression coefficients.

	Measure Model	b 1	b 2	Adj.R 2	Measure Model	b 1	b 2	Adj.R 2
	LR					RF
	Kappa	Linear	0.006*		0.04	Kappa	Linear	0.012***	0.15
		Quadratic	0.613** -0.91***	0.12		Quadratic	1.226*** -1.29***	0.31
	AUC	Linear	0.003*		0.02	AUC	Linear	0.005***	0.15
		Quadratic	0.234	-0.46***	0.06		Quadratic	0.582*** -0.61***	0.31
	ER	Linear	0.008***		0.61	ER	Linear	0.009***	0.73
		Quadratic	0.861*** -0.23***	0.66		Quadratic	1.008*** -0.34***	0.81
	ANN					SVM
	Kappa	Linear	0.007**		0.04	Kappa	Linear	0.009***	0.12
		Quadratic	0.736** -1.07***	0.12		Quadratic	0.927*** -0.72***	0.19
	AUC	Linear	0.003**		0.03	AUC	Linear	0.004***	0.11
		Quadratic	0.357** -0.56***	0.11		Quadratic	0.347*** -0.35***	0.18
	ER	Linear	0.009***		0.59	ER	Linear	0.009***	0.73
		Quadratic	0.993*** -0.31***	0.66		Quadratic	1.032*** -0.40***	0.84
	*p<0.05, **p<0.01, ***p<0.001.

  List of sampling sites with theirLCBD indices, p-values and corrected p-values (cor-p). Corrected p-values that are significant at the 0.05 level are in bold (applied to the global communities). * indicates the hotspot location shared by different component communities. The first letter of each site name represents the country (T: Thailand, L: Laos, C: Cambodia, V: Vietnam); the second letter represents the channel of the river (M: Mekong, B: Bassac, T: Tributary). For example, TM1: site number 1 located on the Mekong in Thailand. R: richness (number of taxa).

	Taxonomic Taxonomic Taxonomic Taxonomic Taxonomic Appendix T3 Site	Global communities	Annelid communities	Species Species Species Species Species Crustacean communities	Mollusk communities	Insect communities
	group group group group group code	R	Order/Clade Order/Clade Order/Clade Order/Clade Order/Clade LCBD	p	Family Family Family Family Family cor-p	LCBD	Species Species Species Species Species p	cor-p	LCBD	p	occurrence occurrence occurrence occurrence occurrence cor-p	LCBD	p	cor-p	LCBD	p	cor-p
	Taxonomic group Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Annelid Crustacean Amphipoda (10) Order/Clade Aphanoneura Haplotaxida (15) --------------Lumbriculida Phyllodocida (4) ---Scolecida (3) --Spionida (2) -Unknown (4) ---Crustacean -Crustacean -Crustacean -Crustacean -Crustacean -Crustacean -Crustacean Isopoda (3) Crustacean -Crustacean -Crustacean Tanaidacea (7) Crustacean -Crustacean -Crustacean -Crustacean -Crustacean -Crustacean -Insect Coleoptera (9) Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect Diptera (36) Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect Ephemeroptera (32) Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect Hemiptera (8) Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect Plecoptera (3) Insect -Insect -Insect Trichoptera (21) Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Insect -Mollusc Arcoida (2) Mollusc -Mollusc Mytiloida (3) Mollusc -Mollusc -Mollusc Unionida (18) Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc -Mollusc Caenogastropoda (50) Ampullariidae (3) Family Aeolosomidae Naididae (15) --------------Lumbriculidae Nephtyidae Nereididae (3) --Opheliidae Orbiniidae (2) -Spionidae (2) -Unknown (4) ---Corophiidae (7) --Parathelphusidae (2) Parathelphusidae sp. Species Aeolosoma bengalense Aulodrilus prothecatus Branchiodrilus semperi Branchiura sowerbyi Chaetogaster langi Chaetogaster limnaei limnaei Chaetogaster sp. Dero denticulata Dero pectinata Dero sp. Dero sp.1 Dero sp.2 Dero sp.3 Limnodrilus hoffmeisteri Naididae sp. Pristina sp. Lumbriculidae sp. Nephtys polybranchia Namalycastis abiuma Namalycastis longicirris Neanthes caudata Polyophthalmus pictus Orbinia johnsoni Scoloplos sp. Polydora sp. Prionospio sp. Oligochaeta sp. Polychaeta sp. Polychaeta sp.1 Polychaeta sp.2 Corophium intermedium Corophium minutum Palaemonidae larva -Somanniathelphusa germaini Potamidae Ranguna cochinchinensis Unknown Decapoda larva Anthuridae (2) Cyathura carinata -Cyathura truncata Corallanidae Tachaea chinensis Apseudidae (2) Apseudes vietnamensis -Apseudidae sp. Gammaridae Gammarus sp. Haustoriidae (2) Eohaustorius sp. -Eohaustorius tandeensis Hyalidae (2) Hyale hawaiensis -Hyale sp. Amphizoidae Amphizoa sp. Dytiscidae (2) Hyphydrus sp. -Neptosternus sp. Elmidae (4) Cloeon sp. -Elmidae sp. -Heterlimnius sp. -Stenelmis sp. Haliplidae (2) Haliplidae sp. -Haliplus sp. Ceratopogonidae (3) Bezzia sp. -Culicoides sp. -Dasyhelea sp. Chaoboridae Chaoborus sp. Chironomidae (23) Ablabesmyia sp. -Chironomidae sp. -Chironomus attenuatus -Chironomus sp. -Chironomus sp.1 -Clinotanypus sp. -Clypeocaetis sp. -Cricotopus sp. -Cryptochironomus sp. -Tipulidae sp. Baetidae (6) Acentrella sp. -Baetiella sp. -Baetis sp. -Centroptilum sp. -Heterocloeon sp. -Procloeon sp. Behningiidae Behningiidae sp. Caenidae (4) Caenis sp. -Caenoculis sp. -Caenodes sp. -Cercobrachys sp. Ephemeridae (4) Afromera sp. -Eatonigenia sp. -Ephemera sp. -Hexagenia sp. Heptageniidae (6) Cinygmina sp. -Cladopelma sp. -Epeorus sp. -Heptagenia sp. -Heptageniidae sp. -Thalero sphyrus sp. Leptoplebiidae (3) Choropterpes sp. -Leptophlebia sp. -Traverella sp. Oligoneuriidae (2) Chromarcys sp. -Pentagenia sp. Palingeniidae (2) Anagenesia sp. -Palingeniidae sp. Polymitarcyidae (2) Ephoron sp. -Povilla sp. Potamanthidae Potamanthus sp. Prosopistomatidae Prosopistoma sp. Belostomatidae Diplonychus rusticus Corixidae (4) Micronecta sp. -Corixa sp. -Microtendipes sp. -Ophiogomphus sp. -Orientogomphus sp. -Progomphus sp. Libellulidae Libellula sp. Petaluridae Tachopteryx sp. Perlidae (3) Etrocorema sp. -Perla sp. -Phanoperla sp. Dipseudopsidae Dipseudopsis sp. Ecnomidae Economus sp. Glossosomatidae Glososoma sp. Goeridae Goera sp. Hydropsychidae (3) Cheumaatopsyche sp. -Hydropsyche sp. -Macronema sp. Hydroptilidae (3) Agraylea sp. -Orthotrichia sp. -Oxyethira sp. Leptoceridae (4) Leptocerus sp. -Leptoceridae sp. -Oecetis sp. -Nectopsyche sp. Limnephilidae (2) Farula sp. -Limnephilus sp. Molannidae Molanna sp. Philopotamidae (2) Chimarra sp. -Philopotamidae sp. Psychomyiidae Psychomyiidae sp. Rhyacophilidae Rhyacophila sp. Arcidae (2) Scaphula pinna -Scaphula sp. Mytilidae (3) Limnoperna siamensis -Limnoperna sp. -Sinomytilus harmandi Unionidae (18) Ensidens ingallsianus ingallsianus -Hyriopsis Hyriopsis bialatus -Hyriopsis Limnoscapha desowitzi -Corbicula cyreniformis -Corbicula fluminea -Corbicula gustaviana -Corbicula lamarckiana -Corbicula larva -Corbicula leviuscula -Corbicula moreletiana -Corbicula sp. -Corbicula tenuis Pisidiidae Afropisidium clarkeanum Pila ampullacea Mollusc --Pila polita Mollusc --Pila scutata Mollusc -Assimineidae (2) Cyclotropis bollingi Mollusc --Cyclotropis sp. Mollusc -Bithyniidae (3) Bithynia siamensis Mollusc --Bithynia sp. Mollusc --Wattebledia siamensis Mollusc -Buccinoidae (3) Clea scalarina Mollusc --Clea helena Mollusc --Clea sp.1 Mollusc -Cochliopidae Cochliopa riograndensis Mollusc -Hydrobiidae Kareliania sp. Mollusc -Pachychilidae (2) Adamietta housei Mollusc --Brotia sp. Mollusc -Pomatiopsidae (8) Hubendickia crooki Mollusc --Hubendickia sp. Mollusc --Hydrorissoia sp. Mollusc --Lacunopsis sp. Mollusc --Pachydrobia brevis Mollusc --Pachydrobia sp. Mollusc --Pachydrobiella sp. Mollusc --Paraprososthenia sp. Mollusc -Stenothyra annandalei Mollusc -Stenothyra glabrata Mollusc -Stenothyra jiraponi Mollusc --Stenothyra koratensis holosculpta LT6 31 0.0103 1 1 0.0113 0.875 1 0.0012 -LT5 6 0.0113 1 1 0.0052 0.964 1 0.0012 -LT4 25 0.0118 0.998 1 0.0128 0.786 1 0.0012 Stenothyridae (9) Mollusc --Trochotaia trochoides Mollusc Heterobranchia (2) LT7 28 0.0121 0.997 1 0.0126 0.800 1 0.0422 1E-04 0.006 Species occurrence 29 11 10 28 20 10 12 18 30 46 25 33 11 25 11 21 14 20 17 10 12 11 20 10 14 14 11 42 12 10 18 1 1 1 1 1 1 Ellobiidae (2) Melampus nucleolus Mollusc -LT8 31 0.0152 0.669 1 0.0162 0.435 1 0.0394 1E-04 0.006 -Melampus fasciatulus Mollusc Hygrophila (3) LT9 12 0.0137 0.911 1 0.0052 0.966 1 0.0012 1 1 Lymnaeidae (2) Lymnaea swinhoei Mollusc --Lymnaea viridis LT10 19 0.0175 0.156 1 0.0244 0.009 0.561 0.0012 1 1 11 Mollusc -Planorbidae Gyraulus sp. Mollusc Heterostropha Pyramidellidae Morrisonietta spiralis Mollusc Neritimorpha (2) Neritidae (2) Neritina rubida Mollusc --Neritina violacea Mollusc Unknown (2) Unknown (2) Gastropoda larva Mollusc --Gastropoda sp. Site code R Global communities Annelid communities Crustacean communities LT11 41 0.0138 0.915 1 0.0198 0.201 1 0.0012 1 1 LCBD p cor-p LCBD p cor-p LCBD p cor-p LT12 32 0.0152 0.656 1 0.0169 0.387 1 0.0012 1 1 CB1 30 0.0165 0.337 1 0.0153 0.516 1 0.0379 1E-04 0.006 TM1 9 0.0200 1E-04 0.012 0.0192 0.236 1 0.0012 1 1 CM1 6 0.0185 0.041 1 0.0165 0.414 1 0.0012 1 1 TM2 21 0.0210 1E-04 0.006 0.0192 0.239 1 0.0379 1E-04 0.006 CM2 8 0.0174 0.166 1 0.0167 0.400 1 0.0012 1 1 TM3 12 0.0158 0.519 1 0.0165 0.414 1 0.0012 1 1 CM3* 12 0.0228 1E-04 0.006 0.0256 5E-04 0.032 0.0376 1E-04 0.006 TM4 12 0.0209 1E-04 0.006 0.0192 0.236 1 0.0012 1 1 CM4 16 0.0141 0.874 1 0.0120 0.833 1 0.0354 1E-04 0.006 TM5 12 0.0146 0.785 1 0.0165 0.416 1 0.0012 1 1 CT1 11 0.0114 1 1 0.0165 0.422 1 0.0012 1 1 TT1 18 0.0165 0.353 1 0.0158 0.472 1 0.0012 1 1 CT2 15 0.0124 0.993 1 0.0165 0.411 1 0.0012 1 1 TT2 10 0.0153 0.630 1 0.0165 0.412 1 0.0012 1 1 CT3 26 0.0132 0.968 1 0.0117 0.853 1 0.0368 1E-04 0.006 TT3 16 0.0176 0.136 1 0.0221 0.107 1 0.0012 1 1 CT4 10 0.0149 0.737 1 0.0167 0.411 1 0.0012 1 1 TT4 8 0.0133 0.953 1 0.0139 0.682 1 0.0012 1 1 CT5 16 0.0142 0.852 1 0.0116 0.858 1 0.0012 1 1 TT5 10 0.0200 0.001 0.065 0.0192 0.242 1 0.0012 1 1 CT6 6 0.0130 0.978 1 0.0165 0.414 1 0.0012 1 1 TT6 23 0.0165 0.341 1 0.0192 0.233 1 0.0379 1E-04 0.006 CT7 25 0.0172 0.202 1 0.0131 0.766 1 0.0385 1E-04 0.006 TT7 27 0.0144 0.823 1 0.0130 0.775 1 0.0012 1 1 CT8 16 0.0191 0.014 0.756 0.0125 0.814 1 0.0394 1E-04 0.006 TT8 17 0.0188 0.027 1 0.0192 0.239 1 0.0012 1 1 CT9 11 0.0122 0.997 1 0.0165 0.422 1 0.0012 1 1 VB1 74 0.0162 0.420 1 0.0190 0.252 1 0.0341 1E-04 0.006 CT10 17 0.0100 1 1 0.0119 0.849 1 0.0012 1 1 VB2 18 0.0168 0.286 1 0.0128 0.791 1 0.0353 1E-04 0.006 CT11 25 0.0112 1 1 0.0086 0.960 1 0.0394 1E-04 0.006 VB3 52 0.0195 0.007 0.226 0.0196 0.210 1 0.0354 1E-04 0.006 CT12 29 0.0139 0.894 1 0.0135 0.732 1 0.0386 1E-04 0.006 VB4 34 0.0204 5E-04 0.012 0.0218 0.114 1 0.0388 1E-04 0.006 LM1 39 0.0159 0.484 1 0.0118 0.851 1 0.0428 1E-04 0.006 VB5 66 0.0157 0.533 1 0.0189 0.247 1 0.0356 1E-04 0.006 LM2 16 0.0150 0.707 1 0.0147 0.585 1 0.0012 1 1 VM1 53 0.0189 0.017 0.896 0.0142 0.639 1 0.0360 1E-04 0.006 LM3 19 0.0148 0.743 1 0.0140 0.674 1 0.0012 1 1 VM2 67 0.0177 0.128 1 0.0186 0.269 1 0.0360 1E-04 0.006 LM4 7 0.0146 0.790 1 0.0052 0.965 1 0.0012 1 1 VM3 19 0.0166 0.318 1 0.0135 0.728 1 0.0400 1E-04 0.006 LM5 24 0.0129 0.980 1 0.0116 0.859 1 0.0012 1 1 VM4 33 0.0193 0.009 0.468 0.0212 0.140 1 0.0390 1E-04 0.006 LM6 23 0.0186 0.034 1 0.0192 0.224 1 0.0424 1E-04 0.006 VM5 61 0.0176 0.126 1 0.0180 0.310 1 0.0358 1E-04 0.006 LM7 8 0.0131 0.972 1 0.0167 0.405 1 0.0012 1 1 VM6 49 0.0190 0.002 0.104 0.0192 0.242 1 0.0012 1 1 LT1 41 0.0154 0.615 1 0.0147 0.582 1 0.0012 1 1 VT1 19 0.0190 0.017 0.896 0.0117 0.856 1 0.0012 1 1 LT2 24 0.0197 0.004 0.205 0.0208 0.159 1 0.0012 1 1 VT2 10 0.0140 0.882 1 0.0158 0.467 1 0.0012 1 1 LT3 15 0.0150 0.718 1 0.0244 0.009 0.527 0.0428 1E-04 0.006 VT3 8 0.0165 0.344 1 0.0165 0.415 1 0.0012 1 1	0.0105 0.0107 0.0119 0.0185 0.0136 0.0042 0.0185 Mollusk communities 1 1 0.999 1 0.965 1 0.104 1 0.825 1 1 1 0.098 1 0.0154 0.618 1 LCBD p cor-p 0.0175 0.247 1 0.0152 0.646 1 0.0042 1 1 0.0142 0.771 1 0.0216 1E-04 0.006 0.0224 1E-04 0.006 0.0107 0.999 1 0.0228 1E-04 0.006 0.0207 0.001 0.065 0.0144 0.752 1 0.0130 0.877 1 0.0115 0.986 1 0.0198 0.016 0.750 0.0134 0.850 1 0.0211 2E-04 0.011 0.0123 0.933 1 0.0102 0.999 1 0.0131 0.875 1 0.0218 1E-04 0.006 0.0169 0.358 1 0.0208 9E-04 0.047 0.0231 1E-04 0.006 0.0202 0.007 0.312 0.0228 1E-04 0.006 0.0137 0.831 1 0.0164 0.427 1 0.0208 0.001 0.051 0.0110 0.997 1 0.0135 0.834 1 0.0107 0.999 1 0.0167 0.384 1 0.0118 0.969 1 0.0197 0.025 1.000 0.0097 1 1 0.0189 0.069 1 0.0189 0.065 1 0.0141 0.78 1 0.0155 0.602 1 0.0195 0.032 1 0.0135 0.843 1 0.0160 0.519 1 0.0107 1 1 0.0125 0.92 1 0.0201 0.007 0.320 0.0186 0.087 1 0.0209 4E-04 0.022 0.0180 0.165 1 0.0042 1 1 0.0202 0.006 0.27 0.0197 0.022 1.000 0.0208 5E-04 0.027 0.0224 1E-04 0.006 0.0107 0.999 1 0.0129 0.891 1 0.0214 1E-04 0.006	0.0087 0.0101 0.0116 0.0100 0.0138 0.0112 0.0176 0.0130 Insect communities 0.999 1 0.986 1 0.918 1 0.99 1 0.716 1 0.949 1 0.29 1 0.806 1 LCBD p 0.0152 0.558 1 cor-p 0.0186 0.201 0.0233 0.012 0.702 1 0.0237 0.0238 0.005 0.317 0.008 0.474 0.0098 0.991 0.0165 0.406 1 1 0.0265 0.0207 0.073 1 1E-04 0.006 0.0158 0.477 0.0153 0.542 1 1 0.0108 0.967 0.0179 0.262 1 1 0.0120 0.891 0.0146 0.625 1 1 0.0136 0.733 0.0206 0.075 1 1 0.0167 0.378 0.0109 0.962 1 1 0.0154 0.535 0.0214 0.051 1 1 0.0103 0.982 0.0162 0.444 1 1 0.0146 0.623 0.0105 0.973 1 1 0.0231 0.0199 0.119 1 0.012 0.696 0.0121 0.878 0.0189 0.182 1 1 0.0088 1 0.0181 0.239 1 1 0.0112 0.939 0.0204 0.095 1 1 0.0123 0.866 0.0248 0.002 0.124 1 0.0146 0.62 0.0157 0.495 1 1 0.0158 0.48 0.0197 0.129 1 1 0.0154 0.528 0.0175 0.301 1 1 0.0139 0.713 0.0144 0.648 1 1 0.0083 1 0.0218 0.038 1 1 0.0190 0.161 0.0185 0.212 1 1 0.0115 0.932 0.0218 0.04 1 1 0.0124 0.868 0.0148 0.599 1 1 0.0174 0.309 0.0150 0.568 1 1 0.0152 0.55 1 0.0168 0.368 1
	Crustacean Insect Insect Mollusc Mollusc	-----			-----		Corophium sp. Diamesinae sp. Sigara sp. Indonaia pilata Stenothyra koratensis koratensis
	Crustacean Insect Insect Mollusc Mollusc	-----			--Delphacidae --	Grandidierella lignorum Einfeldia sp. Megamelus sp. Physunio cambodiensis Stenothyra labiata
	Crustacean Insect Insect Mollusc Mollusc	-----			--Gerridae --		Grandidierella vietnamica Glyptotendipes sp. Gerridae sp. Physunio micropterus Stenothyra mcmulleni		11 12
	Crustacean Insect Insect Mollusc Mollusc	-----			--Naucoridae --	Kamaka sp. Goeldichironomus sp. Naucoris sp. Pilsbryoconcha exilis compressa Stenothyra moussoni		27
	Crustacean Insect Insect Mollusc Mollusc	--Lepidoptera --		--Pyralidae --		Monocorophium sp. Orthocladiinae sp. Pyralidae sp. Pilsbryoconcha exilis exilis Stenothyra sp.
	Crustacean Insect Insect Mollusc Mollusc	--Odonata (22) --		Melitidae (2) -Aeshnidae -Thiaridae (6)	Melita sp. Parachironomus sp. Aeshna sp. Pilsbryoconcha lemeslei Melanoides sp.		11
	Crustacean Insect Insect Mollusc Mollusc	-----			--Calopterygidae (2) Agrion sp. Melita vietnamica Polypedilum sp. -Pseudodon cambodjensis -Melanoides tuberculata		44 12
	Crustacean Insect Insect Mollusc	----			Oedicerotidae ---	Perioculodes sp. Procladius sp. Calopteryx sp. cambodjensis Neoradina prasongi
	Crustacean Cumacea Insect -Insect -Mollusc -Mollusc -		Diastylidae -Corduliidae (3) --	Diastylis rathkei Pseudochironomus sp. Cordulia sp. Pseudodon inoscularis cumingi Sermyla tornatella		16
	Crustacean Decapoda (18) Insect -Insect -Mollusc -Mollusc -		Alpheidae ----		Alpheus bisincisus Pseudodiamesa sp. Epitheca sp. Pseudodon vondembuschianus Tarebia granifera
	Crustacean Insect Insect Mollusc	----			Atyidae (3) ---	Caridina nilotica Sergentia sp. Macromia sp. ellipticus Thiara scabra
	Crustacean Insect Insect Mollusc Mollusc	-----			--Gomphidae (14) -Viviparidae (12)	Caridina sp. Smittia sp. Amphylla williamsoni Scabies scobinata Angulyagra polyzonata
	Crustacean Insect Insect Mollusc Mollusc	-----			-----		Caridina verrata verrata Tanypus sp. Aphylla sp. Scabies sp. Angulyagra sp.		11 10
	Crustacean Insect Insect Mollusc Mollusc	-----			Palaemonidae (10) Macrobrachium dienbienphuensis -Tanytarsus sp. -Arigomphus sp. -Trapezoideus exolescens comptus -Anulotaia sp.
	Crustacean Insect Insect Mollusc Mollusc	-----			-Dolichopodidae ---	Macrobrachium equidens Hydrophorus sp. Burmagomphus sp. Uniandra contradens ascia Filopaludina Filopaludina doliaris
	Crustacean Insect Insect Mollusc Mollusc	-----			-Limoniidae ---	Macrobrachium lanchesteri Limnophila sp. Dromogomphus sp. Uniandra contradens tumidula Filopaludina Filopaludina filosa		23
	Crustacean Insect Insect Mollusc Mollusc	-----			-Tabanidae (2) ---	Macrobrachium mekongene Chrysops sp. Gastrogomphus sp. Unionida sp. Filopaludina sp.
	Crustacean Insect Insect Mollusc Mollusc	---Veneroida (15) -		---Corbiculidae (14) -	Macrobrachium pilimanus Tabanus sp. Gomphidae sp. Corbicula arata Mekongia sp.
	Crustacean Insect Insect Mollusc Mollusc	-----			-Tipulidae (5) ---	Macrobrachium rosenbergii Antocha sp. Gomphus sp. Corbicula baudoni Mekongia swainsoni braueri		12
	Crustacean Insect Insect Mollusc Mollusc	-----			-----		Macrobrachium secamense Eriocera sp. Labrogomphus sp. Corbicula blandiana Mekongia swainsoni flavida		15 16
	Crustacean Insect Insect Mollusc Mollusc	-----			-----		Macrobrachium sp. Pedicia sp. Megalogomphus sp. Corbicula bocourti Mekongia swainsoni swainsoni
	Crustacean Insect Insect Mollusc Mollusc	-----			-----		Palaemon curvirostris Tipula sp. Octogomphus sp. Corbicula castanea Sinotaia aeruginosa

Table 1

 1 Mean value (and standard deviation) for richness, abundance and Shannon diversity of macroinvertebrate assemblage and for environmental variables in each cluster.

	Cluster [n]

WT: water temperature, DO: dissolved oxygen, WC: water conductivity, SAW: the surface area of watersheds, Mix_evg. dec: mixed evergreen and deciduous forests. The number of samples [n] in each cluster is indicated between square brackets. * indicates ANOVA and Kruskal-Wallis Test for significant differences among clusters at p < 0.05. Superscripts (IIa, IIb1, IIb2) indicate significant pair-wise comparisons between the corresponding cluster (each column) and superscript-labeled clusters (i.e. IIa, IIb1, IIb2) at p < 0.05.

Table 2

 2 List of indicator taxa (and their indicator values, IndVal) of macroinvertebrate assemblage in each cluster.

				Table 2 (continued)
				Insect
	Cluster I	IndVal	p-value	Anagenesia sp. Caenoculis sp.	0.67 0.52	0.005 0.01
	Annelid Aeolosoma bengalense Aulodrilus prothecatus Chaetogaster langi Chaetogaster limnaei limnaei	0.52 0.67 0.85 0.6	0.01 0.005 0.005 0.005	Caenodes sp. Choropterpes sp. Dipseudopsis sp. Heterocloeon sp. Micronecta sp.	0.74 0.51 0.69 0.52 0.6	0.005 0.005 0.005 0.01 0.005
	Dero pectinata	0.67	0.005
	Dero sp.	0.74	0.005	Cluster IIb1	IndVal	p-value
	Dero sp.1	1	0.005	Mollusc
	Dero sp.2	0.95	0.005	Filopaludina filopaludina filosa	0.45	0.025
	Lumbriculidae sp. Namalycastis longicirris Nectopsyche sp.	0.6 0.9 0.67	0.005 0.005 0.005	Insect Pentagenia sp.	0.62	0.01
	Orbinia johnsoni	0.52	0.01
	Polydora sp.	0.67	0.005	Cluster IIb2	IndVal	p-value
	Crustacean Corophium minutum	0. 8	0.005	Annelid Naididae sp.	0.76	0.005
	Corophium sp.	0.67	0.005	Insect
	Cyathura carinata	0.74	0.005	Gomphidae sp.	0.56	0.01
	Cyathura truncata	0.57	0.005
	Decapoda sp.	0.91	0.005
	Eohaustorius sp.	0.6	0.005
	Eohaustorius tandeensis	0.67	0.005
	Gammarus sp.	0.6	0.005
	Grandidierella lignorum	0.78	0.005
	Grandidierella vietnamica	1	0.005
	Hyale hawaiensis	0.67	0.005
	Hyale sp.	0.85	0.005
	Kamaka sp.	0.6	0.005
	Macrobrachium equidens	0.6	0.005
	Melita sp.	0.82	0.005
	Palaemon curvirostris	0.6	0.005
	Mollusc
	Afropisidium clarkeanum	0.73	0.005
	Angulyagra polyzonata	0.6	0.005
	Angulyagra sp.	0.9	0.005
	Bithynia siamensis	0.67	0.005
	Corbicula baudoni	0.87	0.005
	Corbicula bocourti	0.74	0.005
	Corbicula leviuscula	0.97	0.005
	Corbicula moreletiana	0.86	0.005
	Corbicula sp.	0.95	0.005
	Gastropoda sp.	0.74	0.005
	Hyriopsis bialatus	0.64	0.005
	Limnoperna siamensis	0.99	0.005
	Limnoperna sp.	0.95	0.005
	Lymnaea viridis	0.94	0.005
	Mekongia swainsoni swainsoni	0.67	0.005
	Sinomytilus harmandi	0.9	0.005
	Stenothyra annandalei	0.6	0.005
	Stenothyra glabrata	0.85	0.005
	Trochotaia trochoides	0.52	0.005
	Insect
	Arigomphus sp.	0.67	0.005
	Cricotopus sp.	1	0.005
	Clinotanypus sp.	0.52	0.005
	Monocorophium sp.	0.91	0.005
	Sigara sp.	0.6	0.005
	Cluster IIa	IndVal	p-value
	Annelid
	Oligochaeta sp.	0.99	0.005
	Polychaeta sp.1	0.6	0.005
	Mollusc
	Corbicula sp.	0.88	0.005
	Hubendickia sp.	0.6	0.01
	Kareliania sp.	0.52	0.01
	Scaphula sp.	0.52	0.01
	Stenothyra sp.	0.6	0.005

Table 1

 1 Observed environmental factors, richness (number of taxa) and abundance (number of individuals) of macroinvertebrates across the 63 sampling sites.

	Variables	Unit	Min	Max	Mean	Standard deviation
	Altitude	m	3	546	127	132
	Water temperature	° C	17	31	27	3
	Dissolved oxygen	mg/L	2.7	9.3	7.4	2.6
	Water conductivity	mS/m	3.9	66.6	17.9	11.9
	River width	m	11	1629	467	466
	River depth	m	0.4	15.0	4.8	3.9
	Secchi depth	m	0.2	3.0	0.9	0.6
	pH	-	6.8	8.4	7.6	0.6
	Richness	taxa/sample	6	74	23	16
	Abundance	individuals/sample	13	1997	315	396

  .3 in Part I: Synthesis). The SCBD values are small because the SCBD indices are relative to the total sum of squares in the community composition table and sum to 1. SCBD indices indicate taxa that have the highest variance across sites. The 3 highest SCBD indices belonged to insect taxa: Polypedilum sp. (0.054), Ablabesmyia sp. (0.039), Cryptochironomus sp. (0.037), followed by Corbicula tenuis (mollusk, 0.037), Goeldichironomus sp. (insect, 0.035) and Corbicula leviuscula (mollusk, 0.034).

Table 2

 2 Pearson correlation coefficients between the normalized geographical factors and richness and abundance of each component community among the important taxa. R: richness (number of taxa), A: abundance (number of individuals), LONG: longitude (m), LAT: latitude (m), ALT: altitude (m). Significant relationships are marked with stars.* p < 0.05, ** p < 0.01, *** p < 0.001

		Annelids		Crustaceans	Mollusks		Insects
		R	A	R	A	R	A	R	A
	LONG 0.15	0.06	0.03	0	-0.01	-0.04	-0.03	-0.05
	LAT	-0.58**	-0.36**	-0.70*** -0.25*	-0.65***	-0.35**	0.47*** 0.33**
	ALT	-0.65*** -0.42***	-0.65*** -0.26*	-0.69***	-0.35**	0.41*** 0.25*

Table 3

 3 Contribution of geo-environmental factors to the variation of LCBD indices of global macroinvertebrate communities. The model showing the strongest contribution is in bold.

		Adjusted		Stepwise		p of the
	Model	R 2	AIC	selection	F-stat	model
	Environmental factors
	LCBD ~ WT	--	-551.2	--	0.01	0.956
	LCBD ~ DO	--	-551.5	--	0.3	0.609
	LCBD ~ WC	0.135	-561.4	--	10.7	0.002
	LCBD ~ RW	0.010	-552.9	--	1.6	0.206
	LCBD ~ RD	0.120	-560.2	--	9.5	0.003
	LCBD ~ SD	0.061	-556.2	--	5.1	0.029
	LCBD ~ pH	--	-551.5	--	0.3	0.590
	LCBD ~ all factors	0.235	-563.6	WC+RD+SD	3.7	0.002
	LCBD ~ WC+RD+SD	0.260	-569.4	WC+RD+SD	8.3	<0.001
	Geographical factors
	LCBD ~ LONG	--	-551.8	--	0.6	0.449
	LCBD ~ LAT	0.001	-552.3	--	1.0	0.314
	LCBD ~ ALT	0.010	-552.9	--	1.6	0.207
	LCBD ~ all factors	--	-549.7	--	0.8	0.511

WT: water temperature, DO: dissolved oxygen, WC: water conductivity, RW: river width, RD: river depth, SD: Secchi depth, LONG: longitude (m), LAT: latitude (m), ALT: altitude (m).

Table 4

 4 Contribution of the composition of component communities to the LCBD indices of the global macroinvertebrate communities. The model showing the strongest contribution is in bold. LCBDA, LCBDC, LCBDM, LCBDI are the LCBD indices of each component community (annelids, crustaceans, mollusks and insects, respectively).

	Model

Table 1

 1 Minimum, maximum and mean values and standard deviation for environmental predictors across the 63 sampling sites.

	Predictor	Unit	Min	Max	Mean	Standard deviation
	Altitude	m	3	546	127	132
	Water temperature	• C	16.7	31	27.4	3
	Dissolved oxygen	mg/L	2.7	9.3	7.1	1.3
	Water conductivity	mS/m	3.9	66.6	17.9	11.9
	Distance from the sea	km	82	2597	998	686
	River width	m	11	1629	466	466
	River depth	m	0.4	15	4.8	3.9
	Secchi depth	m	0.2	3	0.9	0.6
	pH	-	6.8	8.4	7.6	0.6

  of the modelling techniques. Kappa: Cohen's Kappa Statistic, AUC: area under the curve, ER: error rate, LR: logistic regression, RF: random forest, ANN: artificial neural network, SVM: support vector machine. Asterisks indicate significance levels of regression coefficients.

	Measure	Regression Model	ˇ1	ˇ2	Adjusted R 2
	LR
	Kappa	Linear	0.006 *		0.04
		Quadratic	0.613 **	-0.91 ***	0.12
	AUC	Linear	0.003 *		0.02
		Quadratic	0.234	-0.46 ***	0.06
	ER	Linear	0.008 ***		0.61
		Quadratic	0.861 ***	-0.23 ***	0.66
	RF
	Kappa	Linear	0.012***		0.15
		Quadratic	1.226 ***	-1.29 ***	0.31
	AUC	Linear	0.005 ***		0.15
		Quadratic	0.582 ***	-0.61 ***	0.31
	ER	Linear	0.009 ***		0.73
		Quadratic	1.008 ***	-0.34 ***	0.81
	SVM
	Kappa	Linear	0.009 ***		0.12
		Quadratic	0.927 ***	-0.72 ***	0.19
	AUC	Linear	0.004 ***		0.11
		Quadratic	0.347 ***	-0.35 ***	0.18
	ER	Linear	0.009 ***		0.73
		Quadratic	1.032 ***	-0.40 ***	0.84
	ANN
	Kappa	Linear	0.007 **		0.04
		Quadratic	0.736 **	-1.07 ***	0.12
	AUC	Linear	0.003 **		0.03
		Quadratic	0.357 **	-0.56 ***	0.11
	ER	Linear	0.009 ***		0.59
		Quadratic	0.993 ***	-0.31 ***	0.66
	* p < 0.05.
	** p < 0.01.
	*** p < 0.001.

Table 2 .

 2 List of alien and native molluscs, their occurrences and abundance recorded for each period. The total number of samples for each period is indicated in brackets.

		Corbicula	Dreissena	Ferrisia	Lithoglyphus	Menetus	Physella	Viviparus Merging all genera
	Instances	130	745	381	30	4	1138	94	2028
	Instances (%)	1.4	8.8	4.5	0.4	<0.1	13.9	1.1	26.4
				Occurrence				Abundance

Table 3 .

 3 Mean (and standard deviation) for environmental predictors and occurrence instances of each class of the response variable. NH 4 + : Ammonium, COD: Chemical Oxygen Demand, Pt: Total Phosphorus, EC: Electrical Conductivity, NO 3 Dissolved Oxygen. The number of sampled sites for each period is shown in square brackets.

	-: Nitrate, NO 2	-: Nitrite, oPO 4 :

Table 2 .

 2 Mean value (and standard deviation) for environmental predictors, and occurrence instances of each class of the response variable. NH4 + : Ammonium, COD: Chemical Oxygen Demand, Pt: Total Phosphorus, EC: Electrical Conductivity, NO3 -: Nitrate, NO2 -: Nitrite, oPO4: Orthophosphate, DO: Dissolved Oxygen, SI: Sinuosity.

									Grouping
	Genus	Corbicula Dreissena Ferrisia Lithoglyphus Menetus Physella Viviparus	all genera
	Instances	130	745	381	30	4	1138	94	2028
	Instances (%)	1.4	8.8	4.5	0.4	<0.1	13.9	1.1	26.4
						Periods
	Variable	Unit	D1: 1991-1995	D2: 1996-2000	D3: 2001-2005	D4: 2006-2010
	NH4 +	mg/L	2.6 (4.8)		1.9 (3.4)		1.6 (2.8)	2.3 (5.1)
	COD	mg/L	55 (43)		39 (42)		34 (34)	36 (41)
	Pt	mg/L	1.0 (1.4)		0.9 (1.4)		0.9 (2.3)	0.8 (1.0)
	EC	µS/cm	1320 (2618)		987 (1149)		998 (1509)	921 (949)
	NO3 -	mg/L	3.3 (4.1)		4.0 (4.7)		3.6 (3.7)	3.0 (3)
	NO2 -	mg/L	0.2 (0.3)		0.2 (0.3)		0.2 (0.2)	0.2 (0.2)
	oPO4	mg/L	0.6 (0.9)		0.5 (1.1)		0.4 (0.7)	0.5 (0.8)
	pH		7.5 (0.6)		7.6 (0.5)		7.7 (0.4)	7.6 (0.4)
	DO	mg/L	7.5 (3.5)		6.8 (3)		6.9 (3.1)	6.6 (2.9)
	SI		1.1 (0.1)		1.1 (0.1)		1.1 (0.1)	1.1 (0.1)
	Slope	m/1000m	1.0 (1.4)		1.6 (2.8)		2.0 (3.6)	1.7 (2.8)
	Response	class
	Native		863		1860		2259		842
	Alien		6		15		63		133
	Overlap		66		146		567		875
	Total instances	935		2021		2889	1850

Acknowledgements

Firstly, I would like to express my sincere thanks and profound gratitude to my thesis directors Prof Dr Sovan LEK, Prof Dr Peter GOETHALS and Dr Pieter BOETS for their scientific support. They have a huge pile of practical, ecological and theoretical knowledge. Their guidance, advice, suggestions and encouragements allowed me to survive during this scientific journey. I thank you all very much for your precious time spent on reading and correcting my work and discussing with me with an open mind. I also highly appreciate the freedom I received and the possibility I got to begin my PhD research project.

I am grateful to my thesis referees, jury members and all reading committees (Prof Philippe

Acknowledgements

The first author is grateful to the Erasmus Mundus LOTUS project (European Commission) and the Forum Belmont TLSCC project for providing the opportunity to do this work. This publication reflects the views only of the authors, and the European Commission cannot be held responsible for any use which may be made of the information contained therein. We would like to sincerely thank the Mekong River Commission who provided the database for our analyses. Many thanks go to L. Tudesque for his help in determining the surface area of watersheds and the percentage of land cover and K.K.Y. Wong for his editorial comments on this manuscript.

Acknowledgements

The first author is grateful to the Erasmus Mundus LOTUS Unlimited project (European Commission) and the Belmont Forum TLSCC projects for providing the opportunity to do this work. This publication reflects the views only of us the authors, and the European Commission cannot be held responsible for any use that may be made of the information contained therein. We are grateful to the Mekong River Commission who provided the database for our analyzes. Many thanks go to Dr. K.K.Y Wong for his editorial comments on this manuscript. Laboratoire Evolution & Diversité Biologique (EDB) is supported by "Investissement d'Avenir" grants (CEBA, TULIP,.

Acknowledgements

The first author would like to sincerely thank Erasmus Mundus LOTUS project (European Commission) and the Forum Belmont TLSCC project for providing the possibility to do this research. This publication reflects the views only of the authors, and the European Commission cannot be held responsible for any use which may be made of the information contained therein. We are grateful to the Mekong River Commission who provided the database for our analyses.

Acknowledgements

We are grateful to the Flemish Environment Agency (VMM) for the opportunity to study their samples. The first author would like to sincerely thank Erasmus Mundus LOTUS project (European Commission) and the Forum Belmont TLSCC project for providing the possibility to do this research. This publication reflects the views only of the authors, and the Commission cannot be held responsible for any use which may be made of the information contained therein.

ACKNOWLEDGEMENTS

We would like to thank the Flemish Environment Agency (VMM) for the opportunity to study their samples. RS is supported by the Erasmus Mundus LOTUS (European Commission) and the Forum Belmont TLSCC project. This publication reflects the views only of the authors, and the European Commission cannot be held responsible for any use which may be made of the information contained therein.

My grateful appreciation is due to lecturers and senior research staff from the Faculty of Science, Royal University of Phnom Penh (RUPP),

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.limno.2017.04.001. Species occurrence Logistic regression Random forest Artificial neural network Support vector machine Macroinvertebrates Habitat suitability Mekong river a b s t r a c t Predictive models are useful to support decision making, management and conservation planning. However, the performance of models varies across techniques and is affected by several factors including species prevalence (i.e. the occurrence rate of each species in the total samples). Here, we analysed and compared the performance of four common modelling techniques based on the species prevalence. The occurrence of macroinvertebrates collected at 63 sites along the Lower Mekong Basin was predicted using Logistic Regression, Random Forest, Support Vector Machine and Artificial Neural Network (ANN). Model performance was evaluated using Cohen's Kappa Statistic (Kappa), area under receiver operating characteristic curve (AUC) and error rate. We found a highly significant quadratic effect of species prevalence on the four modelling techniques' performance. Kappa and AUC were less depended on the species prevalence, making them a better measure. The best performance (Kappa and AUC) was reached when predicting species with an intermediate prevalence (e.g. 0.4-0.6). The four modelling techniques significantly yielded different performances (p < 0.01), of which ANN performed generally better when using the complete prevalence range (i.e. 0.0-1.0) and the lower prevalence range (i.e. <0.1). However, the four techniques similarly performed when predicting species with a higher prevalence range (i.e. ≥0.3). Our results provide useful insights into the application of modelling techniques in predicting species occurrence and how their performance varies for species with different prevalence ranges. We suggest that the selection of appropriate modelling techniques should carefully take into account the species prevalence, particularly in the case of rare and generalist species.

© 2017 Elsevier B.V. All rights reserved.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolmodel.2017. 03.006.

Supplementary material

The following supplementary material is available for this article: Appendix 1. Mean as well as minimum and maximum values of each physical-chemical variable measured at sites where the presence of each genus of alien and native molluscs was recorded. Appendix 2. Confusion matrices showing the observed and predicted classes obtained from Classification Tree models of the four periods.

This material is available as part of online article from:

http://www.aquaticinvasions.net/2017/Supplements/AI_2017_Sor_etal_Supplement.pdf Article 4

1 Title: Optimizing the reliability of classification tree models in predicting alien mollusc occurrence: a hindcasting-and forecasting-based approach Authors: Ratha Sor 1,2,3,* , correspondence, e-mail: sorsim.ratha@gmail.com Pieter Boets 2,4 , e-mail: pieter.boets@oost-vlaanderen.be Sovan Lek 1 , e-mail: sovannarath.lek@univ-tlse3.fr Peter Goethals 2 , e-mail: Peter.Goethals@UGent.be Institutional address:

AUTHORSHIP

RS, PB and PG conceived the study, RS and SL performed the modelling work, RS wrote the first draft of the manuscript, and all authors contributed significantly to revisions.

Supporting information

Fig. S1 Mean and standard error bar indicating the overall performance of CT models based on Kappa and CCI, when using the field data (Fig. 1, 2) and combined data (Fig. 3,4). Each column corresponds to each resampling approach (random and stratified).

Fig. S2

The 288 trees built for the two scenarios using the two types of dataset (field and combined data) and the two resampling approaches (random and stratified). The following slides, from slide number 1 to 10 and from 11 to 20, are the trees using the field and combined data, respectively. For the hindcasting scenario, the models were calibrated using the data period D4 and were validated using the environmental predictor of the data period D4, D3, D2 and D1. For the forecasting scenario, the models were calibrated using the data period D1 and were validated using the environmental predictor of the data period D1, D2, D3 and D4. The calibrated models' complex parameter and their cross-validation error (xerror) are given under each tree. The increased xerror (↑) is highlighted in yellow. The decreased and stable xerrors are respectively labeled with (↓)and (~). The detailed information of each slide title is provided in the following:

-Field_Hind_V-D4(3,2,1)_Random/Stratified: models were calibrated using the field data, for hindcasting scenario (validation based on D4, D3, D2, D1), and using random/ stratified resampling approach.

-Field_Fore_V-D1 [START_REF] De Troyer | Water quality assessment of streams and wetlands in a fast growing East African City[END_REF][START_REF]Summary Statistics for Globally Threatened Species. Table 1: Numbers of threatened species by major groups of organisms[END_REF]4)_Random/Stratified: models were calibrated using the filed data, for forecasting scenario (validation based on D1, D2, D3, D4), and using random/stratified resampling approach.

-Comined_Hind_V-D4(3,2,1)_Random/Stratified: models were calibrated using the combined data, for hindcasting scenario (validation based on D4, D3, D2, D1), and using random/ stratified resampling approach.

Comined_Fore_V-D1 [START_REF] De Troyer | Water quality assessment of streams and wetlands in a fast growing East African City[END_REF][START_REF]Summary Statistics for Globally Threatened Species. Table 1: Numbers of threatened species by major groups of organisms[END_REF]4)_Random/Stratified: models were calibrated using the combined data, for forecasting scenario (validation based on D1, D2, D3, D4), and using random/stratified resampling approach.

Table S1. Summary of the 288 models which were built based on the two types of dataset (field and combined data) and the two resampling approaches.

R: random resampling, S: stratified resampling, CV: cross-validation, Rep: replicate, which refers to a different seed set chosen for the determination of new folds.

Supporting Information

Figure S1 Mean and standard error bar indicating the overall performance of CT models based on Kappa and CCI, when using the field data (Fig. 1, 2) and combined data (Fig. 3,4). Each column corresponds to each resampling approach (random and stratified).

Fig. 1 The overall performance of models using the field data.

Supporting Information

Figure S2 The 288 trees built for the two scenarios using the two types of dataset (field and combined data) and the two resampling approaches (random and stratified). The following slides, from slide number 1 to 20 and from 21 to 40, are the information and trees built using the field and combined data, respectively. For the hindcasting scenario, the models were calibrated using the data period D4 and were validated using the environmental predictor of the data period D4, D3, D2 and D1. For the forecasting scenario, the models were calibrated using the data period D1 and were validated using the environmental predictor of the data period D1, D2, D3 and D4. The calibrated models' complex parameter and their crossvalidation error (xerror) are given under each tree. The increased xerror (↑) is highlighted in yellow. The decreased and stable xerrors are respectively labeled with (↓)and (~). The detailed information of each slide title is provided in the following:

-Field_Hind_V-D4(3,2,1)_Random/Stratified: models were calibrated using the field data, for hindcasting scenario (validation based on D4, D3, D2, D1), and using the random/stratified resampling approach.

-Field_Fore_V-D1 [START_REF] De Troyer | Water quality assessment of streams and wetlands in a fast growing East African City[END_REF][START_REF]Summary Statistics for Globally Threatened Species. Table 1: Numbers of threatened species by major groups of organisms[END_REF]4)_Random/Stratified: models were calibrated using the filed data, for forecasting scenario (validation based on D1, D2, D3, D4), and using the random/stratified resampling approach.

-Comined_Hind_V-D4(3,2,1)_Random/Stratified: models were calibrated using the combined data, for hindcasting scenario (validation based on D4, D3, D2, D1), and using the random/ stratified resampling approach.

-Comined_Fore_V-D1 [START_REF] De Troyer | Water quality assessment of streams and wetlands in a fast growing East African City[END_REF][START_REF]Summary Statistics for Globally Threatened Species. Table 1: Numbers of threatened species by major groups of organisms[END_REF]4)_Random/Stratified: models were calibrated using the combined data, for forecasting scenario (validation based on D1, D2, D3, D4), and using the random/stratified resampling approach.