Ariane Le Gros
email: ariane.le.gros@gmail.com

Philippe Clergeau

Dario Zuccon

Raphaël Cornette

Blake Mathys

Sarah Samadi

Invasion history and demographic processes associated with rapid morphological changes in the Red-whiskered bulbul established on tropical islands

Keywords: contemporary evolution, invasive species, neutral evolution, phylogeography, population genetics, Pycnonotus jocosus

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

Introduction

The invasion history of a species established outside of its native range and the demographic processes following introduction events can impact the genetic diversity and thus the evolutionary trajectory of the newly established populations. For example, propagule pressure (i.e. the number of individuals introduced in a location) is classically positively correlated with genetic diversity in invasive populations [START_REF] Lockwood | The role of propagule pressure in explaining species invasions[END_REF][START_REF] Blackburn | The influence of numbers on invasion success[END_REF]. Multiple introductions can also increase genetic diversity, and admixture between individuals from different source populations can generate new genotypes [START_REF] Curran | Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions[END_REF][START_REF] Kolbe | Admixture determines genetic diversity and population differentiation in the biological invasion of a lizard species[END_REF][START_REF] Roy | Hybrid "superswarm" leads to rapid divergence and establishment of populations during a biological invasion[END_REF]. Conversely, sequential introductions in new locations associated with bottlenecks and rapid range expansion can increase the impact of drift and decrease genetic diversity [START_REF] Clegg | Genetic consequences of sequential founder events by an island-colonizing bird[END_REF][START_REF] Curran | Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions[END_REF][START_REF] Peischl | Expansion load: recessive mutations and the role of standing genetic variation[END_REF]. Changes in genetic diversity will in turn affect the evolutionary potential of an invasive population and thus its ability to spread where it has been introduced [START_REF] Lee | Evolutionary genetics of invasive species[END_REF]Bacigalupe 2008). Knowing the invasion history of populations introduced in new locations (i.e. date of introduction(s), number of introductions, propagule pressure, source population(s), demographic history, admixture) is thus an important element in designing control and prevention strategies against the spread of these populations [START_REF] Estoup | Reconstructing routes of invasion using genetic data: why, how and so what?[END_REF][START_REF] Jacquet | Colonization of the Mediterranean Basin by the vector biting midge species Culicoides imicola : an old story[END_REF][START_REF] Rijal | Reconstructing the invasion history of Heracleum persicum (Apiaceae) into Europe[END_REF].

Recently, it has been suggested that rapid local adaptation is a factor enhancing the invasive success of introduced populations [START_REF] Lee | Evolutionary genetics of invasive species[END_REF][START_REF] Lee | Evolutionary origins of invasive populations[END_REF]. Many studies have identified cases of rapid phenotypic changes in introduced populations (e.g. [START_REF] Losos | Adaptive differentiation following experimental island colonization in Anolis lizards[END_REF][START_REF] Stockwell | Translocations and rapid evolutionary responses in recently established populations of western mosquitofish (Gambusia affinis)[END_REF][START_REF] Kristjánsson | Rapid divergence in a recently isolated population of threespine stickleback (Gasterosteus aculeatus L.)[END_REF][START_REF] Yonekura | A peculiar relationship between genetic diversity and adaptability in invasive exotic species: bluegill sunfish as a model species[END_REF][START_REF] Curran | Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions[END_REF][START_REF] Kooyers | Rapid evolution of an adaptive cyanogenesis cline in introduced North American white clover (Trifolium repens L.)[END_REF]Adachi et al. 2012;[START_REF] Rollins | A genetic perspective on rapid evolution in cane toads (Rhinella marina)[END_REF]. However, it is often unclear whether these phenotypic changes result from natural selection or from other mechanisms such as genetic drift or admixture [START_REF] Colautti | Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation[END_REF]. Indeed, genetic drift can lead to non-adaptive phenotypic changes. Admixture can generate new genotypic combinations and potentially new phenotypes that can differ in fitness from parental forms in the new local conditions. Therefore, in addition to predicting invasive potential, documenting the invasion history and the genetic structure of invasive populations should give insights on the roles of neutral processes in phenotypic changes that can be observed in these populations [START_REF] Gonzalez-Quevedo | Drift, not selection, shapes toll-like receptor variation among oceanic island populations[END_REF].

The Red-whiskered Bulbul (Pycnonotus jocosus, Carolus Linnaeus, 1758) is a popular cage bird native to South-eastern Asia [START_REF] Peters | Check-list of birds of the world[END_REF][START_REF] Clergeau | Fast colonization of an introduced bird: the case of Pycnonotus jocosus on the Mascarene Islands[END_REF]. It has been introduced to and successfully invaded many tropical and subtropical regions all over the world, mainly during the 1960s and 1970s [START_REF] Lever | Naturalised birds of the world (A & C Black[END_REF]: the South-Eastern coast of Australia (New South Wales and South Australia), the United States (California, Florida and Hawaii), and several islands of the Indian Ocean (Comoros, Mauritius, Reunion Island and Seychelles). Population size increased very rapidly in all these locations [START_REF] Carleton | The red-whiskered bulbul in Florida: 1960-71[END_REF][START_REF] Williams | Differential range expansion and population growth of bulbuls in Hawaii[END_REF][START_REF] Clergeau | Fast colonization of an introduced bird: the case of Pycnonotus jocosus on the Mascarene Islands[END_REF] even though the environments are quite diverse in terms of biotic and abiotic conditions, suggesting that this species is able to adapt rapidly to new ecological conditions. The Red-whiskered bulbul is considered a pest in a part of its invasive range due to the damages it causes to orchards [START_REF] Williams | Differential range expansion and population growth of bulbuls in Hawaii[END_REF][START_REF] Mandon-Dalger | Relationships between alien plants and an alien bird species on Reunion Island[END_REF][START_REF] Lever | Naturalised birds of the world (A & C Black[END_REF]) and its negative impacts on the invaded ecosystems where it can compete with, or predate on, endemic species [START_REF] Lever | Naturalised birds of the world (A & C Black[END_REF][START_REF] Linnebjerg | Diet composition of the invasive red-whiskered bulbul Pycnonotus jocosus in Mauritius[END_REF], serve as a reservoir for parasites affecting local fauna [START_REF] Shehata | Changes in the native and introduced bird populations on O'ahu, infections diseases and species replacement[END_REF], or facilitate the dispersal of seeds of invasive plants [START_REF] Carleton | The red-whiskered bulbul in Florida: 1960-71[END_REF][START_REF] Mandon-Dalger | Relationships between alien plants and an alien bird species on Reunion Island[END_REF][START_REF] Linnebjerg | Diet composition of the invasive red-whiskered bulbul Pycnonotus jocosus in Mauritius[END_REF]. Despite these impacts, the only information available on invasive populations is generally the date of first introduction or first observation.

Rapid phenotypic changes have been previously reported in an invasive population of Red-whiskered bulbul. Amiot et al. (2007) showed that in Reunion Island, individuals inhabiting opposite coasts have distinct morphologies. These two coasts present contrasting ecological conditions: the coast exposed to the prevailing wind (windward) is more humid and densely forested than the opposite coast (leeward). Historical sources suggest that the Redwhiskered bulbul was introduced only once on Reunion from Mauritius in 1972 at the south-eastern point of the island (Saint Philippe) and colonized both coasts from this point [START_REF] Clergeau | Fast colonization of an introduced bird: the case of Pycnonotus jocosus on the Mascarene Islands[END_REF][START_REF] Lever | Naturalised birds of the world (A & C Black[END_REF]. Considering the reproductive success and rapid spread of this species on Reunion, Amiot et al. (2007) hypothesized that rapid local adaptation could explain the morphological differences observed between the windward and leeward coasts, and thus the invasive success of the Red-whiskered bulbul on Reunion. In addition, the two coasts are separated by a high mountain range which is thought to constitute an obstacle to the dispersal of Red-whiskered bulbuls, isolating the two populations, and might thus facilitate local adaptation (Amiot et al. 2007). However, some Red-whiskered bulbuls were observed at the opposite position from the known introduction point (Sainte Marie) in 1978, long before the population introduced in Saint Philippe could have reached this location by natural colonization as it only started to spread in the surroundings of St Philippe in 1980 [START_REF] Clergeau | Fast colonization of an introduced bird: the case of Pycnonotus jocosus on the Mascarene Islands[END_REF]. This suggests that bulbuls were transported by humans across Reunion or that unknown introduction events took place. With such an uncertain historical context it is not possible to exclude the possibility that the phenotypic differences described by Amiot et al. (2007) result from the invasion history of the population rather than from natural selection.

In this article, we wanted to test the hypothesis that the morphological differences observed between the two coasts of Reunion are due to neutral evolutionary processes. We used three approaches to test this hypothesis. First, we studied the morphology of Red-whiskered bulbuls established in two other tropical islands very similar to Reunion in terms of size, geography and climate: Mauritius and Oahu (Hawaii). Like Reunion, these two islands are divided by a mountain range and have a windward coast and a leeward coast. Redwhiskered bulbuls are common birds in these three islands. Historical sources suggest that the Red-whiskered bulbul was introduced on Mauritius, probably from India, in 1892, whereas Oahu was colonized approximately at the same period as Reunion (1965) but from an unknown source [START_REF] Lever | Naturalised birds of the world (A & C Black[END_REF]. We predicted that if the morphological differences observed between leeward and windward coasts of Reunion are adaptive, we would find the same morphological pattern on Mauritius and Oahu. On the contrary, if these morphological differences are due to neutral processes, we expected to find different morphological patterns on each island. Secondly, we investigated with molecular data the invasion history of Red-whiskered bulbuls established on these three islands. In particular, we used an approximate Bayesian computation approach widely used to reconstruct the invasion history of invasive species (e.g. [START_REF] Jackson | Ancestral origins and invasion pathways in a globally invasive bird correlate with climate and influences from bird trade[END_REF][START_REF] Jacquet | Colonization of the Mediterranean Basin by the vector biting midge species Culicoides imicola : an old story[END_REF][START_REF] Rijal | Reconstructing the invasion history of Heracleum persicum (Apiaceae) into Europe[END_REF][START_REF] Wei | Population genetic structure and approximate Bayesian computation analyses reveal the southern origin and northward dispersal of the oriental fruit moth Grapholita molesta (Lepidoptera: Tortricidae) in its native range[END_REF]. We predicted that if multiple introductions from different origins, founder effects or genetic drift are the cause of the morphological patterns observed, we would find evidence of multiple introductions and a neutral genetic structure matching the morphological patterns observed. Finally, we compared the morphology of Redwhiskered bulbuls between the invasive populations and the source populations identified with molecular data. We predicted that if the morphological patterns found on the islands are due to neutral processes, each population would have diverged differently from its source population whereas if the morphological patterns are the result of adaptation to local conditions, all the populations experiencing similar conditions would have diverged towards a similar morphology.

Material and methods

Samples

Tissue samples and morphological data were collected from the whole native range of the Red-whiskered bulbul and from three invasive populations (Reunion, Mauritius and Oahu). Sampling live individuals from the whole native range was unrealistic considering it covers a very large region in Asia. Instead toe-pads from 27 museum specimens covering evenly the whole native range and all recognized subspecies were obtained from the collections of several museums (figure 1b and table S1). Morphological measurements and pictures were taken by the same person (ALG) on 191 adult specimens from the native range in the collection of the British Natural History Museum. For invasive populations, feathers and morphological data were collected on wild-caught individuals. Morphological measurements and pictures were taken by the same person (ALG) and only on adult individuals. In 2002 and 2003, about 400 Red-whiskered bulbuls were captured on Reunion Island during a control program andND2). Support values indicated at the nodes are the posterior probability (threshold 0.50) and the bootstrap support (threshold 50%) obtained from the maximum likelihood analysis, and the Bayesian approach respectively. b) Geographic distribution of the nine subspecies of Red-whiskered Bulbul in its native range as described in [START_REF] Peters | Check-list of birds of the world[END_REF]. Symbols show the geographic origin of the specimens used in the phylogeographic analysis. The color and form of the symbols indicate the phylogenetic clade to which they belong.

organized by FDGDON (Fédération Départementale des Groupements de Défense contre les Organismes Nuisibles de la Réunion). They were captured in 12 sites spread along the windward and the leeward coasts (figure 3b and table S2). In 2013, mist-nets were used to capture 50 adult Redwhiskered bulbuls on Mauritius and 45 on Oahu. On Mauritius, bulbuls were caught in three biological stations administrated by an NGO, the Mauritian Wildlife Foundation. One site was located on the windward coast (M-IA), one site was on the leeward side of the island (M-B) and the last one was in the mountain range delimiting the windward and leeward coasts but on the leeward side (M-C; figure 3b and table S2). On Oahu, the bulbuls were caught in two sites, one on each coast (figure 3b and table S2). Wild-caught individuals were sexed using the PCR-based protocol of [START_REF] Griffiths | A DNA test to sex most birds[END_REF], whereas for museum specimens we relied on the information available on specimen labels.

DNA extraction, amplification and sequencing

In order to prevent contamination from modern samples, DNA extraction and PCR amplification from museum samples were conducted in a separate laboratory dedicated to ancient DNA work. In addition, equipment and benches were decontaminated before and after use by application of UV light and bleach. Finally, negative controls were added during amplifications. Total genomic DNA was extracted from toe-pads or feathers with the QIAamp DNA Micro Kit (Qiagen) following the manufacturer instructions for the blood and tissue samples. The digestion volume was doubled, with the final concentration of 2mg/mL for Proteinase K and 50 mM for dithiothreitol. For the phylogeographic analysis, a region was amplified in the two mitochondrial genes COI and ND2 (655 bp and 564 bp respectively). For the fresh tissue samples, the two genes were amplified in one fragment whereas in the case of museum samples, short overlapping fragments (ca. 200 bp) were amplified with internal primers (table S3). The amplification protocols are described in Appendix S1 of the Supporting Information. Positive PCR products were sequenced in both directions on an Applied Biosystems 3130XL DNA sequencer. The sequences were assembled in CODONCODE 4.0 (CODONCODE Corporation). For the population genetics study, ten polymorphic microsatellite loci were selected from the literature (table S4) and amplified following the protocol described in Appendix S2 of the Supporting Information. Samples were genotyped on an Applied Biosystems 3130XL DNA sequencer. Genotypes were scored with GeneMapper 4.0 (Applied Biosystems) and checked manually.

Phylogeographic analysis

In addition to the 27 museum specimens, the COI and ND2 sequences of three additional specimens from the native range were retrieved from Genbank (table S1). The COI and ND2 genes were also sequenced in a total of 34 individuals caught on both coasts in the three invaded islands: Reunion (8 individuals), Mauritius (8 individuals), and Oahu (18 individuals, table S1). Seven individuals from an invasive population located around Sydney (Australia) were also added to see whether there is only one source for all invasive populations worldwide (table S1). Finally, following the results of [START_REF] Moyle | Phylogenetic relationships of the bulbuls (Aves: Pycnonotidae) based on mitochondrial and nuclear DNA sequence data[END_REF], Pycnonotus sinensis, Pycnonotus barbatus and Pycnonotus cafer were selected as out-groups, and sequences for one individual of each species were retrieved from Genbank (table S1). The COI and ND2 genes were concatenated in a partitioned dataset analyzed under the Bayesian inference and the maximum likelihood criteria. These two procedures are described in Appendix S3 of the Supporting Information.

Morphological data

Traditional morphological measurements and geometric morphometric data were used to study the morphology of 191 individuals from the native range and more than 300 individuals from invasive populations (Reunion: 231; Mauritius: 50; Oahu: 45). Five traditional measurements were recorded: beak length, beak width, beak depth, folded wing length and central tail feathers length.

Log-shape ratios [START_REF] Mosimann | New statistical methods for allometry with application to Florida red-winged blackbirds[END_REF] were used in order to allow the study of morphological variables independent of size. Geometric morphometrics was used to describe more precisely the conformation of the beak [START_REF] Zelditch | Geometric morphometrics for biologists: a primer[END_REF]. Pictures in lateral view of the beak of each individual were taken in standardized conditions. TPSDIG 2 [START_REF] Rohlf | TPSDIG[END_REF]) was used to digitize four landmarks (anatomical points) and 20 sliding semi-landmarks (non-homologous points) from these pictures (figure S1, [START_REF] Gunz | Semilandmarks: a method for quantifying curves and surfaces[END_REF]). All pictures were digitized by the same person (ALG) and the repeatability of the digitization process was tested using a principal component analysis (PCA) on three repetitions taken on five specimens chosen randomly in the same subspecies (figure S2).

Morphometric analyses

Statistical analyses were done with R 2.15.3 (R Core Team 2013) and using the libraries ade4 [START_REF] Thioulouse | ADE-4: a multivariate analysis and graphical display software[END_REF], Ape [START_REF] Paradis | APE: Analyses of Phylogenetics and Evolution in R language[END_REF], Hotelling [START_REF] Curran | Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions[END_REF] and Rmorph (Baylac 2012). All the analyses were done separately for the traditional measurements taken on the whole body, and for geometric morphometric data on the beak. In order to assess if there are morphological differences between individuals with different origins, we defined two levels of geographic origin: the global geographic origin (the four regions defined by phylogeographic clade in the native range, Reunion, Mauritius, and Oahu), and nested into this variable, the sample sites where the individuals were captured. Size-Analyses of covariance (ANCOVAs) were performed on body size and on beak size with global geographic origin and sample site as explanatory variables. Sex and the interactions between geographic groups and sex were added as co-factors. Tuckey's HSD post hoc tests were then used to test for pairwise differences between global geographic origins, sample sites and coasts within each island. Conformation-PCAs were performed on the two data sets (whole body and beak) and the principal component axes representing 95% of the total variance were kept as conformation variables for the following analyses. Before assessing an effect of origin on the conformation of individuals, we checked for differences in allometry between geographic groups using ANCOVAs performed for each conformation variable and for each data set with size, global geographic origin, sample site, sex and all possible interactions between them as explanatory factors. If some interactions between size and geographic groups were significant, the conformation variables were corrected for size by regressing size on it. Then, the effect of geographic groups on the conformation of individuals was assessed with multivariate analyses of covariance (MANCOVAs) performed for each data set, with conformation variables as response variables, and global geographic origin and sample site as explanatory variables. Sex and the interactions between geographic groups and sex were added as co-factors. Hotelling's T-squared post hoc tests were used to test for pairwise differences in conformation between global geographic origins, sample sites and coasts within each island. The threshold of acceptance of the null hypothesis was adjusted following the Bonferroni correction for global geographic origins and coasts, and following the Holm correction for sample sites (as the number of pairwise tests was elevated in the last case).

Neutral genetic structure

The ten microsatellite loci were amplified in a total of 480 individuals (Reunion: 385, Mauritius: 50 and Oahu 45). The presence of null alleles was assessed with FREENA [START_REF] Chapuis | Microsatellite null alleles and estimation of population differentiation[END_REF]). Sample sites with fewer than 20 individuals were excluded from the analysis to increase its detection power. Mean number of alleles, Shannon's information index, observed heterozygosity, unbiased expected heterozygosity, and fixation index were assessed over all loci and for each sample site with GENEALEX 6.5 [START_REF] Peakall | GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update[END_REF]. Deviation from Hardy-Weinberg equilibrium and linkage disequilibrium between pairs of loci were tested for each sample site with GENEPOP 4.2.1 [START_REF] Rousset | Genepop'007: a complete re-implementation of the genepop software for Windows and Linux[END_REF] using default parameter values.

Pairwise FST values were calculated between all sample sites with FREENA [START_REF] Chapuis | Microsatellite null alleles and estimation of population differentiation[END_REF]. The Bayesian clustering approach implemented in STRUCTURE 2.3.3 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies[END_REF] was used to describe the genetic structure in the data set. Ten runs were performed for each value of K (number of clusters) from 1 to 17 (burn-in period: 1.10 5 , 1.10 6 iterations). The admixture model and the assumption of correlated allele frequencies were chosen.

Isolation by distance and migration on Reunion Island

In order to assess the causes of the genetic structure found on Reunion Island, isolation by distance (IBD) was investigated. The linearized genetic differentiation between each sampling site: Fst/(1-Fst), was used as genetic distance. Fst values were computed with FREENA (1.10 4 permutations). Two kinds of geographic distances were used. First, the mountain range was not considered as an obstacle and the Euclidian distance between sites was used. Secondly, the mountain range was considered as an obstacle and the distances between sites were calculated following the coastline. In the first case, the distances were logtransformed as movements in two dimensions were considered. In the second case, the raw distances were used directly as only movements on one direction were considered [START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF]. IBD between sample sites was tested first on the whole island, and then within the groups delimited by STRUCTURE. This approach allows to account for the effect of "by chance" geographic separation of groups that are also genetically differentiated [START_REF] Perez | Effects of landscape features and demographic history on the genetic structure of Testudo marginata populations in the southern Peloponnese and Sardinia[END_REF]. The P-values of the correlation coefficients were computed using a Mantel test (1.10 5 permutations).

In addition, actual migration rates between the groups delimited by STRUCTURE were assessed with BAYESASS 1.3 [START_REF] Wilson | Bayesian inference of recent migration rates using multilocus genotypes[END_REF]. The parameters were chosen to fit with the author's recommendations: mixing parameters for allele frequencies: 0.12; inbreeding coefficients: 0.15 and migration rates: 0.07; number of iterations = 1.10 7 ; burnin = 1.10 6 ; sampling frequency: 100. The sample site R-W2 was not used in this analysis as its individuals could not be confidently assigned to one of the two genetic clusters defined by STRUCTURE.

Invasion scenarios

The approximate Bayesian computation (ABC) method implemented in DIYABC 2.0.3 [START_REF] Cornuet | Inferring population history with DIY ABC: a user-friendly approach to approximate Bayesian computation[END_REF]) was used on the microsatellite dataset to investigate the invasion history of the Red-whiskered bulbul in Reunion. This method can be used to compare potential invasion scenarios proposed by the user [START_REF] Estoup | Reconstructing routes of invasion using genetic data: why, how and so what?[END_REF]. We tested whether the two genetic clusters identified by STRUCTURE on Reunion (RB and RO) were founded by a single introduction event or independent ones. Four scenarios were compared. In the first two scenarios, a single introduction event happened on Reunion and the population then split in two (figure 4a, scenarios A and B). These two scenarios simulate the hypotheses of IBD or individuals transported by human within the island to explain the formation of the two clusters. In the two other scenarios, the two clusters were founded by independent introduction events from Mauritius (figure 4a, scenarios C and D). For each pair of scenarios, both possible orders for the foundation of the clusters were tested. Mauritius (M) was used as the source of all the introduction events on Reunion as historical data suggested that bulbuls were introduced to Reunion from there and as our molecular data (mitochondrial sequences and microsatellites) support this hypothesis. The priors for the demographic parameters were chosen according to historical and biological data available with ranges large enough to allow flexibility in the model (prior set 1, table S5). The microsatellite loci were separated in three groups according to their repeat motif. The mutation model parameters were left to default values following [START_REF] Gotzek | Global invasion history of the tropical fire ant: a stowaway on the first global trade routes[END_REF]. Within sample summary statistics used were means of number of alleles, genetic diversity, size variance and Garza-Williamson's M. Among samples summary statistics used were the means of the number of alleles, genetic diversity, size variance Fst values, shared allele distances and (dµ)² distances.

For each scenario, one million data sets were simulated, and the posterior probability of each scenario was assessed using a polychotomous logistic regression on the 1% simulated data sets closest to the observed data set. The posterior distributions of demographic parameters were estimated under the most probable invasion scenario using a local linear regression on the 1% closest data sets. In order to evaluate the robustness of our inferences a second analysis was run with an alternative set of priors (prior set 2, table S5) and the posterior probabilities of all scenarios were estimated using the 0.1% and 1% closest simulated data sets for both sets of priors. The ability of the program to correctly choose the true scenario was evaluated by analyzing new data sets simulated from the different scenarios. For each scenario, 100 data sets were simulated using parameter values drawn from the prior set 1. Posterior probabilities of scenarios were estimated for each simulated test data set using the 1% closest data sets. These probabilities were used to compute type I and II errors in scenario choice. The data from the sample site R-W2 were not included in this analysis as the individuals could not be confidently assigned to one of the two genetic clusters defined by STRUCTURE.

Finally, the program BOTTLENECK 1.2 [START_REF] Piry | BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data[END_REF]) was used to look for signatures of recent bottlenecks in the invasive populations of Oahu, Mauritius and Reunion (RO and RB). The two-phase model (TPM) was used with and a proportion of single-step mutations set to 70 % with a variance of 30 among multiple steps. One tail Wilcoxon's signed rank tests were used to assess statistical significance of heterozygosity excess in each population.

Results

Phylogeographic analysis

The output of MrModelTest give the HKY+Γ and the GTR+I models as the best fit for the COI and ND2 genes, respectively. The phylogenetic analysis recovers four main clades (figure 1a). The four clades are coherent with the geographic origin of the specimens but match rather poorly with the current subspecies subdivision (figure 1b). The first clade contains the individuals from Western and Southern India (subspecies P. j. abuensis and fuscicaudatus). The second clade is composed of individuals from the lowlands of Eastern India to Western Burma (some but not all P. j. emeria, monticola and pyrrhotis) and Andaman Islands (P. j. whistleri). The third clade is composed of individuals from the Himalayan region (some but not all P. j. monticola and pyrrhotis). Finally, the individuals from the Indochinese peninsula (South-eastern Burma, China, Laos, Malaysia, Thailand and Vietnam) fall together in a fourth group (P. j. jocosus, hainanensis, pattani and some but not all emeria and monticola). All the individuals from Reunion and Mauritius share the same haplotype and fall in the Eastern Indian clade. The majority of Oahu individuals also share one haplotype and fall in the Eastern Indian clade. However, two individuals (out of 18) fall into the Indochinese clade. Finally, the individuals from Australia all belong to the Indochinese clade (figure 1a).

Morphometric analyses

Size-Global geographic origin and sample site have a significant effect on size in the ANCOVAs performed on both data sets (table 1). Sex also has a significant effect on size (females are smaller and have a smaller beak than males, figures S3a and S4a, table 1). For both data sets, post hoc tests show that 1) in most cases individuals from invasive populations are significantly different from the individuals of their source clade; and 2) there are significant differences in size between the three invasive populations (table 2, figures S3b and S4b). There are also significant differences in size between some sample sites but not all in the whole body data set (figure S3c). In the beak data set, there are no significant differences in size between any pair of sample sites within each island (table S6, figure S4c).

When we compare the size of the individuals between the windward and leeward coasts, there are significant differences in body size between the coasts in the three islands and no differences in beak size (table 3). Conformation-There are significant differences in allometry between geographic groups in the whole body data set (table S7). The morphological variables were therefore corrected for size in this data set. In the whole body data set, there is a significant effect of global geographic origin and sample site on conformation (table 1, figure 2a andc). Sex also has a significant effect on conformation but not the interactions between geographic groups and sex (table 1). In the beak data set, there is a significant effect of the triple interaction between sex, global geographic origin and sample site (P=1.95.10 -2 ), therefore the MANCOVA x: Non significant differences. females for beak conformation; table 2, figure 2a andb). Within islands, there are differences between some sample sites, but not all, in both data sets (figure 2c, d ande, table S8). When we compare the conformation of individuals between the windward and leeward coasts on each island, there are significant differences in conformation in Reunion and Mauritius in the whole body data set and also in Reunion in the beak data set but only for males (table 3). On Oahu, the two coasts are not significantly different in both data sets.

Neutral genetic structure Amplification of the microsatellite loci was successful with only 1.8% of data missing over all loci and individuals. Null alleles were detected at locus TG05-046 with FREENA. The data from this locus were therefore not used. None of the sample sites significantly deviate from Hardy-Weinberg equilibrium nor present linkage disequilibrium (table 4). As expected in invasive populations the genetic diversity is not very high: the mean number of alleles per locus is usually between two and three and the observed and expected heterozygosity are under 0.50 (table 4). The genetic diversity is a little higher in Mauritius than in the two other islands (table 4). Genetic differentiation is low between the two sites in Oahu (FST=0.04) and the three sites in Mauritius (FST ≤0.01). On Reunion Island, two groups of sites with low levels of differentiation (FST ≤0.02) can be distinguished: the sites R-W5, R-W4, R-W3 on one side and all the other sites in the other group (except the site R-W2 which is little differentiated compared to all the other sites on Reunion, table 5). When the sites of different islands are compared, Mauritius is less differentiated from Reunion than from Oahu (table 5).

The log likelihood of the simulations run with STRUCTURE increases sharply until K=4, and starts decreasing afterwards, whereas the delta K presents a first peak at K=2 and a second one at K=4 (figure S5). For K=2, the individuals from Reunion are separated into two clusters: the sites R-W3, R-W4 and R-W5 (orange cluster, figure 3a) and the other sites (blue cluster, figure 3a). These clusters do not match with the windward and leeward coasts. Individuals from Oahu fall in the orange cluster whereas each individual from Mauritius is assigned to both clusters in approximately equal proportions. For K=4, all individuals from Oahu are grouped in a single cluster as well as individuals from Mauritius. The individuals from Reunion are separated into the same two clusters obtained for K=2 (figure 3a andb).

Isolation by distance and migration on Reunion Island

When all sampled sites are considered, the correlation between genetic and geographic distances, whatever the type of distances, is x: Non significant differences. N: Mean number of individuals, Na: number of alleles, I: Shannon's diversity Index, Ho: observed heterozygosity, uHe: unbiased expected heterozygosity, and F: fixation index per sampling site and over all loci. HWE: proportion of loci deviating from Hardy-Weinberg equilibrium and LD: proportion of pair of loci showing linkage disequilibrium per sampling site.

Table 5 Pairwise FST based on microsatellite loci between sample sites in invasive populations. Euclidian distances P=0.31, coastline distances P=0.47, figure S6). Actual migration rates of 4% per generation from the orange to the blue cluster and 5% in the other way were estimated with BAYESASS.

Site

O-K O-PC M-B M-C M-IA R-W5 R-W4 R-W3 R-W2 R-W1 R-Intro R-L1 R-L2 R-L3 R-L4 R-L5 O-PC 0.

Comparison of invasion scenarios

The scenario with two independent introductions on Reunion (scenario C) was found to have the highest posterior probability. In this scenario, a first introduction formed the blue cluster and a second introduction later formed the orange cluster (figure 4b, table S9)). The 95% confidence interval shows that scenario C can be confidently discriminated from the three other ones (figure 4b). Similar results were obtained with the second set of priors (table S9). The probability to reject a true scenario is quite high (Type I error for scenario A: 0.52, B: 0.54, C: 0.46 and D: 0.36) but this is mainly caused by the fact that scenarios A and B and scenarios C and D are too similar to be well discriminated.

The type I error drops when pair of scenarios are considered together (Type I error for scenario A & B: 0.26 and C & D: 0.15). The probability to accept a false scenario is around 0.15 (Type II error for scenario A: 0.14, B: 0.14, C: 0.17 and D: 0.18). The estimates of the demographic parameters for scenario C are coherent with historical data and support the hypothesis that both populations in Reunion experienced bottlenecks following their introduction.

The bottleneck length was estimated to be 23 generations for the orange cluster and 31 for the blue cluster, and the effective population sizes during these bottlenecks were estimated to be 30 times lower than actual size for both clusters, table S10). The analysis of heterozygosity excess conducted with BOTTLENECK also detects the signature of a bottleneck in the blue population of Reunion (P=2.44.10 -2 ) and in the populations of Oahu (P=3.91.10 -3 ) and Mauritius (P=1.86.10 -2 ), but not in the orange population of Reunion (P=1.02.10 -1 ).

Discussion

Repeatability of morphological differences between island coasts

Our morphological study shows significant differences in body and beak morphology between the two coasts of Reunion Island and thus confirms the results of Amiot et al. (2007). On average, individuals from the windward coast have shorter wings and tails than individuals of the same size from the leeward coast (figure S7). Wing and tail length can affect maneuverability and longdistance flight capacities. The differences we observed could thus reflect local adaptations to the different kind of vegetation present on each coast. These individuals also have a different beak conformation with wider, deeper and shorter beaks (figure S7), which could potentially result from differences in food items consumed [START_REF] Grant | The adaptive significance of some size trends in island birds[END_REF][START_REF] Herrel | Bite performance and morphology in a population of Darwin's finches: implications for the evolution of beak shape[END_REF].

A stable isotope analysis conducted on the same individuals did not reveal clear diet differences between the two coasts of Reunion [START_REF] Roussel | Stable isotope analysis to examine possible dietmorphology relationships among Red-whiskered Bulbuls introduced on Réunion Island[END_REF]. However, that kind of study would not detect differences between diets based on hard or soft seeds which could impact beak morphology. The individuals from the windward coast also have a larger body size than individuals from the leeward coast. Finally, more subtle differences were also detected in beak conformation with the geometric morphometric approach. Males and females are divided into several groups that do not exactly correspond to the two coasts suggesting that other factors than the environmental differences between island sides are affecting the conformation of the beak.

A morphological differentiation between the two coasts was also detected on Mauritius although only in the whole body dataset. Individuals from the windward side of Mauritius have a conformation close to the one found on the windward coast of Reunion, and those of the leeward side also have longer wings and tails and wider, deeper and shorter beak than in the windward side but this difference is more pronounced than in Reunion (figure S7). In addition, as in Reunion, the individuals from the windward side are larger than those of the leeward side. On Oahu, the individuals from the windward side are also larger than those from the leeward side. However, we did not detect significant differences in conformation between the individuals from the two coasts on this island.

In the three islands, we thus found that individuals on the windward side have a larger body size than those of the leeward side. This suggests that the difference in size between the two environments is adaptive, it could for example reflect an adaptation to different thermoregulation needs (Bergmann's rule). However, it could also result from plastic responses, for example a richer diet on windward coasts which are more humid and where fruits and insects could thus be more abundant. Concerning the conformation of the individuals, the fact that we find the same morphological pattern on Reunion and Mauritius indicates that the differences in body conformation we observed could be adaptive and driven by similar selective forces. The fact that we did not find this pattern on Oahu could be due to less marked differences in ecological conditions between the two sides of Oahu than in Reunion and Mauritius. Alternatively, the population on Oahu might lack the evolutionary potential to adapt to the different environments. Indeed, we found that the genetic diversity is lower on Oahu than on the two other islands. In conclusion we found a repeatable difference in body size between windward and leeward coasts which suggests that body size evolution on each type of coast is adaptive and driven by selective forces that are common to the three islands. The differences in conformation between coasts were also similar on Reunion and Mauritius suggesting again an adaptive evolution of morphology associated to the type of coast although this was not observed on Oahu.

Repeatability of morphological divergence from source populations

Based on mitochondrial sequences, the Redwhiskered bulbul can be divided into four main clades corresponding to four distinct geographic areas of its native range. Our phylogenetic analysis shows that the population introduced on Mauritius is native to a region embracing Eastern India and Western Burma. Our analysis also supports the historical data reporting that the populations of Reunion come from Mauritius as the individuals from both islands have similar haplotypes and are also close at neutral loci. On Oahu two distinct haplotypes from different geographic clades coexist (Eastern India/Western Burma and Indochinese peninsula) suggesting multiple introduction events, but it is possible that all individuals introduced on Oahu come from a region where both genetic clades co-exist (e.g. Myanmar). If several introductions took place, the Indochinese origin seems less important as this haplotype was found only in two individuals out of 18. We did not find evidence for a double origin in our microsatellite data set. However, admixture between individuals from different origins is not always detected in both mitochondrial and nuclear sequences [START_REF] Bradbury | Transatlantic secondary contact in Atlantic Salmon, comparing microsatellites, a single nucleotide polymorphism array and restriction-site associated DNA sequencing for the resolution of complex spatial structure[END_REF][START_REF] Patten | Biased introgression of mitochondrial and nuclear genes: a comparison of diploid and haplodiploid systems[END_REF]. In conclusion we can say that the population of Oahu originates from the same broad geographic region as the population of Mauritius but might also have in lesser proportions some origins in the Indochinese peninsula.

The sizes of the individuals and their beak have diverged between invasive populations and their source in the native range but this divergence is not the same for the three invasive populations. For body size, individuals from Reunion are not significantly different from those from Eastern India (and smaller than those from Mauritius); the individuals from Mauritius are larger than those from Eastern India; and the individuals from Oahu are larger than those from Eastern India but not significantly different from those from Indochina. Concerning beak size, the individuals from Reunion have a smaller beak than those from Eastern India (and Mauritius); and the individuals from Mauritius and Oahu have a bigger beak than those from Eastern India and Indochina. This inconsistency in body size and beak size divergence suggests that either the divergence in size is not due to natural selection or that the selective forces acting on body and beak size are different on each island although their geography, climate and vegetation seem very similar.

Conversely, the body and beak conformations of the individuals from all the invasive populations are significantly different from the conformation of individuals from the native range (except for beak conformation on Oahu). In addition, the three invasive populations have diverged in a similar direction from the native range: toward shorter wings and tails and shorter, wider and deeper beaks (figure 2). This suggests that common selective forces are acting on the phenotypes of the three invasive populations and could indicate adaptation to insular ecological conditions [START_REF] Grant | The adaptive significance of some size trends in island birds[END_REF]. This adaptation can result from natural selection and/or phenotypic plasticity. Complementary approaches such as common garden experiments, reciprocal transplants or the study of evolution in genes associated with morphology would give further insight into the role of natural selection and phenotypic plasticity in the morphological changes we observed, and whether they are adaptive [START_REF] Merilä | Climate change, adaptation, and phenotypic plasticity: the problem and the evidence[END_REF].

Neutral genetic structure, invasion history and morphological changes

On Mauritius and Oahu, we did not detect a neutral genetic structure. The differences we observed in body size and conformation between the two coasts of Mauritius and in body size between the two coasts of Oahu cannot therefore be explained by random processes and neutral evolution, otherwise one would expect to see differences in neutral loci between groups that are morphologically distinct. In the case of Mauritius and Oahu, morphological differentiations between coasts are therefore due either to natural selection and/or phenotypic plasticity.

On Reunion, there is a neutral genetic structure with two genetic clusters: one located on the Northern half of the windward coast and the other occupying the rest of the island. Isolation by distance is unlikely to explain this structure as we did not detect it in any cluster even in the blue cluster which embraces the largest part of the island, and as we found evidence for migration between the two clusters. Conversely, the comparison of invasion scenarios supports the hypothesis that the two clusters result from independent introduction events. Signatures of bottlenecks were detected both by the ABC approach and the analysis of heterozygosity excess in the blue cluster and by the ABC approach in the orange cluster. The genetic structure on Reunion is therefore probably due to founder effects rather than in situ isolation and subsequent drift. These founder effects could also explain the differences in morphology found between the two coasts of Reunion. However, the neutral genetic structure does not exactly match with the morphological pattern which is more complex. When body conformation is compared between sites, the sites of each coast are more or less similar between themselves and different from the sites of the other coast (except for R-L5, figure 2c). This is also true for body size (figure S3). For beak conformation, two groups of sites are detected in both males and females but they do not match the coasts or the neutral genetic clusters (figure 2d ande). It is thus possible that several combined mechanisms (i.e. neutral processes, natural selection and/or phenotypic plasticity) are responsible for the morphological patterns observed on Reunion.

Finally, concerning the morphological divergence between invasive populations and their sources, the analysis of the microsatellite data shows that genetic diversity is not very high in the populations of the three islands and that they are genetically distinguishable even if they come from the same region of the native range. In addition, a striking lack of genetic diversity is also observed at the mitochondrial loci (one shared haplotype on Reunion and Mauritius, and two haplotypes on Oahu). Finally signatures of bottlenecks were detected in each invasive population (although not by the heterozygosity excess analysis for the orange cluster on Reunion). This suggests impacts of founder effects and bottlenecks on the genetic diversity of the invasive populations which could explain a part of the morphological divergence we observed between the invasive populations and their sources. In conclusion, neutral evolution can be confidently excluded as a factor explaining the morphological differentiations observed between the windward and leeward coast on Mauritius and Oahu. Conversely, it is possible that it explains, at least in part, the morphological differentiation observed between coast on Reunion, and the morphological divergence observed between the invasive populations and the native range.

Conclusion

Phenotypic changes in natural populations, and especially in invasive populations, can stem from different causes such as founder effects, drift, admixture, natural selection even on very short time scales, and phenotypic plasticity. In this article we show that it is necessary to study the historical context in which such phenotypic changes happen to understand the mechanisms involved and especially the role of neutral processes in these changes. Using this approach in the case of the invasive Redwhiskered bulbul, we were able to rule out the hypothesis of neutral evolution to explain phenotypic differences between windward and leeward coasts on Mauritius and Oahu whereas we showed that neutral evolution can explain, at least in part, phenotypic differences between coasts on Reunion and morphological divergence between invasive populations and their sources. Our comparative morphological study was complementary to this approach as it showed that body size and conformation differences between coasts on the three islands as well as divergence in conformation between invasive populations and their native ranges were repeated and thus likely due to adaptive processes. By combining these two approaches we can conclude that both neutral and rapid adaptive processes shaped the morphological divergence between invasive populations and their sources as well as the morphological differences between coasts on Reunion whereas rapid adaptive processes alone are responsible for the morphological differences between coasts on Mauritius and Oahu. Precise knowledge of the causes of phenotypic changes can be useful in the elaboration of control strategies for invasive species introduced in new environments or of conservation strategies for endangered species facing changes in their environment [START_REF] Santamaría | Evolution in biodiversity policy -current gaps and future needs[END_REF]. On the other hand, incomplete knowledge can lead to inadequate strategies and have dramatic consequences in both cases. This is why we argue that the role of neutral evolutionary processes should always be assessed alongside other evolutionary forces.

body R-W5 R-W4 R-W3 R-W2 R-W1 R-Intro R-L1 R-L2 R-L3 R-L4 R-L5 R-W4 x R-W3 2.98e-02 x R-W2 x x x R-W1 x x x x R-Intro x x 2.69e-03 x x R-L1 x x 4.77e-02 x x x R-L2 x x 5.94e-03 x x x x R-L3 x x 9.90e-03 x x x x x R-L4 5.54e-03 2.69e-02 3.14e-02 x x x x x x R-L5 x x <1.00e-05 x x x x x x x R-L6 x x 1.01e-02 x x x x x x x x Beak R-W5 R-W4 R-W3 R-W2 R-W1 R-Intro R-L1 R-L2 R-L3 R-L4 R-L5 R-W4 x R-W3 x x R-W2 x x x R-W1 x x x x R-Intro x x x x x R-L1 x x x x x x R-L2 x x x x x x x R-L3 x x x x x x x x R-L4 x x x x x x x x x R-L5 x x x x x x x x x x R-L6 x x x x x x x x x x x
x Non-significant differences. Nx: number of individuals in the population x (constant in time), Nxb: number of individuals in the population x during the bottleneck following introduction (constant during the whole bottleneck), ti: number of generations between present and an introduction event or a split i, tbi: number of generations during the bottleneck following an introduction event i.

Fig. 1

 1 Fig. 1 Phylogeographic analysis based on the mitochondrial genes COI and ND2 sequenced in a sub-sample of individuals from the native range and invasive populations. a) Majority rule consensus tree obtained from the mixed-model Bayesian analysis of the concatenated dataset (COIand ND2). Support values indicated at the nodes are the posterior probability (threshold 0.50) and the bootstrap support (threshold 50%) obtained from the maximum likelihood analysis, and the Bayesian approach respectively. b) Geographic distribution of the nine subspecies of Red-whiskered Bulbul in its native range as described in[START_REF] Peters | Check-list of birds of the world[END_REF]. Symbols show the geographic origin of the specimens used in the phylogeographic analysis. The color and form of the symbols indicate the phylogenetic clade to which they belong.

Fig. 2

 2 Fig. 2 Morphometric analyses based on traditional measurements and geometric morphometric data collected on the individuals from the native range and Mauritius, Reunion and Oahu. a) and b) Morphospaces defined by the first two axes of the PCAs on conformation variables. The correlation between the variables and the axes are represented for the PCA performed on the whole body data set (a). The extreme conformations associated with the axes are represented for the PCA performed on the beak data set (b, red conformation: positive end of the axis; blue conformation: negative end of the axis). c), d) and e): neighborjoining trees calculated with the Euclidian distance between the centroid of each sample site. Colors refer to the geographic location of individuals and sample sites.

Fig. 3

 3 Fig. 3 Neutral genetic structure in the invasive populations based on the STRCUTURE analysis of the microsatellite dataset. a) Cluster assignments of each individual obtained with STRUCTURE for a number of cluster (K) equal to 2, 3 and 4. For each value of K the results of the 10 runs were pooled together using CLUMPP (Jakobsson & Rosenberg 2007). Each vertical line represents a single individual and individuals are grouped by sample sites and island. b) Average cluster assignments of individuals for each sample site calculated with CLUMPP in the case where for K=4. The white arrow indicates the location of the introduction of Redwhiskered bulbuls on Reunion Island in 1972.

Fig. 4

 4 Fig. 4 Comparison of invasion scenarios conducted with DIYABC on the microsatellite dataset to infer invasion history on Reunion. a) Topologies of the four scenarios compared with DIYABC. M: Mauritius, H: Oahu, RB: Reunion blue cluster, RO: Reunion orange cluster. The thick lines represent population of constant size. The thin lines represent the duration of bottlenecks (b1 and b2). The time scale is given on the right, t1 and t2 are the number of generations elapsed between the present and introduction events or splits, tb1 and tb2 represent the length of the bottlenecks b1 and b2 in number of generations. b) Posterior probabilities of the four scenarios calculated with the logistic regression approach implemented in DIYABC on the 1% closest data sets. Error bars represent 95% confidence intervals.

  ) R-W5 R-W4 R-W3 R-W2 R-W1 R-Intro R-L1 R-L2 R-L3 R-L4 RP-values under the adjusted significance threshold corrected following the Holm correction.

Fig. S3 :

 S3 Fig. S3: boxplots of the body size of individuals according to a) their sex, b) their geographic origin, and c) their sample site.

Fig. S6 :

 S6 Fig.S6: Plots of the pairwise genetic distance between sample sites on Reunion against the pairwise geographic distance between sample sites. The geographic distance between sample sites are either Euclidian distances (top graphic) or disatnces calculated following the coastline (bottom graphic). Gey dots represent pairwise comparisons between sites belonging to different clusters whereas orange and blue dots represent pairwise comparisons between sites belonging to the orange cluster and the blue cluster respectively.

Table 1

 1 Effect of geographic origin, sample site and sex on body size and beak size of the individuals (ANCOVAs, upper part of the table) and on body conformation and beak conformation of the individuals (MANCOVAs, lower part of the table).

	Size		Df	F value	P value
		Geo.origin:Site:Sex 14, 410	0.77	0.71
	Whole body	Geo.origin:Sex Geo.origin:Site	6, 424 14, 430	0.49 6.59	0.82 3.759e-12*
		Geo.origin	6, 430	39.33	< 2.20e-16*
		Sex	1, 430	80.45	< 2.20e-16*
		Geo.origin:Site:Sex 14, 427	1.23	0.25
		Geo.origin:Sex	6, 441	0.33	0.92
	Beak	Geo.origin:Site	14, 447	3.25	6.06e-05*
		Geo.origin	6, 447	52.16	< 2.20e-16*
		Sex	1, 447	41.56	2.96e-10*
	Conformation	Df	F value	P value
		Geo.origin:Site:Sex 14, 410	0.71	0.92
	Whole body	Geo.origin:Sex Geo.origin:Site	6, 424 14, 430	0.71 2.69	0.80 5.23e-08*
		Geo.origin	6, 430	25.29	< 2.20e-16*
		Sex	1, 430	58.69	< 2.20e-16*
	Beak	Geo.origin:Site	14, 263	2.46	3.41e-08*
	(Males)	Geo.origin	6, 263	9.13	< 2.20e-16*
	Beak	Geo.origin:Site	14, 165	1.56	7.35e-03*
	(Females)	Geo.origin	6, 165	7.63	< 2.20e-16*
	* P-values under the significance threshold.

Table 2

 2 Pairwise comparisons of body and beak size (Tuckey's HSD tests, upper part of the table) and pairwise comparisons of body and beak conformation (Hotelling's T-squared tests, lower part of the table) between the different geographic groups.

	Size		Western India	Himalaya Indochina	Eastern India	Reunion Mauritius
		Himalaya	x
		Indochina	<1.00e-05 <1.00e-05
	Whole body	Eastern India Reunion	x x	x x	<1.00e-05 <1.00e-05	x
		Mauritius	<1.00e-05	x	<1.00e-05 <1.00e-05 <1.00e-05
		Oahu	2.64e-02	1.94e-03	x	4.50e-03 <1.00e-05 <1.00e-05
		Himalaya	x
		Indochina	2.81e-03	2.85e-03
	Beak	Eastern India Reunion	x <1.00e-05 <1.00e-05 6.23e-04 <1.00e-05 x 1.92e-02
		Mauritius	x	x	<1.00e-05 1.02e-02 <1.00e-05
		Oahu	8.02e-03	x	<1.00e-05 1.36e-04 <1.00e-05	x
	Conformation		Western India	Himalaya Indochina	Eastern India	Reunion Mauritius
		Himalaya	x
		Indochina	x	x
	Whole body	Eastern India Reunion	x 0.00	x 7.87e-14	x 0.00	0.00
		Mauritius	0.00	1.52e-10 2.12e-12	0.00	0.00
		Oahu	0.00	3.37e-09 1.30e-13	0.00	0.00	0.52
		Himalaya	x
		Indochina	x	x
		Eastern India	x	x	x
	Beak (males)	Reunion	5.26e-11	3.26e-06	0.00	8.44e-15
		Mauritius	1.94e-05 5.01e-04 1.80e-09 6.53e-07 3.90e-05
		Oahu	x	x	x	x	5.34e-08	x
		Himalaya	x
		Indochina	x	x
	Beak	Eastern India	x	x	x
	(females)	Reunion	8.60e-10 2.90e-04	0.00	1.15e-12
		Mauritius	x	x	1.05e-07 3.18e-05 2.99e-04
		Oahu	x	x	5.23e-04	x	6.26e-11 3.24e-04

Table 3

 3 Pairwise comparisons of body and beak size (Tuckey's HSD tests, upper part of the table) and pairwise comparisons of body and beak conformation (Hotelling's Tsquared tests, lower part of the table) between the windward and leeward coasts for each island.

	Data Set	Island	P-value
	Size
		Reunion	0.00*
	Whole body	Mauritius	4.10e-06*
		Oahu	4.35e-03*
		Reunion	x
	Beak	Mauritius	x
		Oahu	x
	Conformation
		Reunion	3.51e-05*
	Whole body	Mauritius	1.78e-03*
		Oahu	x
		Reunion	5.22e-03*
	Beak (males)	Mauritius	x
		Oahu	x
		Reunion	x
	Beak (females)	Mauritius	x
		Oahu	x

Table 4

 4 Genetic diversity at microsatellite loci in the invasive populations for each sampling site.

			N	Na	I	Ho	uHe	F	HWE LD
	Oahu	O-K	9.80 (+/-0.20) 2.10 (+/-0.28) 0.57 (+/-0.11) 0.36 (+/-0.09) 0.39 (+/-0.07) 0.05 (+/-0.11)	0/8	3/28
		O-PC	34.6 (+/-0.27) 2.30 (+/-0.37) 0.59 (+/-0.14) 0.35 (+/-0.07) 0.37 (+/-0.07) 0.02 (+/-0.07)	1/8	1/28
		M-B	16.00 (+/-0.00) 3.20 (+/-0.39) 0.77 (+/-0.14) 0.40 (+/-0.08) 0.43 (+/-0.08) 0.09 (+/-0.07) 1/10 0/45
	Mauritius	M-C	22.00 (+/-0.00) 3.30 (+/-0.45) 0.83 (+/-0.16) 0.52 (+/-0.10) 0.47 (+/-0.08) -0.07 (+/-0.10) 1/10 1/45
		M-IA	11.90 (+/-0.10) 2.80 (+/-0.42) 0.74 (+/-0.15) 0.48 (+/-0.09) 0.45 (+/-0.08) -0.09 (+/-0.06) 0/8	1/36
		R-W5	78.80 (+/-0.47) 3.40 (+/-0.37) 0.74 (+/-0.13) 0.37 (+/-0.07) 0.42 (+/-0.07) 0.11 (+/-0.09) 3/10 4/45
		R-W4	30.60 (+/-0.27) 3.00 (+/-0.37) 0.70 (+/-0.14) 0.40 (+/-0.09) 0.40 (+/-0.08) 0.07 (+/-0.12)	1/9	2/45
		R-W3	33.20 (+/-0.42) 3.10 (+/-0.35) 0.75 (+/-0.15) 0.44 (+/-0.09) 0.43 (+/-0.08) 0.01 (+/-0.08) 2/10 2/45
		R-W2	12.90 (+/-0.10) 2.70 (+/-0.40) 0.74 (+/-0.16) 0.41 (+/-0.09) 0.45 (+/-0.09) 0.04 (+/-0.08)	0/8	0/36
		R-W1	19.90 (+/-0.10) 3.20 (+/-0.36) 0.73 (+/-0.13) 0.40 (+/-0.08) 0.42 (+/-0.07) 0.06 (+/-0.08)	3/9	3/45
	Reunion	R-Intro 13.90 (+/-0.10) 2.90 (+/-0.46) 0.75 (+/-0.15) 0.45 (+/-0.09) 0.45 (+/-0.08) -0.02 (+/-0.09) 1/8	0/36
		R-L1	27.60 (+/-0.87) 2.80 (+/-0.42) 0.73 (+/-0.16) 0.41 (+/-0.08) 0.43 (+/-0.09) -0.01 (+/-0.07) 0/8	4/36
		R-L2	36.70 (+/-0.62) 2.90 (+/-0.50) 0.72 (+/-0.18) 0.41 (+/-0.10) 0.41 (+/-0.09) -0.03 (+/-0.05) 1/8	1/36
		R-L3	31.40 (+/-0.27) 3.30 (+/-0.40) 0.81 (+/-0.13) 0.44 (+/-0.09) 0.47 (+/-0.07) 0.09 (+/-0.12) 3/10 6/45
		R-L4	28.40 (+/-0.22) 3.50 (+/-0.48) 0.80 (+/-0.15) 0.45 (+/-0.09) 0.45 (+/-0.08) 0.05 (+/-0.09) 3/10 3/45
		R-L5	52.80 (+/-0.61) 3.60 (+/-0.50) 0.81 (+/-0.14) 0.42 (+/-0.08) 0.46 (+/-0.08) 0.12 (+/-0.08) 3/10 4/45
		R-L6	10.80 (+/-0.20) 2.50 (+/-0.31) 0.67 (+/-0.14) 0.38 (+/-0.09) 0.42 (+/-0.09) 0.02 (+/-0.11)	1/7	0/36

Table S2 :

 S2 Sample sites names and locations in Reunion, Mauritius and Oahu.

	Island	Site ID Site Name	latitude longitude
		R-W5 Sainte-Suzanne	-20.91 55.60
		R-W4 Saint-André	-20.96 55.65
		R-W3 Bras-Panon	-21.01 55.69
		R-W2 Saint-Benoît	-21.04 55.72
		R-W1 Sainte-Rose	-21.13 55.79
	Reunion	R-Intro Saint-Philippe R-L1 Saint-Pierre	-21.36 55.76 -21.33 55.47
		R-L2	Le Tampon	-21.29 55.52
		R-L3	Saint-Louis	-21.29 55.41
		R-L4	Etang-Salé	-21.27 55.36
		R-L5	Les Avirons	-21.24 55.33
		R-L6	Saint-Leu	-21.17 55.29
		M-B	Bel Ombre field station -20.46 57.43
	Mauritius	M-C	Camp field station	-20.39 57.45
		M-IA Ile aux aigrettes	-20.42 57.73
	Oahu	O-K O-PC Pearl city Kaneohe	21.41 -157.80 21.42 -157.95

Table S3 :

 S3 primers used for amplification of COI and ND2 genes.

	Gene	Primer F		Primer R
			ACGCTTTAACACTCAGCCAT		ACGTGGGAGATAATTCCAAATC
	COI (single fragment)	COI-ExtF	CTTACC	COI-BirdR1	CTG
			GATGAYTATTTTCAACCAAC		CTYATGTTGTTTATTCGRGGGAA
	COI (1st fragment)	COI-ExtF1b	CACAAAGA	COI-R220b	AGC
			TTGGCGGATTYGGAAACTG		TGTGATAGGGCKGGGGGTTTTA
	COI (2nd fragment)	COI-F167b	ACTAGT	COI-R451b	TGTT
			GGTRTCTCCTCAATCTTAGG		GGTAGGATTAGGATATAGACTT
	COI (3rd fragment)	COI-F403b	AGCAAT	COI-R661b	CTGGATG
			AGCTATCGGGCCCATACCCC		TTGAAGGCCTTCGGTTTAGGTG
	ND2 (single fragment)	ND2-ExtF	GAA	ND2-ExtR	A
			AGCTATCGGGCCCATACCCC		TGTCCAGTGTACCATGCGTTGGT
	ND2 (1st fragment)	ND2-ExtF	GAA	ND2-R223	CA
			ACTTCTTGACCCAAGCAACA		GAATAGTAGCGTGATTGGGGGG
	ND2 (2nd fragment)	ND2-F175	GCCTCA	ND2-R405	AATT
			TGCAAGGATCCCCCCTTATy		GGGGTTGTAGGTGATGATGATG
	ND2 (3rd fragment)	ND2-F355	ACTGGA	ND2-R568	GCTA

Table S4 :

 S4 Description of microsatellite markers used in the study: name, number of alleles found in our study, sequence of the primers used for amplification, type of motif according to the literature, fluorescent label used for the amplification and original reference.

Table S6 :

 S6 Pairwise comparisons of body size (upper part of the table) and beak size (lower part of the table) between sample sites on Reunion conducted with Tuckey's HSD tests.

Table S7 :

 S7 Test for allometry: effect of geographic origin, sample site and sex on the conformation variables for the whole body data set (ANCOVAs, upper part of the table) and for the beak data set (ANCOVAs, lower part of the table). Conformation variables: first axes representing 95% of the total variance of the PCAs performed on measurements for the whole body data set and on geometric morphometric coordinates for the beak data set.

	Whole	Geo.origin:site:sex	Geo.origin:sex:	Geo.origin:site:	Geo.origin:	Size:	Geo.origin:	Geo.origin:
	body	:size	size	sex	site:size	sex	sex	size
	Axis 1	0.53	0.85	0.77	0.73	0.45	8.04e-02	2.28e-07*
	Axis 2	0.59	0.94	0.78	0.79	0.11	0.27	0.11
	Axis 3	0.31	0.79	0.48	0.23	0.45	0.50	0.29
	Beak	Geo.origin:site:sex :size	Geo.origin:sex: size	Geo.origin:site: sex	Geo.origin: site:size	Size: sex	Geo.origin: sex	Geo.origin: size
	Axis 1	0.83	0.66	0.39	0.89	0.89	0.85	0.87
	Axis 2	0.59	0.87	0.10	0.72	0.49	0.29	5.44e-02
	Axis 3	0.68	7.93e-02	7.97e-02	0.81	1.85e-02*	0.15	0.36
	Axis 4	0.72	0.60	0.26	5.51e-02	4.34e-02*	0.21	0.24

* P-values under the significance threshold. A significant interaction between global geographic origin / site and size (highlighted in bold) indicates that there is a difference in allometry between geographic groups.

Table S8 :

 S8 Pairwise comparisons of body conformation (upper part of the table) and beak conformation (lower part of the table) between sample sites on Reunion conducted with Hotelling's T-squared tests.

Table S9 :

 S9 Posterior probabilities of the scenarios compared with the ABC analysis of the microsatellite data set for each set of priors. These posterior probabilities were assessed using a polychotomous logistic regression on the 0.1 and 1% simulated data sets closest to the observed data set.

		Prior set 1	Prior set 1
	Scenario P 0.1% P 1% P 0.1% P 1%
	A	0.13	0.14	0.09	0.08
	B	0.10	0.11	0.06	0.07
	C	0.57	0.55	0.55	0.58
	D	0.21	0.21	0.30	0.27

Table S10 :

 S10 Probability distributions of the demographic parameters of scenario C based on the 1% of simulated data sets closest to the real data set.

	Parameter mean	median	mode	q025	q050	q250	q750	q950	q975
	NM	1.60e+003 1.44e+003 1.10e+003 5.25e+002 6.43e+002 1.06e+003 1.96e+003 3.04e+003 3.50e+003
	NRO	9.20e+003 8.70e+003 2.16e+003 8.24e+002 1.17e+003 4.14e+003 1.39e+004 1.87e+004 1.94e+004
	NRB	9.92e+003 9.81e+003 1.74e+003 9.77e+002 1.43e+003 4.91e+003 1.48e+004 1.89e+004 1.95e+004
	t2	2.58e+001 2.55e+001 2.24e+001 8.32e+000 1.03e+001 1.85e+001 3.29e+001 4.19e+001 4.42e+001
	tb2	2.33e+001 2.36e+001 2.49e+001 2.96e+000 5.28e+000 1.45e+001 3.20e+001 4.08e+001 4.28e+001
	NROb	2.80e+002 2.14e+002 1.04e+002 2.45e+001 3.72e+001 1.10e+002 3.84e+002 7.63e+002 8.64e+002
	t1	3.93e+001 4.09e+001 4.76e+001 2.06e+001 2.36e+001 3.43e+001 4.59e+001 4.93e+001 5.00e+001
	tb1	3.07e+001 3.25e+001 4.17e+001 6.43e+000 9.51e+000 2.32e+001 3.97e+001 4.61e+001 4.72e+001
	NRBb	3.24e+002 2.65e+002 1.37e+002 3.51e+001 5.18e+001 1.48e+002 4.47e+002 8.10e+002 9.02e+002

Acknowledgements

We thank the LabEx BCDiv and the ATM "Formes possibles, Formes réalisées -Exploration du vivant" which financed our field missions and molecular analyses. We thank the British Natural History Museum for giving us access to their collection of Red-whiskered bulbuls. We thank the Australian Museum, the British Natural History Museum, the California Academy of Sciences, the Field Museum of Natural History, the Muséum National d'Histoire Naturelle, the University of Michigan Museum of Zoology and the Victoria Museum, for the samples they lent us. We also thank Frederic Austerlitz, Paul Verdu and Arnaud Estoup for their advice on ABC analyses. Finally, we thank the molecular systematic platform (SSM) of the Muséum National d'Histoire Naturelle, its staff, and particularly Josie Lambourdière, which was of great support for the setting up of microsatellites amplification protocols and for genotyping.

P. jocosus NA Mauritius Bel Ombre 2013 Le Gros B07 KX529932 KX530000 P. jocosus NA Mauritius Bel Ombre 2013 Le Gros B08 KX529933 KX530001 P. jocosus NA Mauritius Bel Ombre 2013 Le Gros B09 KX529934 KX530002 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC01 KX529924 KX529992 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC02 KX529925 KX529993 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC03 KX529926 KX529994 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC04 KX529927 KX529995 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC07 KX529928 KX529996 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC08 KX529929 KX529997 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC09 KX529930 KX529998 P. jocosus NA USA (Hawaii) Pearl City 2013 Le Gros PC10 KX529931 KX529999 P. jocosus NA USA (Hawaii) Kailula 2013 Le Gros Ka01 KX529923 KX529991 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K01 KX529914 KX529982 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K02 KX529915 KX529983 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K03 KX529916 KX529984 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K04 KX529917 KX529985 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K05 KX529918 KX529986 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K06 KX529919 KX529987 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K07 KX529920 KX529988 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K08 KX529921 KX529989 P. jocosus NA USA (Hawaii) Kaneohe 2013 Le Gros K09 KX529922 KX529990 P. jocosus NA Australia NA 2009 AM O.72754 KX529895 KX529964 P. jocosus NA Australia New South Wales, Richmond 2007 AM O.71663 KX529894 KX529963 P. jocosus NA Australia New South Wales, Figtree 1987 AM O.59906 KX529892 KX529961 P. jocosus NA Australia New South Wales, Sydney 1971 AM O.43790 KX529891 KX529960 P. jocosus NA Australia New South Wales, Sydney 1991 AM O.62732 KX529893 KX529962 P. jocosus NA Australia New South Wales 1985 MV B926 KX529959 KX530027 P. jocosus NA Australia New South Wales 1985 MV B924 KX529958

Data Accessibility

The morphometric and microsatellite data supporting this manuscript are archived in Dryad (doi:10.5061/dryad.563n9). The DNA sequences have been deposited in GenBank (see table S1 for access numbers).

Author Contributions

All authors participated in the writing of this manuscript. ALG collected all the data, with the help of BM for Oahu. ALG did the phylogeographic analysis with the help of DZ, the morphometric analyses with the help of RC, and the population genetic analyses. PC and SS supervised this study.

Supplementary information

Appendix S1 -Amplification of mitochondrial genes

The following reagent quantities were used to amplify the COI and ND2 genes in a single fragment (primers, table S3): 13.94µL of Milli-Q water, 2µL of polymerase buffer 10x (Qiagen), 1.5µL of MgCl2 (Qiagen, 25mM), 1µL of Dimethyl sulfoxide (DMSO), 0.8µL of dNTPs mix (1.7mM each), 0.32µL of each primer (10µM), 0.12µL of DNA polymerase (Qiagen, Taq 5 units/µL). Cycling conditions: 94°C, 5 min.; (94°C, 40 sec.; 55°C, 40 sec.; 72°C, 60sec.) x 40 cycles; 72°C, 5 min.

For the amplification of short fragments (primers, table S3), illustra™ Hot Start Mixes (GE Healthcare) were used with the following reagent quantities: 19µl of Milli-Q water, 1µL of MgCl2 (Qiagen, 25mM), 1µL of each primer (10µM

Appendix S2 -Amplification of microsatellite loci

The ten microsatellite loci were amplified for each individual in two separate multiplex and tagged with fluorescent forward primers (dyes: 6-FAM, VIC, NED, PET; Applied Biosystems, table S4). PCR amplifications were done using the following reagent quantities: 1.25µL of the primer mix (1µM of each primer and TE buffer), 4µL of RNase-free water (Qiagen), 6.25µL of 2x Type-it Multiplex PCR Master Mix (Qiagen) in a final volume of 11.5µL. The following cycling conditions were used: 95°C, 5 min.; (95°C, 30 sec.; 57°C, 90 sec.; 72°C, 30sec.) x 25 cycles; 60°C, 30 min.

Appendix S3 -Phylogenetic analyses

The Bayesian inference was conducted with MrBayes 3.1.2 [START_REF] Ronquist | MrBayes 3: Bayesian phylogenetic inference under mixed models[END_REF].

In order to account for the potential differences in nucleotide substitution models between the data partitions corresponding to the two genes, a mixed model approach was implemented. MrModelTest 2.3 [START_REF] Nylander | MrModeltest v2. Program distributed by the author[END_REF]) and PAUP* [START_REF] Swofford | Phylogenetic Analysis Using Parsimony (*and Other Methods)[END_REF] were used to obtain the models best fitting the data, according to the AIC criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF]. Uniform interval priors were selected for the parameters, except for base frequencies, which were assigned a Dirichlet prior [START_REF] Huelsenbeck | MRBAYES: Bayesian inference of phylogenetic trees[END_REF]. Two independent runs of four incrementally heated Metropolis-coupled MCMC chains were run for 10 million generations. Sampling was done every 1000 generations, yielding 20000 trees. The online version of AWTY [START_REF] Nylander | AWTY (are we there yet?): a system for graphical exploration of MCMC convergence in Bayesian phylogenetics[END_REF] was used to assess the convergence of the MCMC chains and to estimate the "burn-in" length (2000 trees). Maximum likelihood searches of the partitioned dataset were conducted with RAxML v. 7.0.3 [START_REF] Stamatakis | RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models[END_REF]) using a GTR+Γ+I model and a random starting tree. The α-shape parameters, GTR-rates, and empirical base frequencies were estimated and optimized for each partition. Nodal support was estimated using 100 bootstrap replicates. Nx: number of individuals in the population x (constant in time), Nxb: number of individuals in the population x during the bottleneck following introduction (constant during the whole bottleneck), ti: number of generations elapsed between the present and an introduction event or a split i, tbi: number of generations elapsed during the bottleneck following an event i.

Supplementary tables

Supplementary figures


Bonabe De Rouge
email: r.bonabederouge@efficacity.com

Stabat Pascal

Marchio Dominique

A comparison of integration solutions for a gas Stirling micro-cogeneration system in residential buildings

Micro-Combined Heat and Power (µCHP) is a technology allowing to produce distributed electricity (electrical power below 36 kWel). The main advantage of µCHP is to avoid thermal energy losses while producing and distributing electricity in order to use primary energy rationally. The electricity can be consumed locally in the building or in the neighborhood area while thermal energy can be recovered to satisfy space heating (SH) and domestic hot water (DHW) demands. Stirling Engine (SE) is one of the various µCHP technologies that can be integrated in individual or small multi-family residential buildings. The aim of this paper is to study the impact of storage configurations, hydraulic schemes and control strategies on the performance of a µCHP system. This study is carried out on an individual residential building with various electrical, DHW and SH loads. A dynamic model was developed to represent a 1 kWel gas Stirling µCHP system (producing up to 24 kWth thanks to an auxiliary burner) and has been experimentally calibrated and validated. The model and storage, emitter and distribution network models are implemented in DYMOLA/Modelica environment in order to get a dynamic representation of the whole installation coupled to the building. The comparison between the different cases is shown in terms of primary energy consumption and analysis of operating conditions. The results allow to identify the best SE-µCHP designs according to building conditions.

INTRODUCTION

Reduction of energy use in buildings is a challenge to reduce energy consumption and lower greenhouse gases emissions. Many actions are made to reduce the energy demand of buildings such as improving wall insulation and air tightness. A lot of work can also be done to make the production of electricity and thermal energy for buildings more efficient. Micro-Combined Heat and Power systems allow to produce simultaneously on-site electrical and thermal power to fulfill the needs of buildings. It results in high primary energy efficiency electricity production unlike centralized power plants. The produced electricity can be selfconsumed if there is a concomitant electrical need or exported to the electrical network. However, the thermal power needs to be totally consumed in the building. An auxiliary burner (AB) is usually integrated in SE-µCHP systems for peak heat demand. The sizing of a microcogeneration system is thereby critical as it is needed to produce as much electricity as possible (to reduce the primary energy demand of the building), for a given thermal need.

Furthermore, most of the micro-cogeneration systems must deal with slow dynamics. Solid-Oxyde Fuel Cells for example may take several hours to reach full capacity [START_REF] Beausoleil-Morrison | Experimental Investigation of Residential Cogeneration Devices and Calibration of Annex 42 Models[END_REF] and start-stop cycles need to be avoided. For micro Stirling Engines (SE), such as the one studied in this paper, the full capacity of the system is reached around 15 minutes after start-up [START_REF] Bouvenot | Gas Stirling engine micro CHP boiler experimental data driven model for building energy simulation[END_REF]. Decreasing number of start-ups during the year can result in increasing global efficiency. Furthermore, in France, feed-in tariff for exported electricity to the grid is lower than purchase price and requires today a high investment cost for an additional electric meter. Thus, it seems more interesting to increase self-consumption of produced electricity. Designers can rely on many possibilities to get optimal performance from µCHP systems: indeed, storage tanks can be used to store thermal energy for space heating and domestic hot water needs. Control can also be adjusted to longer operating phases and reduce the need for auxiliary burner to increase the share of the µCHP in thermal energy production and so produce more electricity. Simulations and experiments also showed that both electrical and thermal needs [START_REF] Dorer | Performance assessment of fuel cell micro-cogeneration systems for residential buildings[END_REF] influence the global primary energy savings of the µCHP system and the profitability of the installation regarding pricing conditions of electricity export [START_REF] Alanne | Techno-economic assessment and optimization of Stirling engine micro-cogeneration systems in residential buildings[END_REF]. Authors have often studied the impact of building and occupant needs regarding the technology of micro-cogeneration or the optimal size for storage [START_REF] Dorer | Energy and CO 2 emissions performance assessment of residential micro-cogeneration systems with dynamic whole-building simulation programs[END_REF] but they rarely explored the impact of hydraulic storage configuration on primary energy savings or profitability. In addition, it was rarely possible to detail the complexity of controls of such systems. For such a study, yearly dynamic simulation seems to be the best approach. This paper relies on detailed dynamic modeling of systems, buildings, controls and occupants. Implementation of manufacturer control rules has been made and three hydraulic patterns are explored.

METHODS

This section presents the main hypothesis and the methods to model the µCHP system integrated in the building with distribution and emission of thermal energy.

Modeling tool

Dymola [START_REF]Dymola User Manual[END_REF] is a modeling and simulation tool based on Modelica language which is an acausal object-oriented programing language. A graphic interface allows to represent and assemble the different models of components. The strength of using Modelica is that is enables to fully represent interactions between all the systems and its controls with variable timestep solver.

Modeling of buildings

The building is modeled thanks to a resistor-capacitance grey box model with two nodes for light (air and furniture) and heavy capacities (walls and floors). Six thermal resistances represents walls and windows conductive resistances, convective resistances and ventilation [START_REF] Berthou | Development of building models for load curve forecast and design of energy optimization and load shedding strategies[END_REF]. This grey-box model is calibrated on air temperature, heating, power and annual heating energy needs thanks to "Design" calibration library [START_REF]Dymola User Manual[END_REF] against the detailed building model of Buildings Library [START_REF] Wetter | Modelica Buildings Library[END_REF]. This simplified building model has been selected instead of a detailed model in order to assure a faster simulation (about 10 times faster). The calibration is made for three levels of building performance, intending to represent high, medium and low levels of consumption for the area of Paris, France. Targets of yearly energy consumption are based on past and actual French thermal regulations. The characteristics of the buildings are given in Table 1. Although it is simplified, this building model enables to represent the building envelope interacting with systems and especially the transient behavior of hot water loop. Heat recovery ventilation 100 m 3 .h -1 Modeling of occupancy Modeling the human behavior seems to be critical in building simulation and especially with cogeneration systems for multiple reasons. First, the humans interact with the temperature setpoint in the building. They also generate heat, consume electricity and domestic hot water. Since the concomitance of electrical and thermal needs is critical for µCHP system, proper modeling is needed. Standard estimations of human behavior as well as averaged daily consumption of electricity seem inadequate for such a study. For residential applications, the human behavior in buildings has been studied such as Wilke [START_REF] Wilke | Probabilistic Bottom-up Modelling of Occupancy and Activities to Predict Electricity Demand in Residential Buildings[END_REF] and Vorger [START_REF] Vorger | Study of the influence of the inhabitants behavior on the energy performance of buildings[END_REF] studies. They both relied on statistical studies of occupant's activities in residential buildings. Thanks to their work, it is possible to represent a number of occupants and activities with a reduced time step of 10 minutes. Two sets of data (low and high) were created for two different families (2 and 4 persons). Scenarios are made with night and day setbacks (20 °C during day with occupation and 16 °C during inoccupation and night).

Modeling of electrical consumption

For each family, a weekly electricity consumption profile was used with a 10 minutes time step which proved to be adapted to represent residential electricity consumption [START_REF] Hawks | Impacts of temporal precision in optimisation modelling of micro-Combined Heat and Power[END_REF]. The weekly profiles are created thanks to bottom-up method based on previously exposed work. The electrical consumption is linked to each activity of the occupant and the energy class of their appliances. Aggregating every consumptions allows to generate a stochastic electrical pattern for a full week. As the average of non-HVAC consumption of electricity in French households is about 2700 kWh/yr [START_REF] Enertech | Mesure de la consommation des usages domestiques de l'audiovisuel et de l'informatique -Projet REMODECE[END_REF], two profiles are chosen below (1600 kWh/yr) and above (4500 kWh/yr) the mean. The study also uses annual lighting scenario for each family, depending on occupancy and time of the year to take into account the reduction of lighting during summer period. HVAC consumptions are added to the previous profiles to represent the consumption of Controlled Mechanical Ventilation (CMV) and of the additional space heating water pump when it is needed, according to the hydraulic configuration. The electrical consumption of the heat generation system is modeled in the µCHP model. An example of the daily aggregated consumption of the "high consumption" household for Tuesday is shown below: Figure 1. An example of daily electricity consumption with 10 minutes time step Modeling of domestic hot water consumption Similarly, domestic hot water profiles are generated with a 10 minutes time step and they are consistent with occupancy and activities. The water is tapped with a flow rate of 1.2 10 -5 m 3 /s for each occupant using hot water. The profiles are generated for each of the two families. The "low consumption" household consumes 89 L/day of water at 60 °C while the "high consumption" one consumes 166 L/day of water at 60 °C. These values can be compared reasonably to those of the multiple studies cited in [START_REF] Lagogue | PACTE ECS Les besoins d'ECS en habitat[END_REF] The Table 2 summarizes the characteristics of the two "families": Table 2 The water exiting the µCHP system is directly sent to the emitters. A storage is present for DHW needs. The auxiliary burner can help the SE only for SH needs as it is assumed that the engine is powerful enough to always manage DHW storage. Weather dependent setpoint (WDS) for space heating water is computed. DHW storage temperature is constant. Multiple sizes of DHW storages are tested with an insulation of 4 cm.

 Configuration 2 (C2)

A common storage is used for both DHW and SH needs.

The DHW passes through the tank to exchange with the stored water. The stored water circulates across the emitters. The AB can help the SE for both DHW and SH needs as they are not differentiated in the tank. An additional pump is needed and WDS is not possible as the stored water always needs to be at a high temperature to prepare DHW. Multiple sizes of storage are tested. It should be noted that the storage tank is very well insulated (13 cm). In this case, DHW and SH storage are separated. The 3-way valve controls which of the tanks is filled with hot water coming from the µCHP system. The AB can help the SE only for SH needs as it is assumed that the engine is powerful enough to always manage DHW storage. An additional pump is needed and WDS is possible as DHW and SH are physically separated. Multiple sizes of SH storage will be studied (DHW tank volume is fixed) with an insulation of 4 cm for each tank.

For start-up and stop phases, calibrated exponential correlations model evolution of fuel, thermal and electrical power as a function of time. Delays for electrical and thermal production after the beginning of the start-up are also taken into account. The auxiliary burner is modeled thanks to state of the art modulating boiler [START_REF] Bonabe De Rougé | Modélisation de systèmes de chauffages en résidentiel[END_REF]. It is also based on experimentation and performances depends on the same parameters and part load ratio (PLR). The auxiliary burner and the reference case boiler have a nominal power of 24 kWth. Finally, auxiliary power for the pump and controls electrical consumptions are modeled depending on operating conditions.

Primary energy and reference configuration

The primary energy factor are 1 for natural gas and 2.58 for electricity the value according to the French thermal regulation. The "reference" configuration used to compute primary energy savings is based on a condensing boiler with a 160 L DHW storage and a WDS for SH.

RESULTS

The different configurations can be compared thanks to various criteria. For microcogeneration, self-consumption (share of the produced electricity which is self-consumed) is a significant criterion as it is today in the French context better to consume electricity on-site than export it to the grid. Self-production (share of the electricity needs which is directly covered by the micro-cogeneration) is also an important criterion as it means that the building is not importing electricity from the grid for a part of its needs.

Figure 5 shows self-consumption and self-production ratios for three levels of building insulation and three hydraulic configurations and for "high" occupancy profile. For "low" occupancy profile, two hydraulic configurations and two levels of insulation are tested. Tank volumes are set to the "default" values. It can be seen that self-consumption mainly depends on the electricity load curve of the buildings. The simulations with "Low" electrical profile leads to slightly lower self-consumption ratio (42-47 %) than "High" profile (61-73 %). Figure 5 also shows the best self-consumption ratio is achieved for the building with higher thermal needs with hydraulic configuration No. 3. For both "Middle" and "Low" building cases, the best self-consumption ratios are reached for configuration No. 1. Detailed analysis of the simulations results shows that for these cases, more thermal storage leads to drastically shift the thermal load curve of the building so that thermal and electricity needs are less concomitant. In the best case, for the "Low" occupancy profile with "High" thermal needs, self-production ratio can reach 34 %.

Figure 7 shows Primary energy savings (PES) for various tank sizes, hydraulic configurations building insulations and occupancy profiles. X-axis represents E/H demand ratio of the building. The lowest the E/H demand ratio is (on the left) the highest the thermal needs are and the lowest the electricity needs are. "Low" occupancy profiles are shown in dotted lines. Figure 7 shows that the "best" point for primary energy savings (PES) seems to be achieved Presence of thermal storage for space heating (C2 and C3) allows to reach higher PES for low thermal needs (on the right of the chart) as it allows to longer the Stirling Engine operating cycles. However, for building with high thermal needs, thermal storage for SH is not needed because it does not help to meet higher savings. Configuration No. 2 shows low storage heat losses (see Figure 6) because of its higher insulation. For low thermal needs, high insulation and high storage capacity are primordial because losses can accounts for a large part of the total thermal needs and high storage volume lengthen operating cycles.

For "Middle" and "Low" building cases, using properly sized thermal storage for space heating demand enables to lower the use of the auxiliary burner as can be seen on figure 9 and consequently raises the share of the Stirling Engine in total heat production.

For "High" building case with configuration No. 2 as both SH and DHW needs are to be handled in the same tank and storage temperature setpoint is always higher than for other configurations (no WDS), the SE is not powerful enough when the tank is discharged and auxiliary burner is often needed.

In every cases, configuration without SH storage leads to increased numbers of Stirling Engine ON/OFF cycles as expected. This could lead to less efficient production but, as return temperature is lower than for storage configurations (2 and 3), global efficiency remains correct: both electrical and thermal efficiency increase for low temperature of returning water. If the number of operating cycles is critical for life expectancy of the µCHP, C2 may be preferred as the larger volume enables the longest operating cycles.

DISCUSSION

This paper shows that design of hydraulic integration clearly influences performances and operating conditions of a SE-µCHP system. First, it seems that DHW storage should be chosen close to the daily needs of hot water as increasing its size does not really enhance performance but grow thermal losses. In the case of this 1 kWel SE, the best configuration for buildings with high heating needs seems to be without SH storage. Thermal needs are large enough to allow long operating cycles and storage only creates heat losses. However, for buildings with low thermal needs, SH storage leads to higher primary energy savings as it allows longer operating cycles and then more electricity to be produced. It also appears that the best primary energy savings and self-production are reached for buildings for which E/H demand ratio is close to the nominal E/H ratio of the µCHP. From a strict economic point of view, solutions without thermal storage may still be preferred if the additional cost for storage tank cannot be refunded during the lifecycle of the µCHP system: this is why a complete economic analysis should be carried on. Parameters that could possibly influence the results are the settings of system controls. For example, larger hysteresis bands may be tested to enable longer operating cycles. Management of auxiliary burner could be explored as well as storage tank piping configurations. Exhaustive sensitivity analysis of the results over occupant's electrical and thermal needs should be assessed. Finally, the presented methodology should be extended to explore the best couples of various µCHP technologies over different buildings types and occupations.

2

 2 . Description of the two scenarios of occupants Level of consumption High Low Number of persons in the household 4 Consumption of non-HVAC electricity except lighting 4500 kWh/yr 1600 kWh/yr Quantity of tapped water at 60 °C per day 166 L/day 89 L/day Hydraulic configurations and controls Three hydraulic configurations are compared in this paper. The term "hydraulic configuration" means number of thermal storages and their type (SH or DHW storage). These configurations are based on µCHP manufacturers. All storage tanks are selected in some manufacturers' catalogues to get real data such as geometry and insulation. The emitters are always dimensioned for 60/50 °C temperature and nominal flow rate is adjusted for maximum heating needs. They are modeled with an equivalent radiator dynamic model and thermostatic radiator valve (TRV). The auxiliary burner can modulate its thermal power from 3.5 to 20 kW. It can deliver thermal power only if the Stirling engine has been working for more than 5 minutes. Two hysteresis cycles determine the temperature at which the Stirling engine and the auxiliary burner will start and stop.  Configuration 1 (C1)

Figure 2 .Figure 3 .

 23 Figure 2. Configuration with DHW storage only

Figure 4 .

 4 Figure 4. Configuration with two storages

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Self-consumption and self-production for various configurations

Figure 8 .Figure 9 .

 89 Figure 8. Auxiliary burner consumption for various configurations

Table 1 .

 1 Description of the three levels of performance for the use-case building

	Level of consumption	High	Medium	Low
	Area	90 m² (+ unheated basement and attic)
	SH annual consumption	172 kWh.m -2 .yr -1	75 kWh.m -2 .yr -1	35 kWh.m -2 .yr -1
	Type of windows	Simple glazing	Double glazing	Triple glazing
	Maximum heating needs (-7 °C)	7000 W	3700 W	2700 W
	Ventilation +	Single flow 150	Single flow 100
	infiltration	m 3 .h -1	m 3 .h -1

Table 3

 3 summarizes the studied sizes of tanks for each hydraulic configuration. "Default" value is shown in bold characters.

Table 3 .

 3 Description of the tested hydraulic configurations

	Configuration/Tank sizes in litres	DHW	SH
	1	160/200/400	-
	2	450/600/750
	3	160	100/200/400
	𝑃 𝑓𝑢𝑒𝑙 = 𝑃 𝑓𝑢𝑒𝑙 𝑛𝑜𝑚 + 𝑎 (𝑇 𝑐𝑤,𝑖 -𝑇 𝑐𝑤,𝑖 𝑛𝑜𝑚 ) + 𝑏 (𝑚̇𝑐 𝑤 -𝑚̇𝑐 𝑤 𝑛𝑜𝑚 )	(1)
	𝑄 ̇𝐻𝑋 = 𝑄 ̇𝐻𝑋 𝑛𝑜𝑚 + 𝑐 (𝑇 𝑐𝑤,𝑖 -𝑇 𝑐𝑤,𝑖 𝑛𝑜𝑚 ) + 𝑑 (𝑚̇𝑐 𝑤 -𝑚̇𝑐 𝑤 𝑛𝑜𝑚 ) 𝑃 𝑔𝑟𝑜𝑠𝑠 = 𝑃 𝑔𝑟𝑜𝑠𝑠 𝑛𝑜𝑚 + 𝑒 (𝑇 𝑐𝑤,𝑖 -𝑇 𝑐𝑤,𝑖 𝑛𝑜𝑚 ) + 𝑓 (𝑚̇𝑐 𝑤 -𝑚̇𝑐 𝑤 𝑛𝑜𝑚 )	(2)

Modeling of micro-cogeneration system

The 1 kWel + 7 kWth Stirling Engine is modeled thanks to the work of Bouvenot and Andlauer

[START_REF] Bouvenot | Gas Stirling engine micro CHP boiler experimental data driven model for building energy simulation[END_REF]

[START_REF] Andlauer | Optimisation systémiques de micro-cogénérateurs intégrés au bâtiment[END_REF]

. The semi-empirical model relies on experimental calibration and validation done at the laboratory of INSA Strasbourg. It was initially implemented as a TRNSYS type and it has been translated to Dymola environment for the purpose of this study. It represents the transient and stationary performance of the Stirling Engine according to its operating conditions, especially the mass flow rate of cooling water (𝑚̇𝑐 𝑤 ) return cooling water temperature (𝑇 𝑐𝑤,𝑖 ). Stationary performances are given by correlations for fuel inlet power (𝑃 𝑓𝑢𝑒𝑙 ), delivered thermal (𝑄 ̇𝐻𝑋 ) and electrical (𝑃 𝑔𝑟𝑜𝑠𝑠 ) power, such as:

ACKNOWLEDGEMENT

The authors would like to thank Efficacity for its support in the achievement of this project.