Aline Muyle
email: amuyle@uci.edu

Hélène Martin
email: helene.martin@tzone.org

Niklaus Zemp

Maeva Mollion

Sophie Gallina

Alexandre Raquel Tavares

Da Silva

Thomas Bataillon

Alex Widmer

Sylvain

Maéva Mollion

Alexandre Silva

Sylvain Glémin

Pascal Touzet
email: pascal.touzet@univ-lille1.fr

Gabriel Ab Marais
email: gabriel.marais@univ-lyon1.fr

Dioecy is associated with high genetic diversity and adaptation rates in the plant genus Silene

de 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

Out of 261,750 angiosperm species, approximately 15,600 (~6%) are dioecious [START_REF] Renner | The relative and absolute frequencies of angiosperm sexual systems: dioecy, monoecy, gynodioecy, and an updated online database[END_REF].

Dioecy is thus rare compared to the dominant sexual system in angiosperms, hermaphroditism, in which individuals carry bisexual flowers. Nonetheless, dioecy can be found in ~40% of angiosperm families. The current view is that many independent (871 to ~5000) and mostly recent transitions occurred from hermaphroditism to dioecy in angiosperms [START_REF] Renner | The relative and absolute frequencies of angiosperm sexual systems: dioecy, monoecy, gynodioecy, and an updated online database[END_REF]). These events probably followed two main evolutionary paths: one through a gynodioecious intermediate (with female and hermaphrodite individuals), called the gynodioecy pathway, and another one through a monoecious intermediate (with a single individual carrying separate female and male flowers) called the monoecy or paradioecy pathway [START_REF] Barrett | The evolution of plant sexual diversity[END_REF]. In few cases, dioecy might have evolved through alternative pathways, such as heterostyly (where individuals have different style and stamen lengths, [START_REF] Barrett | The evolution of plant sexual diversity[END_REF].

Species richness was originally found to be lower in dioecious compared to non-dioecious sister clades, suggesting that dioecy might lead to an increased extinction rate in angiosperms and defining the evolutionary dead-end hypothesis [START_REF] Heilbuth | Population genomic consequences of separate sexes in plants[END_REF]. Theoretical work has explored this idea and shown that an increased extinction rate could result from two main evolutionary handicaps.

First, in dioecious species, only females contribute to seed production and dispersal, which could result in individuals being less dispersed spatially and competing more for local resources compared to populations with hermaphrodites [START_REF] Heilbuth | The consequences of dioecy for seed dispersal: modeling the seed-shadow handicap[END_REF]. This should cause dioecious species to have smaller geographic distributions and species abundance compared to non-dioecious species.

Second, in dioecious species with insect pollination, male intra-sexual competition would result in males being more attractive than females for pollinators, making reproduction and population size sensitive to variation in pollinator density [START_REF] Vamosi | When looks can kill: the evolution of sexually dimorphic floral display and the extinction of dioecious plants[END_REF]. This should lead to fluctuating population sizes in dioecious species and in turn affect levels of genetic diversity. It has been proposed that these handicaps could be the reason why dioecy is associated with a specific suite of traits in angiosperms (fleshy fruits, wind-pollination, inconspicuous flowers and woody growth) as well as ecological preferences because dioecious species are preferentially found in the tropics and on islands [START_REF] Vamosi | Phylogenetic analysis of the ecological correlates of dioecy in angiosperms[END_REF]. Assuming that pollination mode evolves more slowly than dioecy, dioecy could be more frequent among wind-pollinated species because wind pollinated-species are not subject to the pollinator handicap. Similarly, assuming fruit type evolves more slowly than dioecy, dioecy could be more prevalent in species with fleshy fruits because seeds can be dispersed over longer distances by fruit-eating animals, thus limiting the seed-dispersal handicap.

The evolutionary dead-end hypothesis has recently been challenged [START_REF] Kafer | On the rarity of dioecy in flowering plants[END_REF], in part because the sister clade analysis used in early work assumed that speciation events and trait transitions coincide [START_REF] Heilbuth | Population genomic consequences of separate sexes in plants[END_REF]. However, if dioecy evolved after the split between the two diverging lineages and was not coincident with the split, then species richness is not expected to be equal between sister clades even when diversification rates are similar (Kafer and Mousset 2014).

Using a corrected version of the sister clade analysis on an updated dataset, dioecious clades have been found to have actually higher species richness than non-dioecious sister clades in angiosperms (Kafer et al. 2014). Using a phylogenetic method to study trait evolution and species diversification on several genera with dioecious species [START_REF] Maddison | Estimating a binary character's effect on speciation and extinction[END_REF], dioecy was not associated with an increased extinction rate [START_REF] Sabath | Dioecy does not consistently accelerate or slow lineage diversification across multiple genera of angiosperms[END_REF][START_REF] Duan | Evolution of sexual systems and growth habit in Mussaenda (Rubiaceae): Insights into the evolutionary pathways of dioecy[END_REF]. Another hypothesis to explain the rareness of dioecy in angiosperms posits that reversion from dioecy to other sexual systems is frequent and that dioecy is more evolutionary labile than previously thought [START_REF] Kafer | On the rarity of dioecy in flowering plants[END_REF].

This hypothesis was supported by a study of 2145 species from 40 genera in which the authors found evidence of frequent reversions from dioecy to other breeding systems [START_REF] Goldberg | Macroevolutionary synthesis of flowering plant sexual systems[END_REF]).

However, these works relied on macro-evolutionary and comparative approaches and the population functioning of dioecious versus non-dioecious species has so far been poorly explored.

Investigating consequences both at the macro-evolutionary and population scale has been fruitful in the study of the evolution of selfing from outcrossing, another frequent breeding system transition in angiosperms. Selfing was also proposed to be an evolutionary dead-end, which has been supported by several macroevolutionary analyses [START_REF] Goldberg | Species Selection Maintains Self-Incompatibility[END_REF][START_REF] De Vos | Heterostyly accelerates diversification via reduced extinction in primroses[END_REF][START_REF] Zenil-Ferguson | Interaction among ploidy, breeding system and lineage diversification[END_REF]. In parallel, various negative genetic consequences of selfing (i.e. low genetic diversity, accumulation of deleterious mutations) have been reported in several species, even in those of recent origin (see [START_REF] Glémin | Genome Evolution in Outcrossing vs. Selfing vs. Asexual Species[END_REF] for a review). Interestingly, negative genetic consequences have been detected despite the short terms advantage of selfing provided by reproductive assurance, allowing higher colonization potential and larger species range [START_REF] Randle | Can differences in autonomous selfing ability explain differences in range size among sister-taxa pairs of Collinsia (Plantaginaceae)? An extension of Baker's Law[END_REF][START_REF] Grossenbacher | Geographic range size is predicted by plant mating system[END_REF].

If the seed-dispersal handicap of dioecy operates, it should reduce the ability of dioecious plants to spread and increase local competition, which should lead to smaller geographic distributions and lower species abundance in dioecious compared to non-dioecious species. These expectations can be tested using estimates of geographic distributions and species abundance in dioecious and non-dioecious species. If the pollinator handicap of dioecy exists, it should cause fluctuations in dioecious population sizes. The effective population size (N e ) is given by the harmonic mean of the varying population sizes over time [START_REF] Wright | Size of a population and breeding structure in relation to evolution[END_REF]. The harmonic mean tends to be dominated by the smallest population sizes that the population goes through. Fluctuating population sizes are therefore expected to decrease the long-term N e of a species. If the pollinator and the seed-dispersal handicap of dioecy exist, they should cause dioecious species to have smaller N e compared to non-dioecious species. This expectation can be tested by studying genetic diversity Population genomic consequences of separate sexes in plants to estimate the long-term N e of dioecious and non-dioecious species. A reduced effective population size should then lead to a decrease in the efficacy of selection. This expectation can be tested by comparing closely related dioecious and non-dioecious species. A study that used the ratio of nonsynonymous to synonymous substitution rate (d N /d S ) as a proxy to selection efficacy found that one out of two dioecious Silene clade had a reduced efficacy of selection for some genes, but positive selection increased in other genes [START_REF] Kafer | Patterns of molecular evolution in dioecious and non-dioecious Silene[END_REF]. However, transitions to dioecy can happen after speciation, making d N /d S biased because it is computed on entire branches. We therefore decided to use instead a population genetics approach to confront predictions of the dead end hypothesis in Silene. We relied on statistics that allow quantifying the strength of natural selection, such as rates of adaptive and deleterious nonsynonymous substitutions. If dioecy is not a dead-end, as suggested by macro-evolutionary analyses, we expect no effect (or even a positive effect) of dioecy on genetic diversity and the efficacy of selection. Note that if extinction due to the handicaps of dioecy is too rapid to allow for populations to survive long enough to be sampled, our approach would not be able to detect the genetic footprints expected under the dead-end hypothesis.

However, as mentioned above, similar approaches in selfing species have consistently detected the predicted effects, even in very recent species such as Capsella rubella that originated about 20,000 years ago from the self-incompatible C. grandiflora [START_REF] Brandvain | Genomic Identification of Founding Haplotypes Reveals the History of the Selfing Species Capsella rubella[END_REF][START_REF] Slotte | The Capsella rubella genome and the genomic consequences of rapid mating system evolution[END_REF].

In this study, we focused on 13 species of the large genus Silene (~700 species) in the Caryophyllaceae family, in which dioecy has evolved independently at least twice [START_REF] Desfeux | Evolution of reproductive systems in the genus Silene[END_REF]. The Silene genus is considered a model system to study the evolution of sexual systems in angiosperms, thanks to its reasonably well-described sexual systems and phylogeny [START_REF] Bernasconi | Silene as a model system in ecology and evolution[END_REF]. Dioecy in Silene is likely to have evolved through the gynodioecy pathway [START_REF] Fruchard | The evolution of sex determination in plants[END_REF]. If dioecy is associated with the handicaps described above, dioecious Silene species should suffer from them because they do not have any of the traits named previously that likely compensate or prevent these handicaps. Indeed, Silene are short-lived herbs, many of them (but not all) are insect pollinated, they do not have fleshy fruits and they live in temperate regions. The Silene genus is thus adequate to test for the dead-end hypothesis of dioecy.

We sampled dioecious species from two groups in which dioecy evolved independently (5 species in the Melandrium section and 2 species in the Otites section) along with several hermaphroditic and gynodioecious close relatives for each group and one outgroup (Dianthus chinensis). Silene species can have very large genome sizes [START_REF] Bernasconi | Silene as a model system in ecology and evolution[END_REF]) and no complete reference genome is available in the genus (only a fourth of the Silene latifolia genome has been assembled to date [START_REF] Papadopulos | Rapid Y degeneration and dosage compensation in plant sex chromosomes[END_REF]. We therefore used a population genomics approach relying on RNA-seq data, which has proven successful in organisms that lack a reference genome [START_REF] Cahais | Reference-free transcriptome assembly in non-model Population genomic consequences of separate sexes in plants animals from next-generation sequencing data[END_REF][START_REF] Gayral | Reference-free population genomics from nextgeneration transcriptome data and the vertebrate-invertebrate gap[END_REF][START_REF] Romiguier | Comparative population genomics in animals uncovers the determinants of genetic diversity[END_REF]. Multiple individuals were Population genomic consequences of separate sexes in plants sampled across the species' geographic ranges and sequenced by RNA-seq. A pipeline specifically designed for such data was applied to perform SNP-calling [START_REF] Tsagkogeorga | The population genomics of a fast evolver: high levels of diversity, functional constraint, and molecular adaptation in the tunicate Ciona intestinalis[END_REF][START_REF] Gayral | Reference-free population genomics from nextgeneration transcriptome data and the vertebrate-invertebrate gap[END_REF][START_REF] Nabholz | Transcriptome population genomics reveals severe bottleneck and domestication cost in the African rice (Oryza glaberrima)[END_REF]). The SNPs were then analyzed to estimate genetic diversity and the efficacy of selection in each species.

Results

Phylogeny, geographic range, species abundance and population genetics characterization of the sampled species Our species sampling (13 species in total, see Fig. 1) includes two hermaphroditic species (S. viscosa and S. paradoxa), three gynodioecious species (S. vulgaris, S. nutans W and S. nutans E) and seven dioecious species (section Melandrium: S. heuffelii, S. dioica, S. marizii, S. diclinis, S. latifolia and section Otites: S. pseudotites and S. otites) from subgenera Behenantha and Silene, the main subgenera in the Silene genus [START_REF] Jenkins | A phylogenetic comparative study of preadaptation for invasiveness in the genus Silene (Caryophyllaceae)[END_REF]. The criteria for selecting the nondioecious species are detailed in the Materials and Methods. Briefly, we selected (i) close relatives of both dioecious sections based on existing phylogenies, (ii) species for which the sexual system is well described, (iii) species that had similar geographic ranges and (iv) hermaphroditic species not exhibiting high levels of selfing. The latter criterion was meant to perform a fair comparison between dioecious and hermaphroditic species, not confounded by strong differences in outcrossing rates.

For each species, between 4 and 34 individuals were sampled from different populations covering the geographic range of the species (see Fig. 1 for the number of individuals sampled in every species and Fig. S5 for sample locations). Previous work has shown that reliable estimates of population genetics statistics can be obtained when individuals from populations spanning the whole geographic distribution of a species are sampled, even if the number of sampled individuals is small, i.e. 2 or 3 [START_REF] Cahais | Reference-free transcriptome assembly in non-model Population genomic consequences of separate sexes in plants animals from next-generation sequencing data[END_REF][START_REF] Gayral | Reference-free population genomics from nextgeneration transcriptome data and the vertebrate-invertebrate gap[END_REF][START_REF] Romiguier | Comparative population genomics in animals uncovers the determinants of genetic diversity[END_REF].

Seeds were collected in the wild and plants were grown in controlled greenhouse conditions. Flower buds (and in a few cases, leaves) were sampled for 130 individuals in total and sequenced by RNAseq (see Table S1 for details on samples). De novo reference transcriptomes were assembled separately for each species and ORFs and orthologs were annotated (see Materials and methods and Table S1). All dioecious Silene species sampled here have sex chromosomes. Merging autosomal and sex-linked markers can be problematic in population genetics studies as they have different effective population sizes and selection dynamics [START_REF] Gautier | Using genotyping data to assign markers to their chromosome type and to infer the sex of individuals: a Bayesian model-based classifier[END_REF]. Therefore, only autosomal ORFs -as predicted by the SEX-DETector method (Muyle et al. 2016) -were kept for downstream analyses (see Materials and methods). A phylogenetic tree was built in order to check for congruence with published phylogenies of the Silene genus [START_REF] Jenkins | A phylogenetic comparative study of preadaptation for invasiveness in the genus Silene (Caryophyllaceae)[END_REF][START_REF] Rautenberg | Geographic and phylogenetic patterns in Silene section Melandrium (Caryophyllaceae) as inferred from chloroplast and nuclear DNA sequences[END_REF]. Orthologous ORFs were aligned among Silene species and Dianthus chinensis, which was used to root the tree. A total of 55,337 variable positions, spread over 1,294 contigs, were used to build the tree (see Materials and methods). All samples from a given species grouped together as expected (Fig. 1). The two lineages of S. nutans formed two clearly distinct groups that should be treated as distinct species as shown in previous studies [START_REF] Martin | Phylogeographic pattern of range expansion provides evidence for cryptic species lineages in Silene nutans in Western Europe[END_REF][START_REF] Martin | Lineages of Silene nutans developed rapid, strong, asymmetric postzygotic reproductive isolation in allopatry[END_REF]. S. viscosa samples also formed two highly divergent lineages, even more divergent than the two lineages of S. nutans: a western lineage (S. viscosa W), which includes samples from Sweden and the Czech Republic, and an eastern lineage (S. viscosa E), which includes samples from Bulgaria, Russia and Kyrgyzstan. We treated these lineages as two differentiated taxa and analyzed them separately hereafter. Unlike in previous work [START_REF] Slancarova | Evolution of sex determination systems with Population genomic consequences of separate sexes in plants heterogametic males and females in silene[END_REF], S. vulgaris did not group with S. viscosa and the section Melandrium, but instead could not be reliably grouped with either subgenus (Fig. 1).

A total of 1,088,269 variable nucleotides were called (Table 1) using an RNA-seq pipeline (see Materials and methods). Mean heterozygosity levels and Weir and Cockerham F IT statistics were computed for each species. F IT indicates deviations from Hardy-Weinberg equilibrium caused by processes such as selfing or population structure. F IT values were low to moderate (below 0.2) for most sampled species. This indicates that neither population structure nor selfing are creating strong among-species differences in genetic diversity, as already observed [START_REF] Moyle | Correlates of genetic differentiation and isolation by distance in 17 congeneric Silene species[END_REF]). However, S. viscosa E exhibited a high F IT , a low heterozygosity level and low genetic diversity, which could be due to recent demographic shifts or a high selfing rate (Fig. 2).

Proxies for species abundance were calculated based on the number of occurrences of each species retrieved from the Global Biodiversity Information Facility (GBIF) database. The geographic range of each species was extracted from to the Atlas Florae Europaeae (see Materials and methods, Fig. S6, Fig. 3). Dioecious and non-dioecious species had similar geographic ranges (622.6 versus 784 Atlas dots on average, Wilcoxon one-sided rank sum test p-value=0.52) and GBIF species abundance (30,777 versus 23,968.5 on average, Wilcoxon one-sided rank sum test p-value=0.43). This is not consistent with the dead-end hypothesis which predicts smaller geographic distributions and population sizes for dioecious species. However, some Silene dioecious species have very limited geographic ranges, making the study of their genetic diversity informative to test whether they are genetically endangered.

Genetic diversity is higher in dioecious Silene species, regardless of census population size differences

In order to test whether dioecious and non-dioecious Silene species differ in their effective population size (N e ), synonymous genetic diversity (π S ) was used as a proxy for N e with π S =4N e μ (the mutation rate μ was assumed constant among species). N e is a direct measure of the long term mean intensity of drift in populations. Under the evolutionary dead-end hypothesis, dioecious species are expected to have lower N e and therefore lower values of π S compared to non-dioecious species (see Introduction). Table 1 and Fig. 2 showed that dioecious species have higher π S values than non-dioecious species (Wilcoxon one-sided rank sum test W = 35, p-value = 0.026), which does not support the dead-end predictions. However, the 13 species studied here have very different census population sizes (N) with some species being cosmopolitan (e.g. S. latifolia, S. vulgaris) while others are rare or even endemic (e.g. S. diclinis). Genetic diversity can be affected by N, although the relationship can be complex such that N and N e are quite different [START_REF] Kalinowski | Relationship of Effective to Census Size in Fluctuating Populations[END_REF]. In order to make a fair comparison of dioecious and non-dioecious species, we accounted for the effect of N when looking at the effect of sexual systems on genetic diversity. Two different proxies for N were used:

(i) the species abundance in the GBIF dataset, and (ii) the geographic range of each species according to the Atlas Florae Europaeae. We used a linear model to jointly assess the effects of N and sexual systems on π S . GBIF data and Atlas geographic range (estimates of N) were significantly positively associated with π S (pvalues = 0.030 and 0.012, t-values = 2.529 and 3.054, R 2 = 0.390 and 0.483 respectively, and see Fig. 3). This shows that genetic diversity is strongly influenced by N in our Silene sample set. This also suggests that, in the case of the Silene genus, GBIF abundance and Atlas geographic ranges are sensible proxies for N. After controlling for the effect of N on genetic diversity, dioecious species Population genomic consequences of separate sexes in plants had a significantly higher π S than non-dioecious species (linear model pairwise t-test p-values = 0.142 and 0.048, t-values = 1.595 and 2.247 using GBIF data and Atlas geographic range estimates respectively). Adding F IT to the linear model correlating genetic diversity and N did not significantly improve our linear model (anova p=0.102 and p=0.097 for GBIF and Atlas geographic range estimates respectively, see Materials and Methods for details). This suggests that the increase in genetic diversity observed under dioecy is not attributable to inbreeding evasion. These effects remained after correcting for phylogenetic relationships among species using generalized least square models in the R package ape [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF]: N was significantly positively correlated to π S (p-values = 0.012 and 0.03, t-values = 3.05 and 2.53 for Atlas and GBIF respectively) and dioecious species had a significantly higher π S compared to nondioecious species when using Atlas geographic range estimates (p-value = 0.048, t-value = 2.25), but not when using GBIF estimates (p-value = 0.14, t-value = 1.6). In this analysis, different sets of genes were used in different species. As this might affect the results, the analysis was repeated on two smaller sets of genes, one shared by every species in the subgenus Behenantha and the second shared by every species in the subgenus Silene. Similar results were found using these shared gene subsets, i.e. GBIF data and Atlas geographic range have a significant positive effect on π S (p-values = 0.021 and 0.012, t-values = 2.725 and 3.074, R 2 = 0.429 and 0.488, respectively). And dioecious species had higher π S compared to other species after correcting for N (linear model pairwise t-test p-values = 0.127 and 0.045, t-values = 1.662 and 2.283 using GBIF data and Atlas geographic range respectively, see Fig. S1).

Selection is more efficient in dioecious Silene species, regardless of differences in N e

We then tested whether selection efficacy differs in dioecious and non-dioecious Silene species.

Selection should be more efficient in species with larger N e . The ratio of non-synonymous to synonymous substitution rate, d N /d S , could not be used here because the studied species are closely related and some shifts in sexual systems might have occurred too recently to allow for fixed differences to evolve [START_REF] Desfeux | Evolution of reproductive systems in the genus Silene[END_REF][START_REF] Bernasconi | Silene as a model system in ecology and evolution[END_REF]. Rather, the ratio of nonsynonymous over synonymous polymorphism, π N /π S , was used, because it is a more appropriate statistic to study the efficacy of selection on a short evolutionary timescale [START_REF] Mugal | Why time matters: codon evolution and the temporal dynamics of dN/dS[END_REF]. One expects increased π N /π S values when purifying selection is reduced and hence the dead-end hypothesis would predict elevated π N /π S in dioecious species. However, we found no statistical association between sexual systems and π N /π S (Wilcoxon rank sum test p=0.295, Table 1), perhaps reflecting the opposing effects that purifying and positive selection have on π N /π S , or perhaps reflecting the fact that breeding systems do not have an effect on π N /π S .

We thus turned to a more sophisticated approach to quantify selection efficacy and used a modified version of the McDonald-Kreitman framework, Implemented in the method polyDFE [START_REF] Tataru | Inference of Distribution of Fitness Effects and Proportion of Adaptive Substitutions from Polymorphism Data[END_REF][START_REF] Tataru | polyDFE: Inferring the Distribution of Fitness Effects and Properties of Beneficial Mutations from Polymorphism Data[END_REF]. This method estimates the distribution of fitness effects (DFE) of both deleterious and beneficial mutations using polymorphism data only, thus avoiding the problem of shifts in population sizes. Silene species underwent a post-glacial age population expansion in Europe and for this reason it would be biased to use methods that also rely on divergence [START_REF] Keightley | Joint inference of the distribution of fitness effects of deleterious mutations and population demography based on nucleotide polymorphism frequencies[END_REF][START_REF] Galtier | Adaptive Protein Evolution in Animals and the Effective Population Size Hypothesis[END_REF]. Indeed, other methods to estimate the DFE use both polymorphism and divergence data and usually return estimates of proportion of adaptive substitution that are assuming that population size remained constant since divergence with the outgroup. S. viscosa E was not included because its diversity was too low due to selfing and the polyDFE model did not apply. PolyDFE was used to estimate the DFE (Fig. S2) for each species with a gamma distribution for deleterious mutations and exponential distribution for beneficial mutations (see estimated parameter ranges in Fig. S3). From this, the predicted proportion of adaptive substitutions (α), the predicted rate of adaptive substitutions (ω A ), and the predicted rate of deleterious substitutions (ω NA ) were estimated. Under the dead-end hypothesis, dioecious species are expected to have a higher ω NA , a lower α and a lower ω A because of a faster accumulation of deleterious mutations and a lower rate of adaptation due to genetic drift, when compared to nondioecious species. Because N e is known to influence the efficacy of selection, we correct for differences in N e (using π S as a proxy) between species in order to make a fair comparison of dioecious and non-dioecious species. A linear model was used to jointly assess the effects of N e (using π S as a proxy) and of the sexual system on α. π S was significantly positively correlated to α (p-value = 0.003, t-value = 4.07, Fig. 4). The same result has previously been observed in animal species [START_REF] Galtier | Adaptive Protein Evolution in Animals and the Effective Population Size Hypothesis[END_REF] where the correlation between α and π S was solely due to the link between π S and ω NA . Because α = ω A / (ω A + ω NA ), when selection is more efficient (higher N e and higher π S ), the rate of deleterious substitutions ω NA is reduced, which in turn increases α. However, in the present dataset, the correlation between α and π S was both due to higher adaptive evolutionary rates and lower deleterious evolutionary rates when effective population size is high. Indeed, π S was significantly positively linked to ω A (linear model p-value = 0.011, t-value = 3.19) and significantly negatively correlated to ω NA (linear model p-value = 0.00017, t-value = -6.15, see Fig. 4). These effects remained after correcting for phylogenetic relationships among species using gls models in the R package ape (Paradis and Schliep 2019): π S was significantly positively correlated to α and ω A (p-values = 0.026 and 0.0054, t-values = 4.11 and 3.64 for α and ω A respectively) and π S was negatively correlated to ω NA (p-value = 0.0015 and t-value = -4.49).

Estimates of α, ω A and ω NA were then compared between dioecious and non-dioecious species, after controlling for the effect of N e (using π S as a proxy). Non-dioecious species had a significantly lower α and ω A and a significantly higher ω NA compared to dioecious species (linear model pairwise t-test p-values = 0.015, 0.00031 and 0.019, t-values = -2.99, -5.65 and 2.87 for α, ω A and ω NA respectively). These effects remained significant after correcting for phylogenetic relationships among species using gls models in the R package ape (p-values = 0.032, 0.0007 and 0.058, t-values = -2.53, -5.03 and 2.17 for α, ω A and ω NA respectively). Dioecious species thus experience fewer fixations of slightly deleterious mutations due to genetic drift and more adaptive evolution, which is opposite to what we expect under the dead-end hypothesis. The analyses were repeated using orthologous genes shared within the subgenus Behenantha or within the subgenus Silene as different sets of genes in different species could affect the results. The conclusions remained unchanged; differences were still in the same direction although statistical power was reduced in some cases (see Fig. S4 legend for details).

Discussion

In the set of Silene species studied here, dioecy is associated with higher genetic diversity, more efficient purifying selection and increased adaptive evolution. These results contradict the dead-end hypothesis and are consistent with macro-evolutionary analysis that rejected the dead-end from different lines of evidence. High levels of genetic diversity in dioecious species had previously been observed in Silene dioica [START_REF] Giles | Genetic differentiation in Silene dioica metapopulations: estimation of spatiotemporal effects in a successional plant species[END_REF], Silene chlorantha [START_REF] Lauterbach | Genetic population structure, fitness variation and the importance of population history in remnant populations of the endangered plant Silene chlorantha (Willd.) Ehrh. (Caryophyllaceae)[END_REF],

Population genomic consequences of separate sexes in plants

Silene otites [START_REF] Lauterbach | Population genetics and fitness in fragmented populations of the dioecious and endangered Silene otites (Caryophyllaceae)[END_REF], Ficus [START_REF] Nazareno | Dioecy, more than monoecy, affects plant spatial genetic structure: the case study of Ficus[END_REF], Antennaria dioica [START_REF] Rosche | Sex ratio rather than population size affects genetic diversity in Antennaria dioica[END_REF]) and a comparison of 12 monoecious and 11 dioecious plants [START_REF] Nazareno | Dioecy, more than monoecy, affects plant spatial genetic structure: the case study of Ficus[END_REF]. Some of these studies have shown that the sex ratio strongly influences genetic diversity through its influence on the outcrossing rate. We do not have access to sex ratios of the species sampled here but their study would be insightful to better understand patterns of diversity in Silene species in future. Our results clearly show that dioecious Silene species are not at risk of extinction given their high genetic diversity. Observed rates of adaptive non-synonymous substitution patterns suggest that dioecious species undergo more adaptation than non-dioecious species, however, the underlying causes of this pattern remain to be determined.

The Silene genus is a good model system to test for the dead-end hypothesis because dioecious Silene species do not have any of the ecological traits thought to possibly compensate for the handicaps of dioecy (see Introduction). In other genera in which dioecy is ancient, it is likely that enough time has passed for species to evolve compensating mechanisms to counteract the handicaps of dioecy, if these truly exist. This is less likely to be an issue in the Silene genus where dioecy has evolved fairly recently: ~11 Mya old in the Melandrium section [START_REF] Krasovec | The Mutation Rate and the Age of the Sex Chromosomes in Silene latifolia[END_REF] and even more recently in section Otites [START_REF] Slancarova | Evolution of sex determination systems with Population genomic consequences of separate sexes in plants heterogametic males and females in silene[END_REF].

Our approach, however, also has its limitations. First, our set of species represents only two independent events of dioecy evolution and a small sample of the whole Silene genus (~700 species). However, we did not find a significant phylogenetic signal on genetic diversity in our species set (Pagel's λ=0.848, p=0.133), suggesting that each species can be considered as independent even though some species share the same transition event to dioecy. Nonetheless, other genera should be studied in the future to assess whether our findings can be generalized to other angiosperm species, especially to groups in which dioecy evolved through different routes (see

Introduction).

A second limitation is that sampling was limited by the number of Silene species for which the sexual system is well described and for which the geographic distribution is similar among species, especially for hermaphrodites (see Materials and Methods for detailed criteria used to select species). None of the Silene hermaphrodite species have geographic distributions as wide as the widest dioecious species range. However, dioecious species that have geographic ranges comparable to those of the hermaphroditic species we sampled still have higher diversity and selection efficacy compared to hermaphrodites (Fig. 3 and4), which is conservatively in disagreement with the dead-end hypothesis. For two species out of thirteen, namely S. otites (dioecious) and S. paradoxa (non-dioecious), we were not able to sample individuals across the whole geographic range of the species and their genetic diversity might be underestimated compared to others.

Population genomic consequences of separate sexes in plants

A third limitation comes from the comparison of dioecious and non dioecious species which implies comparing species with sex chromosomes to species with only autosomes. Sex-linked genes, which are located on the non-recombining region of sex chromosomes, are known to suffer from lower N e , which could bias comparisons of dioecious to non dioecious species. Because sexlinked genes in Silene represent <10% of genes (discussed in [START_REF] Muyle | Rapid de novo evolution of X chromosome dosage compensation in Silene latifolia, a plant with young sex chromosomes[END_REF], we chose to subset our study to autosomal genes only. Our results therefore overlook the potential negative effects of sex chromosomes on selection efficacy in dioecious species.

A fourth limitation of our study is that, if the extinction dynamics of dioecious species is extremely fast, there might not be enough time for deleterious genetic footprints to accumulate in dioecious species before they go extinct. The micro-evolutionary approach taken here can only detect handicaps of dioecy that do not lead to immediate extinction. Such a study, as others testing evolutionary dead-end hypotheses, suffers from a "survivor bias" as only non-extinct species can be analyzed. However, some of the dioecious species studied here have small geographic distributions and one is even endemic. These species should be strongly affected by the handicaps of dioecy, if there are any. But even these species exhibited high levels of diversity and selection efficacy.

In the Silene genus, all hermaphrodites are self-compatible and therefore able to selffertilize. Outcrossing promotes genetic diversity by increasing effective population size and minimizing the effects of linkage [START_REF] Nordborg | Linkage disequilibrium, gene trees and selfing: an ancestral recombination graph with partial self-fertilization[END_REF][START_REF] Glémin | Genome Evolution in Outcrossing vs. Selfing vs. Asexual Species[END_REF]). Dioecy could be beneficial in Silene by preventing the negative effects of selfing that occur in hermaphroditic and gynodioecious species [START_REF] Charlesworth | A Model for the Evolution of Dioecy and Gynodioecy[END_REF]. One selfing hermaphrodite subspecies (S. viscosa E) could not be included in the analysis of selection efficacy because its genetic diversity was too low. However, two hermaphrodite species that are not heavy selfers (Figure 2) were included in our comparison of dioecious and non-dioecious species. However, the fact that significant differences in terms of selection efficacy between dioecious and non-dioecious species remain, even without this selfing species, highlights that the advantages of dioecy that we observe do not result from comparing selfing and outcrossing species. What is more, the increased genetic diversity observed in dioecious species could not be explained by differences in inbreeding levels in our dataset.

However, the main difference between dioecious species and hermaphrodites is the ability to occasionally self-fertilize, which can strongly affect the demography of a species. Indeed, in hermaphroditic self-compatible species, a single individual can start a new population [START_REF] Schoen | Intraspecific variation in population gene diversity and effective population size correlates with the mating system in plants[END_REF]. This reproductive assurance provides an ecological advantage to hermaphroditic selfcompatible species over dioecious species. However, the other side of the coin is that hermaphroditic self-compatible species can experience recurrent bottlenecks, even if the mean selfing rate is close to zero. This reduces genetic diversity and N e and makes selection less efficient.

Population genomic consequences of separate sexes in plants

On the other hand, dioecious populations require both males and females to be present in sufficient density in a location in order to reproduce. Therefore, dioecious population sizes cannot decline below a given threshold at risk of declining to extinction, a phenomenon called the Allee effect [START_REF] Courchamp | Allee Effects in Ecology and Conservation[END_REF]. A counter-intuitive prediction is thus that, conditioning on survival, the Allee effect promotes genetic diversity in dioecious populations, especially during colonization [START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF][START_REF] Romiguier | Comparative population genomics in animals uncovers the determinants of genetic diversity[END_REF]. Silene species are known to have colonized Europe northward after the last ice age and have experienced population expansion. The Allee effect could have played a strong role during this expansion with self-compatible hermaphrodites colonizing new habitats fast with as few as a single individual and undergoing strong bottlenecks and diversity loss. Dioecious species that survived, on the other hand, would have expended slowly, maintaining sufficiently dense populations throughout the expansion, resulting in the maintenance of genetic diversity. In established dioecious species, the Allee effect could provide a long-term safeguard against genetic depletion.

Our results suggest that established dioecious Silene species are genetically healthy. Our work revisits the fascinating question of why dioecy is rare in angiosperms, since the dead-end hypothesis does not seem to apply, at least to our sample of species. However, our dataset might suffer from a "survivor bias" and we cannot exclude that potential handicaps of dioecy could lead to very early extinction during the emergence of new dioecious populations from hermaphrodite lineages. At a macro-evolutionary scale, this rapid dead-end would lead to an observed low rate of transition towards dioecy. However, current estimates do not show low transition rates to dioecy in angiosperms [START_REF] Goldberg | Macroevolutionary synthesis of flowering plant sexual systems[END_REF], which goes against the possibility of an early dead-end of dioecy.

An alternative hypothesis to explain the rarity of dioecy is that rates of transition to dioecy could be low in angiosperms and/or dioecy could be more evolutionary labile than previously thought and easily lost and reversed to non-dioecious states depending on environmental conditions [START_REF] Goldberg | Macroevolutionary synthesis of flowering plant sexual systems[END_REF][START_REF] Kafer | On the rarity of dioecy in flowering plants[END_REF]. The combination of a low transition rate to dioecy and a high reversibility to non-dioecy together would lead to a low probability to observe dioecious angiosperms in nature, as they would only be transiently existing. This would make dioecy rare without the need for it to be an evolutionary dead-end. Larger and fully resolved angiosperm phylogenies with annotated breeding systems will be required to precisely estimate more precisely transition and reversion rates to and from dioecy.

Materials and methods

Sampling and sequencing

Population genomic consequences of separate sexes in plants

We investigated patterns of polymorphism and divergence of three sexual systems (hermaphroditism, gynodioecy and dioecy) in two subgenera of the Silene genus, subgenus Behenantha and subgenus Silene (sample sizes and sexes are indicated in Fig. 1). For subgenus Behenantha, we sampled five dioecious species (S. latifolia, S. diclinis, S. dioica, S. marizii and S. heuffelli), one gynodioecious species (S. vulgaris), and one hermaphroditic species (S. viscosa). For the subgenus Silene, we sampled two dioecious species (S. otites and S. pseudotites), two gynodioecious sub-species (S. nutans W -Western lineage, and S. nutans E -Eastern lineage, as defined in [START_REF] Martin | Phylogeographic pattern of range expansion provides evidence for cryptic species lineages in Silene nutans in Western Europe[END_REF]) and one hermaphroditic species (S. paradoxa). Two samples of Dianthus chinensis were added to root the phylogeny. Some of the individuals were already used in previous studies (S1).

The two known dioecious sections in the genus Silene were sampled: section Melandrium (all five species of the section were selected) and section Otites (2 closely related dioecious species were selected). The selection of the non-dioecious species was done using different criteria. First, close relatives of the two dioecious sections were selected based on available phylogenies [START_REF] Desfeux | Evolution of reproductive systems in the genus Silene[END_REF][START_REF] Jenkins | A phylogenetic comparative study of preadaptation for invasiveness in the genus Silene (Caryophyllaceae)[END_REF][START_REF] Slancarova | Evolution of sex determination systems with Population genomic consequences of separate sexes in plants heterogametic males and females in silene[END_REF]. We focused on those for which the sexual system is well described. The Silene genus is a large genus comprising ~700 species but only ~60 have a well-described sexual system (dioecious species included, [START_REF] Jürgens | Pollen grain numbers, ovule numbers and pollen-ovule ratios in Caryophylloideae: correlation with breeding system, pollination, life form, style number, and sexual system | SpringerLink[END_REF]. Highly selfing species were avoided by focusing on hermaphroditic species that are perennial, insect pollinated (or at least not autogamous/cleistogamous) and with a high pollen/ovule ratio [START_REF] Jürgens | Pollen grain numbers, ovule numbers and pollen-ovule ratios in Caryophylloideae: correlation with breeding system, pollination, life form, style number, and sexual system | SpringerLink[END_REF]. These criteria left us with only a handful of hermaphroditic species, of which S. viscosa and S. paradoxa were selected for seed availability.

Seeds were sampled across the geographic distribution of each species (Fig. S5) and were sown and grown in controlled greenhouse conditions. In case plants did not flower, young leaves were sampled instead (details on individual names, sex, origin and sampled tissue can be found in Table S1). Total RNA was extracted through the Spectrum Plant Total RNA kit (Sigma, Inc., USA) following the manufacturer's protocol and treated with DNAse. Libraries were prepared with the TruSeq RNA sample Preparation v2 kit (Illumina Inc., USA). Each 2 nM cDNA library was sequenced using a paired-end protocol on an Illumina HiSeq2000 sequencer, producing 100 bp reads. Demultiplexing was performed using CASAVA 1.8.1 (Illumina) to produce paired-sequence files containing reads for each sample in Illumina FASTQ format. RNA extraction and sequencing were done by the sequencing platform in the AGAP laboratory, Montpellier, France (http://umragap.cirad.fr/). The Melandrium and Dianthus samples were sequenced as described in [START_REF] Zemp | Evolution of sexbiased gene expression in a dioecious plant[END_REF][START_REF] Zemp | Evolution of sexbiased gene expression in a dioecious plant[END_REF] at the Quantitative Genomics Facility (QGF, ETH Zürich, Switzerland) and by GATC Biotech AG (Konstanz, Germany).

Population genomic consequences of separate sexes in plants read mapping quality threshold of 10 and minimum base quality threshold of 20. The output obtained was two unphased sequences for each individual and each ORF.

Orthologous ORFs were identified using OrthoMCL [START_REF] Li | OrthoMCL: identification of ortholog groups for eukaryotic genomes[END_REF]) Basic Protocol 2 with default parameters. Best reciprocal blast hits were considered orthologs among species and orthologous ORFs were aligned for all individuals using MACSE [START_REF] Ranwez | MACSE: Multiple Alignment of Coding SEquences accounting for frameshifts and stop codons[END_REF]. Species from the subgenus Behenantha were aligned separately with S. paradoxa, considered as an outgroup for this section and species from the subgenus Silene were aligned with S. viscosa considered as an outgroup for this section.

Phylogenetic reconstruction and phylogenetic signal

Orthologous ORFs among species were obtained using best reciprocal blast hits among the 14 species (all species from the Silene genus and Dianthus chinensis). For each orthologous ORF, the nine species were added one by one into the alignment using MACSE [START_REF] Ranwez | MACSE: Multiple Alignment of Coding SEquences accounting for frameshifts and stop codons[END_REF]). The data were filtered so that a maximum of 20% missing data per site and per species was allowed, for a total of 55,337 variable nucleotides across 1,294 ORFs. We used the maximum number of ORFs available and did not restrict our analysis to autosomal ORFs identified in dioecious species. The phylogenetic tree was reconstructed using the GTRGAMMA model in RaxML [START_REF] Stamatakis | RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models[END_REF].

Bootstrap support values were generated using 100 replications.

Phylogenetic signal and its significance were computed using the R function phylosig in the phytools R package [START_REF] Revell | phytools: an R package for phylogenetic comparative biology (and other things)[END_REF].

Data filtering

Analyses were carried out on autosomal genes only, in order to avoid potential confounding effects of sex chromosomes and sexual systems on estimates of selection efficiency. Autosomal genes were predicted for each dioecious species separately by SEX-DETector (Muyle et al. 2016) thanks to the study of a cross (parents and progeny) sequenced by RNA-seq [START_REF] Martin | Evolution of young sex chromosomes in two dioecious sister plant species with distinct sex determination systems[END_REF] and Muyle et al. in prep). 11,951 contigs were inferred as autosomal in S. latifolia, 7,695 for S. diclinis, 11,675 for S. heuffelli, 8,997 for S. marizii, 10,698 for S. dioica, 2,035 in S. otites and 6,234 in S. pseudotites.

The data were further filtered in order to remove ORFs with frameshifts. Only biallelic sites were kept for analyses (sites with two or less alleles). Only ORFs longer than 100 codons and with a maximum of 20% of missing individuals of the focal species in the alignment were retained, after removing incomplete codons and gaps. We obtained final datasets ranging from 893 retained ORFs for S. otites up to 6,413 retained ORFs for S. viscosa E (see Table 1). Finally, to verify that observed trends are not due to differences in studied ORFs among species, we performed the same

Fig. 1 .

 1 Fig. 1. Sampled species and their phylogenetic relationships. Dioecious species are shown in red, gynodioecious species in blue and hermaphroditic species in magenta. Sample sizes and sex are indicated. Bootstrap support values were generated using 100 replicates. Pictures of S. latifolia, S. dioica, S. nutans W, S. nutans E and S. otites are from Maarten Strack van Schijndel. Pictures of S. marizii, S. vulgaris, S. paradoxa and S. viscosa are from Oxelman et al (2013). Picture of S. heuffelii is from Dr Daniel L. Nickrent.

Fig. 2 .

 2 Fig. 2. Levels of individual heterozygosity (proportion of heterozygote positions), Weir and Cockerham F IT statistic and π S within each species. Dots show the mean value and vertical lines mark the 95% confidence interval (obtained by 1,000 bootstraps on contigs). Dioecious species are shown in red, gynodioecious species are shown in blue and hermaphroditic species are shown in magenta.

Fig. 3 .

 3 Fig. 3. Relationship between synonymous polymorphism (π S ) and (A) census population size estimated from the Global Biodiversity Information Facility (GBIF) and (B) geographic range estimated from the Atlas Florae Europaeae. For each species, dots represent the mean π S and vertical lines mark the 95% confidence interval (obtained by 1,000 bootstraps on contigs). Census population sizes estimated from GBIF refer to the number of occurrences, and geographic distributions from the Atlas Florae Europaeae refer to the number of dots (native occurrences) counted on the Atlas maps (see Materials and methods). Dioecious species are shown in red, gynodioecious species are shown in blue and hermaphroditic species are shown in magenta.

Fig. 4 .

 4 Fig. 4. Relationship between synonymous polymorphisms (π S ) and (A) the proportion of adaptive amino-acid substitutions, α, (B) the rate of adaptive non-synonymous substitutions (ω A ) (C) the rate of non-adaptive non-synonymous substitutions (ω NA ). Dots show the mean value, horizontal lines mark the 95% π S confidence interval (obtained by 1,000 bootstraps on contigs) and vertical lines show the y axis 95% confidence interval (obtained on 500 bootstraps on the SFS SNPs). Dioecious species are shown in red, gynodioecious species are shown in blue and hermaphroditic species are shown in magenta.

Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msaa229/5905496 by guest on 25 November 2020

Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msaa229/5905496 by guest on 25 November 2020

Downloaded from https://academic.oup.com/mbe/advance-article/doi/10.1093/molbev/msaa229/5905496 by guest on 25 November 2020

Acknowledgements

We would like to thank Jos Kafer for his help with fieldwork and a number of colleagues for sending us seeds: J. J. Wieringa, Deborah Charlesworth, Bohuslav Janousek, Tatiana Giraud, Honor C. Prentice, Salvatore Cozzolino, Delphine Griveau, Henk Schat, Cristina Gonnelli, Alessia Guggisberger, Sophie Karrenberg, Adrien Favre, Maria Domenica Moccia, Aria Minder, Carmen Rothenbuehler, Massimiliano Gnesotto and Kerstin Kustas. AM, RT and GABM thank Elise Lacroix technician of the Univ. Lyon 1 greenhouse and HM, SGa and PT thank Eric Schmitt for taking care of the Silene plants. We are grateful to Sylvain Santoni and the AGAP platform staff for the RNA extraction and RNA-seq. The numerical results presented in this paper were carried out using the computing cluster of Laboratoire Biométrie et Biologie-Pôle Rhone-Alpin de Bioinformatique (LBBE-PRABI) at Univ. Lyon 1 and the HPC service of the Centre de Ressources Informatiques (CRI) of Univ. Lille 1, for which we thank Stéphane Delmotte and Bruno Spataro.

We thank Nicolas Galtier for help with his program sfs. HM, SGa and PT are grateful to Jonathan Aceituno for technical support in R and to Céline Poux for discussion on phylogenetic inference.

We are thankful to Brandon Gaut for commenting on the manuscript. This work was supported by the Agence Nationale de la Recherche (ANR-11-BSV7-013-03, TRANS) to SGl, GM, and PT, the Swiss National Science Foundation (31003A-182675), and a PhD fellowship from the French Research Ministry to HM.

Data availability

All the RNA-seq reads generated for this study can be downloaded from the European Nucleotide Archive Accession PRJEB39526.

Population genomic consequences of separate sexes in plants

Reference transcriptome assemblies

Reference transcriptomes were assembled de novo for all species (except for S. marizii, S. heuffelli, S. diclinis and S. dioica). As assemblies were carried in different labs, they differ slightly for the two Silene taxa. For the subgenus Behenantha, reads from all individuals were pooled for each species. Steps in the assembly are represented in Fig. S3. 100% identical reads, assumed to be PCR duplicates, were filtered using the ConDeTri trimming software [START_REF] Smeds | ConDeTri -A Content Dependent Read Trimmer for Illumina Data[END_REF].

Illumina reads were then filtered for sequencing adapters and low read quality using the ea-utils FASTQ processing utilities [START_REF] Aronesty | Command-line tools for processing biological sequencing data[END_REF]. Filtered reads were assembled using trinity with default settings [START_REF] Haas | De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[END_REF]. Poly-A tails were removed using PRINSEQ (Schmieder and Edwards 2011) with parameters -trim_tail_left 5 -trim_tail_right 5. rRNA-like sequences were removed using riboPicker version 0.4.3 [START_REF] Schmieder | Identification and removal of ribosomal RNA sequences from metatranscriptomes[END_REF] with parameters -i 90 -c 50 -l 50 and the following databases: SILVA Large subunit reference database, SILVA Small subunit reference database, the GreenGenes database and the Rfam database. Obtained transcripts were further assembled, inside of trinity components, using the program Cap3 [START_REF] Huang | CAP3: A DNA Sequence Assembly Program[END_REF] with parameter -p 90 and custom perl scripts.

Raw data of the five taxa of subgenus Silene were filtered for sequencing adapters using Cutadapt [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF], low read quality and poly-A tails with PRINSEQ [START_REF] Schmieder | Quality control and preprocessing of metagenomic datasets[END_REF]. For each focal species of the subgenus, reference transcriptomes were assembled de novo by combining filtered reads from two individuals from the same population (the two individuals with the highest number of reads) using default settings of Trinity [START_REF] Haas | De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[END_REF]) and an additional Cap3 (Huang and Madan 1999: 3) run. Finally, in both subgenera, coding sequences were predicted using TransDecoder [START_REF] Haas | De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[END_REF] with PFAM domain searches as ORF retention criteria (see website: http://transdecoder.github.io). Assembly statistics are presented in Table S2.

Mapping, genotyping and alignment

Reads from each species were mapped onto their respective reference transcriptome using BWA (Li and Durbin 2009) version 0.7.12 with parameter -n 5 for the subgenus Behenantha and Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] for the subgenus Silene. Mapping statistics can be found in Table S3. For the Melandrium section (S. latifolia, S. marizii, S. heuffelli, S. diclinis and S. dioica), all reads were mapped to the S. latifolia reference transcriptome. SAM files were compressed and sorted using SAMtools (Li et al. 2009). Individuals were genotyped using reads2snp_2.0 [START_REF] Tsagkogeorga | The population genomics of a fast evolver: high levels of diversity, functional constraint, and molecular adaptation in the tunicate Ciona intestinalis[END_REF][START_REF] Gayral | Reference-free population genomics from nextgeneration transcriptome data and the vertebrate-invertebrate gap[END_REF]) with parameters -par 1 -min 3 (10 for the subgenus Silene) -aeb -bqt 20 -rqt 10 that is to say: stringent filtering of paralogous positions, minimum coverage of 3 (or 10) reads for calling a genotype, accounting for allelic expression bias, minimum Population genomic consequences of separate sexes in plants analyses on a smaller set of orthologous ORFs present in all species within a subgenus, so that parameters were estimated on the same set of ORFs. For these tests, 188 orthologous ORFs were analysed for the subgenus Silene, and 298 for the subgenus Behenantha.

Polymorphism, heterozygosity and F IT

Population statistics were computed on the filtered dataset for each ORF and averaged using the PopPhyl tool dNdSpiNpiS_1.0 [START_REF] Gayral | Reference-free population genomics from nextgeneration transcriptome data and the vertebrate-invertebrate gap[END_REF]) with parameter kappa=2. The following statistics were computed for each contig: the mean fixation index F IT [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF], per-individual heterozygosity (proportion of heterozygous positions), per-site synonymous polymorphism (π S ), per-site non-synonymous polymorphism (π N ) and the ratio of non-synonymous to synonymous polymorphism (π N /π S ). Note that we could not compute F IS or F ST as we do not have the required 2 individuals / population. Confidence intervals were obtained by 1,000 bootstraps on contigs.

Estimation of geographic distribution and population size

Geographic distributions of the Silene species were estimated by counting the black dots (native occurrence) on the Atlas maps by [START_REF] Jalas | Distribution of Vascular Plants in Europe. 7. Caryophyllaceae (Silenoideae)[END_REF]. This was partly achieved using the Dot Count program (http://reuter.mit.edu/software/dotcount). Since no map was available for S. pseudotites, dots were counted for Italy and adjacent regions of France and the Balkans using the S. vulgaris map (the Atlas indicates these are the regions where this plant is abundant).

Census population sizes of the Silene species were estimated using raw species occurrence numbers from GBIF (http://www.gbif.org 17th January 2017 Download of S. dioica: http://doi.org/10.15468/ Population genomic consequences of separate sexes in plants ranges and GBIF species abundance were compared among dioecious and non dioecious species using a one-sided Wilcoxon rank sum test in R (R Core Team 2016). Geographic range and GBIF data were correlated to π S and sexual systems or F IT using R linear model function lm as follows: lm(piS ~ Atlas + sexual system) lm(piS ~ GBIF + sexual system) lm(piS ~ Atlas + F IT ) lm(piS ~ GBIF + F IT )

The significancy of the effect of sexual system and F IT were tested thanks to a comparison to a nested model using the anova function in R. The nested models are:

The same models were used and corrected for phylogenetic relationships among species using gls models with Martins and Hansen's correlation structure (corMartins) in the R package ape [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF].

Selection efficacy

The unfolded synonymous and non-synonymous site frequency spectra (SFS, the distribution of derived allele counts across SNPs) were computed for each species with the sfs program provided by Nicolas Galtier [START_REF] Galtier | Adaptive Protein Evolution in Animals and the Effective Population Size Hypothesis[END_REF]. These SFS were used to estimate α (the proportion of adaptive amino-acid substitutions) in each species. We estimated α according to a modified version of the Eyre-Walker and Keightley method (Keightley and Eyre-Walker 2007), which relies on polymorphism data only to estimate α [START_REF] Tataru | Inference of Distribution of Fitness Effects and Proportion of Adaptive Substitutions from Polymorphism Data[END_REF]. polyDFE (https://github.com/paulatataru/polyDFE) was used to estimate the distribution of fitness effects of mutations (DFEM), α, the (per synonymous substitution) rate of adaptive non-synonymous substitution ω A =α d N /d S and the rate of non-adaptive non-synonymous substitutions ω NA =(1-α) d N /d S . For each species separately, four different models were fitted to the SFS:

 Model 1: full DFEM with adaptive mutations (p_b > 0), estimation of an orientation error rate for ancestral alleles in the SFS (eps_an), estimation of demography (r_i ≠ 1) for a total Population genomic consequences of separate sexes in plants of 6 parameters plus the number of demography parameters r_i which depends on the sample size of each species.

 Model 2: same as Model 1 but without ancestral allele orientation error rate (eps_an=0), for a total of 5 parameters plus the number of demography parameters r_i.

 Model 3: same as Model 1 but without adaptive mutations (p_b=0), for a total of 5 parameters plus the number of demography parameters r_i.

 Model 4: same as Model 1 but without adaptive mutations (p_b=0) and without ancestral allele orientation error rate (eps_an=0), for a total of 4 parameters plus the number of demography parameters r_i.

Each model provided estimates for α, ω A and ω NA . These values were averaged using Akaike weights as follows:

where m stands for model number m, α m stands for the estimated α using model m, and ΔAICAIC m =AIC m -AIC min with AIC the Akaike criterion:

Where n m stands for the number of parameters and L m for the likelihood of model m. Model averaged values for ω A and ω NA were obtained using a similar procedure.

For each species, 500 bootstrapped datasets were randomly generated by random sampling SNPs with replacement from the original SFS. Average values among models were computed as explained previously for each bootstrap. The distribution of values for the 500 bootstraps was used to compute approximate 95% confidence intervals for α, ω A and ω NA . For this, 2.5% of values were removed from each side of the bootstrap distribution [START_REF] Efron | An Introduction to the Bootstrap[END_REF]. α, ω A and ω NA values were analysed in a linear model using R lm function (R Core Team 2016). The model tested the effect of π S and sexual system (dioecious, gynodioecious, hermaphroditic) as fixed effect, while taking into account the variance of the estimations with weights computed as follow: lm(α ~ π S + sexual system, weights) lm(ω A ~ π S + sexual system, weights) lm(ω NA ~ π S + sexual system, weights)

The same models (without weights) were used and corrected for phylogenetic relationships among species using gls models with Martins and Hansen's correlation structure (corMartins) in the R package ape [START_REF] Paradis | ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R[END_REF].

Evolution of separate sexes in plants


Maxime Houssin
email: maxime.houssin@irit.fr

Stéphanie Combettes
email: stephanie.combettes@irit.fr

Marie-Pierre Gleizes
email: marie-pierre.gleizes@irit.fr

Berangere Lartigue
email: berangere.lartigue@univ-tlse3.fr

SANDMAN: a Self-Adapted System for Anomaly Detection in Smart Buildings Data Streams

Keywords: Self-Adaptive Multi-Agent System, Anomaly Detection, Smart Buildings

Currently, energy management within buildings is essential to mitigate climate change. To this end, buildings are increasingly equipped with sensors to assist the building manager. Yet, the heterogeneity and the large amount of generated data make this task quite difficult. The SANDMAN multi-agent system, described in this paper, aims to assist in the automatic detection, in constrained time, of several types of anomalies using raw and heterogeneous data. SANDMAN features a semisupervised learning by considering some feedback from an expert in the field. The results show that SANDMAN detects different types of anomalies, is resilient to noise and is scalable.

I. INTRODUCTION

Buildings account for more 30% of worldwide final energy consumption, this figure rising to over 40% in developed countries. However, much of this energy is wasted -30 % in the United States and Europe. Among other reasons, this waste comes from poor building management, such as the failure to quickly detect problems.

The objective of the work presented in this article is to detect anomalies in energy data issued from smart buildings to quickly solve the problems and minimise the energy waste.

With the advent of IoT (Internet of Things), the number of sensors in existing and new buildings have sharply risen due to their lower costs and the obvious benefit of using them for building management. Sensors can easily be added to buildings or replaced by others. Therefore, managing these sensors and the data they generate is a complex system for building managers.

However, an accurate control of these data is necessary for the proper management of buildings, including their energy performance. This is why a tool for automatic detection of anomalies is an important asset for building managers. This tool must process the data in constrained time to allow to act as soon as possible, which is an essential issue in the detection of anomalies. Because of the large amount of data, the spaces for searching for anomalies also become very large: this is the reason why the system must learn to raise anomalies in a semi-supervised manner thanks to the feedback of an expert of the field. The system detects and gives the found anomalies to the expert and the expert gives feedback to allow the system to learn throughout its life.

The article is composed as follows. Section I presents the characteristics required for anomaly detection systems and defines the notions of anomalies. Section II describes the main existing anomalies detection systems. Section III explains the architecture of the proposed multi-agent system, SANDMAN (Semi-Supervised ANomaly Detection with Multi-AgeNt systems). Section IV is devoted to SANDMAN's experiments and results. Finally, Section V concludes and proposes perspectives to the presented work.

A. Required Characteristics for an Anomaly Detection System

An intelligent building equipped with sensors generates a large amount of available data that must be analysed to improve energy management. For example, the SGE (Service de Gestion et d'Exploitation in charge of energy management for the campus of University of Toulouse III-P. Sabatier) manages about 6000 sensors measuring at least one value per hour. Thus, the design of a smart building anomaly detection system must ideally take into account the following characteristics:

• interaction with an expert, • constrained time detection, • detection of several types of anomalies, • use of raw data, • use of heterogeneous data, • openness, • scaling, • genericity.

B. Categorisation of the Anomalies

Each sensor provides timestamped data. For each one, the actual value is the measured data. A nominal value is calculated from a profile based on the values usually provided by the sensor. The location of the sensors is unknown, as well as any other metadata (type of sensor, acquisition frequency, etc.). As the lack of information on the sensors is a common problem [START_REF] Guerrero | Heterogeneous data source integration for smart grid ecosystems based on metadata mining[END_REF], it is necessary to take them all into account in the method of anomaly detection.

Chandola et al. [START_REF] Chandola | Anomaly Detection: A Survey[END_REF] define an anomaly as an unexpected or undesirable behaviour in a system. The authors define 3 kinds of anomalies:

• point: a measure is outside of an acceptable range for the sensor,

• contextual: a measure is inside an acceptable range for the sensor but anomalous in some context (example: high heating consumption during summer), • collective: a collection of measures is anomalous with respect to the entire data set although the individual measures may not be anomalous in themselves. I illustrates the 3 types of anomalies defined above. They are indicated by the letter A in the last column. Row 1 is considered anomalous by the expert. Sensor 4 has a value of -5, which is far from the nominal value of 11. This is a point anomaly.

Row 8 is an anomaly because of the values of sensors 1 to 4. Although taken individually, none of these values would be considered an anomaly, together they generate a collective anomaly.

The last anomaly in row 13 appears at 13:00 and only then. This is due to the fact that the actual values of sensor 3 are disparate from line 9 to line 13. In real life, an expert may have noticed a problem at 13:00 due to the constant accumulation of small differences in the previous 4 lines. This is called a contextual anomaly because it is necessary to take into account the recent past to detect it.

Row 5 is not an anomaly even if the value of 400 in sensor 5 is very far from the nominal value of 65. It is not considered as an anomaly by the expert, because this sensor is either irrelevant for the detection of anomalies or that the variation is not large enough for an expert to consider it as an anomaly.

II. STATE OF THE ART OF ANOMALY DETECTION IN SMART-BUILDINGS

Given the real challenge of detecting anomalies in buildings, the field of research in this area is rapidly expanding. Multiagent system (MAS) has been introduced in building physics over the last ten years or so, but their applications are almost entirely dedicated to optimising system management (heating, cooling, ventilation) or the entire building [START_REF] Labeodan | On the application of multi-agent systems in buildings for improved building operations, performance and smart grid interaction -A survey[END_REF]. No study, to our knowledge, has been conducted for an application in the detection of anomalies in building data, even though MASs have been used to detect anomalies in other domains [START_REF] Forestiero | Self-organizing anomaly detection in data streams[END_REF] [START_REF] Ng | Multi-agent based collaborative fault detection and identification in chemical processes[END_REF].

Methods used in the literature include physical models, which require modeling the system or building and comparing it with measured data. These methods are quickly limited, given the complexity of the systems to be modeled [START_REF] Turner | Residential HVAC fault detection using a system identification approach[END_REF].

Unsupervised and statistical classification methods have been used to detect anomalies in building data [START_REF] Miller | A review of unsupervised statistical learning and visual analytics techniques applied to performance analysis of non-residential buildings[END_REF] [8] [START_REF] Chou | Real-time detection of anomalous power consumption[END_REF]. However, these studies are only based on a small amount of data and do not show that scaling is possible. They are also based on the assumption that the number of anomalies is small relative to the total number of data.

Neural networks are tools widely used in building studies, but rarely for the detection of anomalies [START_REF] Wu | Hierarchical Temporal Memory method for time-series-based anomaly detection[END_REF] [START_REF] Zhu | Fault diagnosis for sensors in air handling unit based on neural network pre-processed by wavelet and fractal[END_REF]. Although neural networks can detect anomalies in real time, a large set of labeled data is required to train the networks. However, this step seems impossible to perform for the size of systems considered in this study and it does not correspond to what actually happens in building management.

The methods of data mining are used to detect anomalies in building data [START_REF] Capozzoli | Fault detection analysis using data mining techniques for a cluster of smart office buildings[END_REF] [13] [START_REF] Xue | Fault detection and operation optimization in district heating substations based on data mining techniques[END_REF]. However, a data pre-processing step is required, as well as an expert for the selection of relevant data during scaling.

III. SANDMAN ARCHITECTURE

A. Definitions of Situation and Sensor Profile

SANDMAN works with all the available sensor types and thanks to the feedback given by a human expert. For example, in Table I, SANDMAN has access to the values of the 5 sensors. We define the notion of situation as the set of values measured by all these sensors over 24 hours and we call current situation the last situation encountered. The last measurement of all the sensors is called instantaneous values of the sensors: this corresponds to the last row in Table I, containing the instantaneous values of the 5 sensors at 23:00. The situations encountered by SANDMAN are stored in a situation history.

Each sensor is associated with a profile, consisting of 24 measured values, one for each hour of the day. This profile indicates the value expected at any time of the day for the sensor in question. It is updated by SANDMAN throughout the life of the system.

B. Detection of Anomalies by SANDMAN

Since all available sensors are used and the expert may only be interested in certain faults, and since each sensor has its own interval and value range, profile learning alone is not sufficient to detect anomalies. Indeed, some sensors may not be relevant for the detection of anomalies. This is the case, for example, of sensors located in a building that the expert is not interested in, or if the expert is only interested in power consumption then an air quality sensor will not be relevant. Lastly, for sensors that are actually useful, depending on the sensors, some deviations from expected values are greater than others. This is due to the fact that their values are not normalised because SANDMAN uses raw data.

To raise or not to raise an anomaly at each new current situation, SANDMAN uses both the measured values of the sensors over 24 hours and the profiles of the associated sensors. To do this, SANDMAN calculates the disparity of a sensor which is the sum of all the differences between the sensor profile and the sensor measured values over the last 24 hours, i.e. for the current situation (Eq.( 1)).

Disparity t s = t-23 ti=t realV alue ti s -nominalV alue ti s (1) with :
• s: sensor s;

• t: time t of the current situation;

• realV alue ti s : real value of sensor s at time t i ; • nominalV alue ti s : nominal value of sensor s at time t i ; The period under consideration is a sliding 24-hour window. If the time of the current situation is 15:00, all data from 16:00 of the previous day to 15:00 of the current day are considered for the profile. SANDMAN then calculates the Degree of Anomaly (DA) of the current situation which is the weighted sum of the measured values of all sensors (Eq. ( 2)).

DA(Situation

t ) = S s=1 Disparity t s * W eight s (2)
with :

• Situation t : Situation at time t;

• S: number of sensors; • Disparity t s : Disparity of the sensor s at time t;

• W eight s : Weight of the sensor s.

Thus thanks to the degree of anomaly, SANDMAN classifies the current situation, i.e. labels it "normal" or "anomalous" by algorithm 1. The threshold used is a static number Algorithm 1 Algorithm of classification or anomalies detection 1: for each new current situation do 2:

if DA < Threshold then SANDMAN returns : "the situation is anomalous" 6:

end if 7: end for chosen arbitrarily (1000 for SANDMAN). The calculation of each weight adjusts to the threshold whatever its value, and SANDMAN classifies the situation (on the correct side of the threshold) according to the value of this weight. Therefore, the DA also depends of the threshold value.

To correctly detect anomalies, the W eights i associated with the sensors must therefore have the correct values. A sensor has only one associated weight that is used to calculate the degree of anomaly for all situations in the situation history.

C. Learning by SANDMAN

Once SANDMAN has classified a situation, i.e. labeled it as "normal" or "anomalous", and received feedback from the expert, it performs a resolution cycle in which it analyses whether or not it should learn. SANDMAN therefore analyses the following cases:

1) For a situation classified as "anomalous" by SANDMAN and "anomalous" by the expert, a) SANDMAN does nothing 2) For a situation classified as "normal" by SANDMAN and "anomalous" by the expert, a) SANDMAN creates and adds a new "anomalous" situation to the situation history b) SANDMAN self-adapts the weights 3) For a situation classified as "anomalous" by SANDMAN and "normal" by the expert, a) SANDMAN updates the sensor profiles b) SANDMAN creates and adds a new "normal" situation to the situation history c) SANDMAN self-adapts the weights 4) For a situation classified as "normal" by SANDMAN and "normal" by the expert, a) SANDMAN updates the sensor profiles. The main steps for updating profiles, adding a situation in the history and self-adjusting weights are described below.

Profile updates. They are only updated if the situation is labelled "normal" by the expert (steps 3a and 4a), because the sensor values in anomalous situations are not reliable and do not represent the expected values. The update formula applied to each sensor profile is as follows (Eq.( 3)) : Each sensor profile modifies one of its 24 values corresponding to the time t of the current situation using the measured value val t sensor of the current situation.

val t prof ile = (1 -λ) *
Adding a situation to the situations history. This operation, performed in steps 2a and 3b, initiates SANDMAN resolution cycles that occur until all situations in the history are correctly and evenly classified. A situation is classified correctly if its degree of anomaly (DA) enables to deduce the correct classification (normal or anomalous) of the situation. Situations have to be classified in a balanced manner: thus the degree of anomaly of the "normal" situation closest to the threshold and the degree of anomaly of the "abnormal" situation also closest to the threshold are, in absolute value, at the same distance from the given threshold. This is illustrated in Figure 2, where the points respectively above and below the threshold T hreshold and closest to the threshold are at equal distance from the threshold. Each situation is assigned a criticality value according to its degree of anomaly and its classification given by the expert. Criticality represents the degree of satisfaction of the situation with respect to its objective. A situation labeled "normal" aims at a degree of anomaly as low as possible with respect to the T hreshold threshold, while a situation labeled "anomalous" aims at a high degree of anomaly with respect to the T hreshold threshold. The criticality Crit is calculated as for a "normal" situation (Eq.( 4)) and for an "anomalous" situation (Eq.( 5)).

Crit Situation t = DA(Situation t ) -T hreshold (4)

Crit Situation t = T hreshold -DA(Situation t ) (5)
Every situation wants to minimise its criticality, a negative criticality meaning that the situation is correctly classified. Since disparities and weights are always positive, situations labeled "normal" and situations labeled "anomalous" have antagonistic objectives, wanting respectively higher or lower weights. When a resolution cycle is performed, each situation sends a message to each weight agent with:

• the desired direction of weight change (lower for "normal" situations, higher for "anomalous" situations); • the criticality of the situation;

• the influence of the weight in the degree of anomaly. The influence Inf lu of a sensor s in a situation t corresponds to the importance of the sensor weight on the sum of the weight of all the sensors for this situation, as shown in eq.( 6):

Inf lu t s = W eight s * Disparity t s DA(Situation t ) (6)
A resolution cycle starts when situations require the weights to adjust. The cycle ends when the weights have self-adjusted.

Self-adapting weights. Each weight is represented by a weight agent. The goal of the multi-agent weight system is that each weight agent finds its value by cooperating with the others. During a resolution cycle, all situations require weight adjustments (steps 2b and 3c) and the weight agents self-adapt as described in the algorithm 2, then each situation calculates its new degree of anomaly using the updated weights.

On receiving the parameters sent by the situations (direction, criticality and influence), the weight agents decide, or not, to update independently and simultaneously their weights. To do so, they follow the algorithm 2. The Weight agent raises its weight 9:

end if 10: end if

The cooperation between the weight agents guarantees that at least one weight agent updates its weight and that each update cycle of the weight agents results in a better general state where the most critical weight agent decreases its criticality [START_REF] Georgé | Cooperation[END_REF]. After each updating cycle of the weight agents, the situations calculate their new criticality by taking into account the new values of the weight agents. This process is repeated until the most critical situations labeled "normal" and "anomalous" have their respective criticalities both equal (within a small margin of error) and negative. This ensures that all situations are correctly classified and that the following situations have a better chance of being classified accordingly by SANDMAN. Indeed, the best way to reduce the rate of false negatives and false positives [START_REF] Chandola | Anomaly Detection: A Survey[END_REF], and thus to have a robust system, is to keep the most critical known situations (saved in the history) as far from the threshold as possible.

IV. EXPERIMENTATIONS

A. Description of the Experimental Framework

The data used to conduct the experiments were generated using the TSimulus time series generator [START_REF]Tsimulus[END_REF]. These time series include a signed value for each sensor for each hour. Each sensor has its own range of values. The values for each sensor are cyclical, with or without noise over a 24-hour period, i.e. the value of a sensor (without noise) at 3 p.m. is the same every day. These simulated data are modified by a human expert to introduce anomalies of different types. SANDMAN processes the data in constrained time and the expert gives his opinion a posteriori in an asynchronous way.

All the experiments were performed on a 4-core processor with a 2.6 GHz frequency. The data sets used in the following experiments feature one value for each sensor every hour over one month, for a total of 744 hours. The number of sensors depends on the experiment. Four separate data files are used to show the ability of SANDMAN to 1) detect point anomalies and mitigate noise, 2) detect collective and contextual anomalies, 3) adapt the profiles of sensors, 4) scale.

In the experiments, we chose to duplicate the one-month data set in order to show the improvement granted by the learning of SANDMAN. Thus, the results are presented after one month and two months, as the number of:

• TP: True Positive, TN: True Negative, • FN: False Negative, FP: False Positive, • t/sit : calculation time for processed situation in ms.

B. Point Anomalies and Noise Mitigation

Point anomalies are the type of anomaly most often detected in the literature, as both univariate and multivariate approaches can detect them.

In this experiment, we study the efficiency of SANDMAN in detecting point anomalies in noisy data with different levels of noise. The data file contains 20 sensors, i.e. 20 x 744 = 14880 data. 58 point anomalies were added manually by an expert. Three experiments have been conducted, with the noiseless data, 1%-level data and 5%-level data.

The noise is simulated by applying a uniform distribution to each sensor with the following formula (Eq [START_REF] Miller | A review of unsupervised statistical learning and visual analytics techniques applied to performance analysis of non-residential buildings[END_REF]): Table II shows the results of the 3 experiments after month 1 and month 2. We can see that SANDMAN improves its detection with time. After month 2, it detects the right anomalies (58 TP) when the data have a low noise level. However, a higher noise level leads to a higher false positive rate (7 FP after month 2 in the case of 5% noise). SANDMAN reduced this false positive rate along time by creating normal situations, and thus by learning lower weights, as shown in Algorithm 2.

V alue noise = V alue noN oise - |M IN (s) -M AX(s)| 2 +rand [0,1] * interval * noise (7)

C. Collective and Contextual Anomalies

In this experiment, we want to study the ability of SAND-MAN in detecting different types of anomalies, such as point, contextual and collective anomalies. Contextual anomalies occur when a sensor value is unusual several times in a row. Collective anomalies occur when several sensor values are unusual at the same time, but not so unusual that only one of them causes a single point anomaly. In this experiment, the data file contains the values of 20 sensors over one month, and the results are presented after months 1 and 2. The 3 types of anomalies (point, collective and contextual) have been manually added in the data file, for a total of 14 anomalies. The data have 1% of noise.

Table III shows the results of this experiment. The numbers of each type of anomaly, as well as the anomalies detected after the month 1 and then after the month 2. After month 2, SANDMAN detected the right anomalies. There are no false positives in the results of this experiment. SANDMAN created and added to the history the situations corresponding to the misclassified situations during the month 1 to classify them without error during the month 2, regardless of the type of anomaly.

D. Updating Profile on Sliding Data

In previous experiments, the profile of each sensor was stored and contributed to the detection of anomalies, but the data was cyclical on a daily basis. In this experiment, the measured values of the sensors are modified so that each measured value is 1% higher each day, so that over a 31-day month, each sensor has measured values that are 31% higher than the control data. The test file contains 58 point anomalies and each of the 20 sensors causes an anomaly at least once. Table IV presents the results of the experiment. The update of the sensor profile is sufficiently reactive to allow SANDMAN to classify anomalies when the measured values of the sensors change over time.

E. Scaling

In this experiment, an increasing number of sensors is used to measure this effect on the computing time. To do this, a data set of 20 separate sensors is duplicated to obtain up to 800 sensors. The data used is the same as for the previous experiment on rolling data without the change in values of 1% per day. Table V shows the calculation time per situation as a function of the number of sensors. The calculation time of month 2 is always shorter than the calculation time of month 1 because the learning of the weights is only carried out in month 1. We also note that the resolution time is proportional to the number of sensors used and that the time difference between the two months is constant. This is due to the fact that most of the execution time comes from reading the raw data from a database, which has a fixed cost per sensor.

V. CONCLUSION AND PERSPECTIVES

For the management of energy-related data in an intelligent building, several types of anomalies must be detected. They can be point, collective, or contextual anomalies. The detection of these different types of anomalies is absent from the state of the art methods for which data pre-processing is mandatory.

We, therefore, proposed SANDMAN, a semi-supervised real-time anomaly detection system that uses raw data as input and classifies anomalies by learning from the expert's feedback. SANDMAN is able to detect three kinds of anomalies in a generic way and is well suited to a growing number of sensors. The feedback from the expert is reduced to a minimum because it is optional and SANDMAN can be operated without prior labelling.

In future works, after an anomaly occurs the system must be able to observe the return to nominal values of the sensors to inhibit instantaneously the anomalies. This is the improvement we are currently developing. SANDMAN will also have to be able to learn several sets of profiles for the sensors, in order to take into account the different behaviours of building users (for example, weekdays and weekend periods).

Fig. 1 .Figure 1

 11 Fig. 1. 2 sets of measured values from the same sensors over 24h

Fig. 2 .

 2 Fig. 2. Example of 8 classified situations

Algorithm 2

 2 Self-adaptation of the Weight agent of the sensor s Require: Direction, Crit, Inf lu of all situations 1: Select the most critical "normal" situation N S 2: Select the most critical "anomalous" situation AS 3: if Crit(N S) > Crit(AS) then

TABLE I EXAMPLE

 I

OF DATA AND EXPERT'S FEEDBACK

Table