Amy Lynd

David Weetman

Susana Barbosa

Alexander Egyir Yawson

Sara Mitchell

Joao Pinto

Ian Hastings

Martin James Donnelly
email: m.j.donnelly@liv.ac.uk

Alexander Egyir-Yawson

Field, genetic and modelling approaches show strong positive selection acting upon an

Keywords: Anopheles gambiae, malaria, insecticide resistance, selection

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

Introduction

The sequencing of the Anopheles gambiae genome has opened up the possibility for genomewide SNP-based association mapping studies that have been successful in identifying positively selected loci in the human genome [START_REF] Sabeti | Detecting recent positive selection in the human genome from haplotype structure[END_REF][START_REF] Bersaglieri | Genetic signatures of strong recent positive selection at the lactase gene[END_REF][START_REF] Sabeti | Genome-wide detection and characterization of positive selection in human populations[END_REF]). The resolution of the association mapping approach is defined by the probability that recombination will have broken down the association between markers and a trait-associated functional polymorphism. Data from extensive re-sequencing of (primarily) detoxification genes in samples from wild populations of An. gambiae revealed a very high frequency of segregating sites [START_REF] Wilding | High, clustered, nucleotide diversity in the genome of Anopheles gambiae revealed by SNP discovery through pooledtemplate sequencing: implications for high-throughput genotyping protocols[END_REF], consistent with high rates of recombination [START_REF] Begun | Levels of naturally occurring DNA polymorphism correlate with recombination rates in D. melanogaster[END_REF][START_REF] Begun | Population genomics: Whole-genome analysis of polymorphism and divergence in Drosophila simulans[END_REF]) and/or a long history of outbreeding. In isofemale lab strains of Drosophila spp. it has been possible to observe selective sweeps around insecticide-resistance associated loci [START_REF] Schlenke | Strong selective sweep associated with a transposon insertion in Drosophila simulans[END_REF][START_REF] Aminetzach | Pesticide resistance via transposiitonmediated adaptive gene truncation in Drosophila[END_REF] but how long these signatures persist in wild populations is unknown. In this paper we use linkage disequilibrium-based haplotype diversity analysis [START_REF] Sabeti | Positive natural selection in the human lineage[END_REF] to investigate the pattern of molecular genetic variation associated with insecticide resistance mutations at the pyrethroid and DDT knock down resistance locus, kdr, in the African malaria mosquito An. gambiae s.s. Furthermore, as a corollary of this indirect genetic approach we demonstrate, using a series of temporal collections, a dramatic increase in kdr frequency in a population of An. gambiae s.s. over a period of approximately 72 generations. Data from these temporal collections are used to estimate the selection and dominance coefficients operating on kdr in the field to illustrate the potential levels of selection necessary to produce the patterns of linkage disequilibrium (LD) we observe.

Insecticide-treated bed nets are the principal method for preventing malaria in sub-Saharan Africa. Currently, pyrethroids are the only class of insecticides licensed for use on nets and there is concern that resistance will compromise control programmes. To date the most commonly recorded resistance mechanism is termed knock down resistance, and results from single base pair mutations in the voltage-gated sodium channel. The sodium channel gene, located within division 20C near the centromere of chromosome 2L, codes for a protein which is the target-site of pyrethroid insecticides. Two alternative single base pair mutations have been found in An. gambiae and these kdr mutations can cause target site insensitivity to pyrethroids as well as cross-resistance to DDT. The substitutions cause amino acid changes at codon 1014 within the transmembrane structure of segment 6 in domain II of the voltagegated sodium channel (numbering according to the housefly para sequence, Genbank X96668). The L1014F mutation, a leucine to phenylalanine change, was first observed in West Africa [START_REF] Martinez-Torres | Molecular characterization of pyrethroid knockdown resistance (kdr) in the major malaria vector Anopheles gambiae s.s[END_REF]) and the same substitution has been observed in a diverse array of insects (Davies et al. 2007a). A second substitution, L1014S, was observed more recently in East African An. gambiae [START_REF] Ranson | Identification of a point mutation in the voltage-gated sodium channel gene of Kenyan Anopheles gambiae associated with resistance to DDT and pyrethroids[END_REF] and involves the adjacent base of the same codon, resulting in a leucine to serine change.

There are two incipient species within the nominal taxon Anopheles gambiae s.s. which are characterized by mutations on the X-chromosome and are termed M and S form. The distribution of the kdr mutation is not uniform either within or between forms, although in general kdr alleles have been found at much higher frequencies in Anopheles gambiae s.s. Sform samples compared to M-form samples [reviewed in [START_REF] Santolamazza | Distribution of knock-down resistance mutations in Anopheles gambiae molecular forms in west and west-central Africa[END_REF]]. The reasons for the differences in distribution remain unclear since little is known about the origins of the kdr mutations and the selection pressures acting upon them in wild populations.

In a sample from Benin, the L1014F was found in tight linkage disequilibrium with two upstream intronic polymorphisms in both M and S-form individuals. The two upstream polymorphisms associated with the L1014F variant were not found in wildtype M-form individuals but were common in wildtype S form individuals, suggestive of an introgression event from S-form to M-form populations [START_REF] Weill | The kdr mutation occurs in the Mopti form of Anopheles gambiae s.s. through introgression[END_REF]. This linkage between kdr and the intronic polymorphisms was not seen in M form individuals from Bioko Island and was thought to indicate de-novo mutation [START_REF] Reimer | Relationship between kdr mutation and resistance to pyrethroid and DDT insecticides in natural populations of Anopheles gambiae[END_REF]. More recently a study of S-form specimens from 15 countries suggested that the L1014F and L1014S mutations have both arisen independently on at least two separate occasions [START_REF] Pinto | Origins and distribution of knockdown resistance mutations in the afrotropical mosquito vector Anopheles gambiae[END_REF].

Samples were obtained from three regions in sub-Saharan Africa; Kenya (East Africa) Anopheles gambiae, S molecular form, kdr L1014S allele present; Ghana (West Africa) both M and S molecular form, kdr L1014F allele present; Gabon (Central Africa) S molecular form, both the L1014S and L1014F kdr alleles present.

These population samples allow us to address a number of questions.

1. Available evidence suggests that the L1014S mutation has high penetrance for a DDTresistant phenotype but lower penetrance for a pyrethroid-resistant phenotype than the L1014F mutation [START_REF] Ranson | Identification of a point mutation in the voltage-gated sodium channel gene of Kenyan Anopheles gambiae associated with resistance to DDT and pyrethroids[END_REF]. DDT was banned in Kenya in 1990 and we can investigate the signature of positive selection associated with weaker selection or recombination and relaxed selection.

2. The populations from central Africa are some of the few locations where both the L1014F and L1014S alleles are observed sympatrically [START_REF] Santolamazza | Distribution of knock-down resistance mutations in Anopheles gambiae molecular forms in west and west-central Africa[END_REF].

Indeed in an earlier study a significant, albeit marginal, L1014F/L1014S heterozygote excess was observed in samples from Libreville, Gabon [START_REF] Pinto | Co-occurrence of East and West African kdr mutations suggests high levels of resistance to pyrethroid insecticides in Anopheles gambiae from Libreville, Gabon[END_REF]. By comparing patterns or linkage disequilibrium around the three alleles we investigate whether the unusually high frequency of the L1014S allele in these populations (63% [START_REF] Pinto | Co-occurrence of East and West African kdr mutations suggests high levels of resistance to pyrethroid insecticides in Anopheles gambiae from Libreville, Gabon[END_REF])) is a result of a recent selective sweep.

3. In many S-form populations in West Africa, including our collections from Ghana, the L1014F allele is close to fixation. In the absence of wild-type alleles we are unable to control for local variation in recombination rates [START_REF] Sabeti | Genome-wide detection and characterization of positive selection in human populations[END_REF]) and it is therefore impossible to ascribe patterns of LD to a positive selection event.

Recently developed approaches such as Cross Population Extended Haplotype Homozygosity have been developed to allow interpopulation comparisons in instances where alleles have proceed to near fixation in some populations [START_REF] Sabeti | Genome-wide detection and characterization of positive selection in human populations[END_REF]) but in our system resistance alleles may have multiple origins, presenting a confounding variable [START_REF] Pinto | Origins and distribution of knockdown resistance mutations in the afrotropical mosquito vector Anopheles gambiae[END_REF]). However, the presence of sympatric M-form populations in southern Ghana [START_REF] Yawson | Species abundance and insecticide resistance of Anopheles gambiae in selected areas of Ghana and Burkina Faso[END_REF][START_REF] Yawson | Ecological zones rather than molecular forms predict genetic differentiation in the malaria vector Anopheles gambiae s.s. in Ghana[END_REF]) allow us to both document the increase in frequency of the same L1014F haplotype, following an introgression event, over a period of 5 years; and to estimate the selection and dominance coefficients associated with the signatures of positive selection.

Materials and methods

Sample sites, DNA extraction and species identification.

Adult female Anopheles gambiae s.s. mosquitoes used in this study were obtained from aspirator and pyrethroid knock-down collections from the field in various geographic locations (table 1). DNA was extracted from single female Anopheles gambiae using either a modified Livak method or a phenol-chloroform method [START_REF] Livak | Organization and mapping of a sequence on the Drosophila melanogaster X and Y chromosomes that is transcribed during spermatogenesis[END_REF][START_REF] Ballinger-Crabtree | Use of genetic polymorphisms detected by the random-amplified polymorphic DNA polymerase chain reaction (RAPD-PCR) for differentiation and identification of Aedes aegypti subspecies and populations[END_REF]. Species identification PCR, was carried out on Anopheles gambiae s.l. according to the protocol [START_REF] Scott | Identification of single specimens of the Anopheles gambiae complex by the Polymerase Chain Reaction[END_REF]. Reactions were then digested with CFO I restriction enzyme for 24 hours at 37°C in order to type Anopheles gambiae s.s mosquitoes to M and S-form [START_REF] Fanello | Simultaneous identification of species and molecular forms of the Anopheles gambiae complex by PCR-RFLP[END_REF], and products visualised under UV light after electrophoresis on a 2% agarose TBE gel with Ethidium Bromide. Kdr genotypes were determined by allele specific PCR, Heated Oligonucleotide Ligation Assay (HOLA) [START_REF] Lynd | A simplified high-throughput method for pyrethroid knock-down resistance (kdr) detection in Anopheles gambiae[END_REF] or Taqman assay [START_REF] Bass | Detection of knockdown resistance (kdr) mutations in Anopheles gambiae: a comparison of two new high-throughput assays with existing methods[END_REF]) depending upon year of collection.

Sodium channel SNP identification

The voltage-gated sodium channel gene is nearly 74Kbp in length and is composed of 35 exons including two duplicate exons (Davies et al. 2007a). Ten regions of the sodium channel were amplified by PCR for direct sequencing. Where possible, primers were designed to bind within exons to produce amplicons that spanned an intron with a maximum size of 1.5Kbp.

Exons (numbering as Davies et al. 2007a) 1-2, 3, 4, 7-9, 13-14, 15-17, 20c, 23-24, 28-30 and 32-33 were selected as targets for sequencing. Primer and amplification details are provided (supplementary table 1, Supplementary Material online). Sequencing for SNP detection was carried out on up to 12 individuals of known kdr genotype from Ghana, São Tomé, Gabon, Angola, Mozambique, Malawi, Kenya, from a susceptible laboratory strain (KISUMU) and from a permethrin tolerant resistant laboratory strain (RSP), both originating from Kenya. PCR products were cleaned using a Mini Elute PCR Purification kit (Qiagen) and then sequenced in both directions. Sequences were aligned using Bioedit software version 7.0.5.2 [START_REF] Hall | BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF]) then manually annotated for polymorphisms and ambiguities.

In addition, seven M-form individuals from Accra, Ghana, homozygous for the L1014F allele were bidirectionally-sequenced across PCR amplicons 13-14, 15-17 and 21 to determine the associated haplotype of the kdr allele in this population.

SNP screening

SNPs discovered through resequencing were screened in the large scale SNP detection study using the SNPStart Primer Extension Kit on the Beckman CEQ™ 8000 Genetic Analysis System. Details of SNPs both included and excluded from the SNP screening are given in supplementary table 2 (Supplementary Material online). Multiplex PCR was carried out to amplify the regions of DNA containing SNPs of interest, including a region of exon 20 and the preceding intron to allow high throughput detection of the kdr mutation and three other well characterized SNPs [START_REF] Weill | The kdr mutation occurs in the Mopti form of Anopheles gambiae s.s. through introgression[END_REF][START_REF] Diabate | The spread of the Leu-Phe kdr mutation through Anopheles gambiae complex in Burkina Faso: genetic introgression and de novo phenomena[END_REF][START_REF] Pinto | Co-occurrence of East and West African kdr mutations suggests high levels of resistance to pyrethroid insecticides in Anopheles gambiae from Libreville, Gabon[END_REF]) (primers and reaction conditions detailed in supplementary table 3, Supplementary Material online).

Products were visualized on a 2% TBE agarose gel. Successfully multiplexed samples were prepared for subsequent SNP extension by ExoI/SAP enzymatic digestion. Interrogation primers were then designed for each individual SNP chosen for investigation according to the manufacturers recommendations (supplementary table 4, Supplementary Material online).

Single base extension to the 3' end of the interrogation primer by a dye terminator molecule, corresponding to the nucleotide found at the SNP location, was carried out using a GenomeLab SNPStart Primer Extension Kit (Beckman Coulter, Amersham, UK). The SAPdigested product was then scored on the Beckman CEQ 8000 Genetic Analysis System.

Data analysis

As reviewed exhaustively by [START_REF] Sabeti | Positive natural selection in the human lineage[END_REF] there are numerous statistical tests of positive selection which differ in their ability to detect selection events on different time scales. For the present SNP dataset it is not possible to use the suite of sequence-based tests that compare synonymous/ non-synonymous differences or detect an excess of rare alleles.

We are therefore fortunate that on the timescales in which the emergence, and selection, of insecticide resistance is likely to occur, estimates of interpopulation divergence (e.g. based on F-statistics) and screens of LD around selected vs wildtype alleles are likely to be the two most powerful analytical approaches. With the sample sizes available in our study, single marker analyses based on F-statistic estimates would perform better as indicators of selection when markers can be typed at a more coarse scale, with consequently enhanced signal: noise ratio. However, with sample size constraints the signal would be difficult to localise. By contrast long range haplotype analyses such as extended haplotype homozygosity (EHH) [START_REF] Sabeti | Detecting recent positive selection in the human genome from haplotype structure[END_REF] analysis performs very well at a fine physical scale in identifying narrow candidate regions [START_REF] Sabeti | Positive natural selection in the human lineage[END_REF].

EHH analysis was carried out to assess the patterns of LD associated with wild type and the two kdr alleles. EHH can be defined as the probability that two random chosen chromosomes carrying the core (eg. the wildtype or kdr allele) haplotype of interest are identical by descent. This approach first identifies core haplotypes surrounding the locus of interest and then examines the decay in LD from these core haplotypes to the surrounding loci. The resulting EHH can be used as evidence of recent positive selection at a locus in haplotypes that have high frequency and high EHH [START_REF] Sabeti | Detecting recent positive selection in the human genome from haplotype structure[END_REF]. EHH analysis requires haplotype information that cannot be empirically determined from the genotype data gathered by the methods used in this study. Therefore haplotypes were inferred using PHASE software version 2.1.1 using default parameters [START_REF] Stephens | A new statistical method for haplotype reconstruction from population data[END_REF][START_REF] Stephens | Accounting for decay of linkage disequilibrium in haplotype inference and imputation[END_REF]. PHASE utilizes a Bayesian coalescent-based approach to determine phase and allows for varying rates of recombination at each SNP interval. The method is based on the idea that an unresolved haplotype is more likely to be the same or be similar to a previous haplotype. This approach was found to outperform other methods available for autosomal human data sets [START_REF] Stephens | A new statistical method for haplotype reconstruction from population data[END_REF][START_REF] Stephens | Accounting for decay of linkage disequilibrium in haplotype inference and imputation[END_REF]. Data were analysed together rather than as separate subpopulations since (a) previous studies found this to be more accurate, and (b) haplotype determination methods of this nature are relatively insensitive to departures from Hardy-Weinberg equilibrium so are fairly robust to population substructuring. This approach is also more conservative than determining haplotypes for individual populations since the latter is liable to lead to an underestimation in differences in haplotype frequencies [START_REF] Stephens | Accounting for decay of linkage disequilibrium in haplotype inference and imputation[END_REF]. Phase reconstruction was executed ten times upon the total data set and differences in counts of best haplotypes were noted.

The estimated haplotypes obtained from PHASE were used as input for EHH analysis implemented by SWEEP version 2.1.1 [START_REF] Sabeti | Detecting recent positive selection in the human genome from haplotype structure[END_REF]. Core haplotypes were selected manually to include only the two adjacent kdr-causing loci. Significance of EHH values is usually assigned through comparison to an empirically generated null distribution from other regions of the genome. However given that we had already identified the causal mutations of interest we were able to make a comparison of patterns of LD around wild-type and resistant cores. The primary advantage of this approach is that it is not subject to the genome-wide variations in recombination rate which can affect the null distribution approach in species lacking detailed recombination maps. Significant differences in EHH values were determined in two ways. (1) Within country samples, at individual SNP positions with non-overlapping 95% confidence intervals. These confidence intervals were calculated at each SNP position using a bootstrapping procedure, carried out in SAS Version 9 software. Re-sampling was carried out 1000 times. (2) Across all SNPs within and among country samples the diversity of the different kdr allele-bearing haplotypes were compared using sign tests, implemented by SPSS 14. Where exact sign test probabilities could not be calculated, a Monte Carlo procedure with 10000 permutations was performed. The sequential Bonferroni procedure was applied to determine statistical significance following correction for multiple testing [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]. Although our data -EHH values at each SNP position -are not independent it is this non-independence caused by LD that will cause departure from the null hypothesis of equality of median EHH values. Therefore the null hypothesis remains that there is no difference in median EHH between kdr and wild type alleles. Bifurcation plots were also created using the SWEEP software. In a bifurcation plot the core haplotype is represented as a black circle. Each SNP, moving out from the core both upstream and downstream of the kdr locus, is a potential site for a bifurcation which would result from the presence of two segregating alleles. Therefore the diagram provides a means of displaying the breakdown in LD at increasing distance from the core haplotypes. The radius of the circle at each node is proportional to the number of individuals with that haplotype.

Calculation of selection and dominance coefficients

The spread of the L1014F allele was modelled using the standard recursive population genetic formula:

(1)

Where p is the frequency of the L1014F allele, p' is the frequency in the next generation, s is the selective coefficient of the resistance mutation, h is the dominance coefficient (1=complete dominance, 0=complete recessivity) , and W is the normalising factor (Maynard-Smith 1998) .

Tracking allele frequencies over time requires three input parameters: initial allele frequency at time zero; s; and h. Estimates of all three unknown parameters were obtained by maximum likelihood assuming a binomial distribution of observed allele frequencies around the predicted frequency. The analysis was performed in R (http://www.r-project.org) using maximum likelihood functions and optimizing routines. The generation time was set at the standard of one generation per calendar month [START_REF] Lehmann | The effective population size of Anopheles gambiae in Kenya: implications for population structure[END_REF].

Results

SNP discovery and screening.

Ten genomic regions of a combined length of ≈6.5Kb of DNA, spanning a region of ≈73Kb of the voltage-gated sodium channel were amplified and sequenced in Anopheles gambiae s.s.

individuals from seven countries across sub-Saharan Africa. A total of 62 potential SNPs were found, of which 14 were exonic (supplementary table 1, Supplementary Material online). Six intronic indels were observed, usually in poly-A or tandem AT repeats (table 1, Supplementary Material online). On average there was one SNP every 106 bp, which represents a low SNP frequency for An. gambiae, but similar to other genes in the same genomic locality (chromosome 2L division 20; [START_REF] Wilding | High, clustered, nucleotide diversity in the genome of Anopheles gambiae revealed by SNP discovery through pooledtemplate sequencing: implications for high-throughput genotyping protocols[END_REF]). Thirty-two SNPs, including the two kdr mutations, were selected for screening in 258 individuals. In S-form individuals the SNP adjacent to the core in the upstream (centromeric) direction was excluded from further analysis as it was found to be monomorphic. Details of the populations and associated kdr genotypes are given in table 1. The genotypic data were resolved into haplotypes with ten runs of the analysis. In only one instance did the replicate runs resolve a novel estimated haplotype, which in a subsequent comparative analysis was found to exert no qualitative effect on the results. Therefore all analyses reported here are based upon the haplotypes resolved in the vast majority of the phasing runs.

Extended haplotype homozygosity (EHH) analysis was carried out to assess the patterns of LD associated with the wild type and the two kdr alleles. The intronic SNPs which have been used to identify the origin of the kdr mutations were the proximate SNPs in the centromeric direction [START_REF] Weill | The kdr mutation occurs in the Mopti form of Anopheles gambiae s.s. through introgression[END_REF][START_REF] Pinto | Origins and distribution of knockdown resistance mutations in the afrotropical mosquito vector Anopheles gambiae[END_REF]). Linkage disequilibrium decay was examined between these core haplotypes and the remaining 29 or 30 SNP loci (for S-or M-forms, respectively).

Only two core-alleles were present in the western Kenyan sample: wild-type and L1014S. In the downstream telomeric direction, EHH decays at a similar rate for both wildtype and L1014S, but there was a marked contrast between alleles in the centromeric direction, with entirely non-overlapping confidence limits from just a few kb away from the core (fig. 1 and2a). In the Gabonese collection the difference between resistance associated alleles and wild type was even more marked with significantly lower EHH in the wildtype in both centromeric and telomeric directions less than 5kb from the core (fig. 2b). Indeed both of the L1014F and L1014S resistance mutations showed little haplotype bifurcation in the Gabon samples over the length of the sodium channel (fig. 1) suggesting a relatively recent origin for both these mutations accompanied by a strong selective sweep. The patterns of LD are most marked around the resistant L1014F haplotype in Ghanaian S-form samples in which the L1014F kdr allele was at very high frequency (fig. 1 and2c), as would be expected given the near fixation of this allele in southern Ghana in the S molecular form (mean frequency = 0.96; 95% CI 0.95-0.97) [START_REF] Yawson | Species abundance and insecticide resistance of Anopheles gambiae in selected areas of Ghana and Burkina Faso[END_REF]). The presence of only two wildtype haplotypes in the sample prevent any meaningful comparison of LD decay but it should be noted that there was complete LD over the entire 64kb length of the sodium channel in the centromeric direction. The wildtype allele, observed in the Ghanaian M-form populations (fig. 1 and2c) showed marked LD, only in the telomeric direction, between exons 20 and 32, the opposite directional asymmetry to the L1014F mutation in Ghana S-form populations. Whilst simulation studies have shown that LD decay may be asymmetric even when rates of mutation and recombination are constant [START_REF] Kim | Detecting a local signature of genetic hitchhiking along a recombining chromosome[END_REF] it is possible that the LD observed in these samples may reflect the presence of one or more hitherto overlooked selectively advantageous mutants, although we cannot rule out recombination with unsampled haplotypes (supplementary table 5, Supplementary Material online). (Davies et al. 2007b) have summarised that there are a number of additional non-synonymous changes observed in a variety of taxa and detailed association mapping studies are presently underway to investigate this phenomenon. Comparing overall levels of EHH for the whole 72.6 kb regions typed, it is interesting to note that median EHH values are statistically indistinguishable for the same allele typed in different populations (table 2) and that a clear hierarchy of evidence for selective sweeps emerged. Median EHH levels were highest for the L1014F resistance mutation, followed by those for the L1014S mutations, with the lowest for the wildtype allele (table 2). The only exception to this pattern was within the Gabonese sample, the only one in which both resistance alleles were present, where median EHH was equal for the two resistance alleles. Nevertheless, despite the possibilities of different origins of the same allele, and local variation in recombination rates, EHH levels across the genomic region investigated suggest some degree of commonality in selection across populations for each allele, although the actual rate of change in LD with distance can be quite complex and dependent on direction from the core (fig. 1 and2).

We investigated temporal change in the frequency of the L1014F allele and associated haplotype in sympatric populations of M-form individuals in a subset of the populations previously described by [START_REF] Yawson | Species abundance and insecticide resistance of Anopheles gambiae in selected areas of Ghana and Burkina Faso[END_REF]. Using the data reported in [START_REF] Yawson | Species abundance and insecticide resistance of Anopheles gambiae in selected areas of Ghana and Burkina Faso[END_REF] we estimated the L1014F allele frequency in M form populations from around Accra, southern Ghana (≈30km diameter collection area) during 2002 (freq L1014F =0.03; 95% CI 0.01-0.05). Additional screening in 2007 and 2008 from the same greater Accra regions revealed that within five years this frequency had reached freq L1014F =0.54 (95% CI 0.49-0.60; (fig. 3).

The data from years 2007 and 2008 are reported here for the first time. Phasing of the SNP genotypes of two M-form individuals with a wild type/L1014F genotype showed that the L1014F-associated haplotype was identical to that found in the S-form. This was confirmed by sequences obtained from seven M-form individuals collected from Accra, Ghana in 2008, which were homozygous for the L1014F allele (supplementary table 6, Supplementary Material online). Therefore the L1014F allele, which has increased in frequency in M-form populations is the same that has been putatively swept towards fixation in sympatric S-form populations. Introgression of kdr alleles between forms has been documented previously [START_REF] Weill | The kdr mutation occurs in the Mopti form of Anopheles gambiae s.s. through introgression[END_REF] and is unsurprising given that in southern Ghana there is a low but temporally stable level of interform matings [START_REF] Yawson | Species abundance and insecticide resistance of Anopheles gambiae in selected areas of Ghana and Burkina Faso[END_REF][START_REF] Yawson | Ecological zones rather than molecular forms predict genetic differentiation in the malaria vector Anopheles gambiae s.s. in Ghana[END_REF]).

Using a maximum likelihood estimation procedure with random starting values for selection coefficient (s), dominance (h) and initial allele frequency (p 0 ) the parameter estimates converged to s=0.163 (sd= 0.052); h= 0.249 (sd= 0.142) and initial frequency p 0 = 0.025 (sd= 0.008) (fig. 3).

Discussion

These data show that there is marked LD around kdr mutations, loci exhibiting high penetrance, and, for L1014F at least, subject to strong recent positive selection. Despite similar median EHH levels, there were differences in the patterns of LD associated with the L1014S mutation in Kenya and Gabon. In Kenyan samples the rate of dissipation of LD around the L1014S core was quite rapid suggesting that the mutation has not been subject to as recent or as strong a selective sweep as the same mutation in Gabon (or indeed as the L1014F mutation in Ghana). This is as predicted if the serine resistance allele was primarily selected by the use of DDT in the latter part of the 20 th century rather than by the more recent use of pyrethroids in agriculture and insecticide control programmes. In Culex mosquitoes the equivalent L1014S mutation gives low levels of knockdown resistance to pyrethroids compared to the L1014F mutation but confers high levels of DDT resistance [START_REF] Martinez-Torres | Molecular characterization of pyrethroid knockdown resistance (kdr) in the major malaria vector Anopheles gambiae s.s[END_REF][START_REF] Ranson | Identification of a point mutation in the voltage-gated sodium channel gene of Kenyan Anopheles gambiae associated with resistance to DDT and pyrethroids[END_REF]. [START_REF] Stump | Dynamics of the pyrethroid knockdown resistance allele in western Kenyan populations of Anopheles gambiae in response to insecticide-treated bed net trials[END_REF] investigated the change in allele frequency of the L1014S allele before and after the commencement of a large-scale ITN project in Asembo Bay Western Kenya, the site of our collections [START_REF] Stump | Dynamics of the pyrethroid knockdown resistance allele in western Kenyan populations of Anopheles gambiae in response to insecticide-treated bed net trials[END_REF]). The frequency of the L1014S allele in the region approximately ten years before bednet introduction was approximately 0.04 (95% CI 0.02-0.08). In 2002, fifteen years after this initial survey and five years after the introduction of nets the frequency of the L1014S allele had increased, non-significantly to only 0.075 (95% CI 0.05-0.12). This suggests that there is little selective advantage for this mutation in the present environment, although it should be noted that in a neighbouring district in Uganda a recent study reported that the L1014S mutation was at a frequency of 0.85 (95% CI 0.83-0.87) [START_REF] Ramphul | Insecticide resistance and its association with target-site mutations in natural populations of Anopheles gambiae from eastern Uganda[END_REF]). An alternative explanation would be that in Uganda there is an epistatic interaction between L1014S and some, as yet unidentified locus, which may affect the selection, and indeed dominance coefficients, and thereby result in a higher L1014S frequency.

The high frequency and marked LD associated with L1014S in Gabon may be a result of the co-occurrence in genotypes, though not haplotypes (supplementary table 5, Supplementary Material online), with L1014F. A recent study from Cameroon showed that whilst L1014F/L1014S heterozygotes were significantly less resistant to permethrin than L1014F homozygotes, L1014F/L1014S heterozygotes were significantly more resistant to all insecticides tested than L1014F/L1014-wild type heterozygotes [START_REF] Reimer | Relationship between kdr mutation and resistance to pyrethroid and DDT insecticides in natural populations of Anopheles gambiae[END_REF].

Repetitive mutation at the 1014 locus could, at least in part, be responsible for the patterns of LD around the kdr locus in the Gabonese data. Indeed, there is evidence for repeated mutations of kdr alleles across the species range of Anopheles gambiae [START_REF] Pinto | Origins and distribution of knockdown resistance mutations in the afrotropical mosquito vector Anopheles gambiae[END_REF]).

However, we argue that one the recent time-scales on which kdr has arisen and spread that it is more parsimonious to assume that recombination is the dominant influence on patterns of LD rather than postulating high rates of repetitive mutation.

Whilst kdr is the best-documented resistance mechanism in An. gambiae there are many other resistance-associated loci. Microarray and recombinant protein expression work has shown that resistant mosquitoes over express a small number of enzymes that catalyze insecticide degradation [START_REF] Ortelli | Heterologous expression of four glutathione transferase genes genetically linked to a major insecticide-resistance locus from the malaria vector Anopheles gambiae[END_REF][START_REF] Müller | Transcription profiling of a recently colonised pyrethroid resistant Anopheles gambiae strain from Ghana[END_REF][START_REF] Chiu | Comparative molecular modeling of Anopheles gambiae CYP6Z1, a mosquito P450 capable of metabolizing DDT[END_REF]Müller et al. 2008b). Linkage disequilibrium-based screens could be a powerful way of identifying regions of the genome carrying the scars of recent selection that regulate such overexpression. However, whether association mapping approaches will effectively identify genes subject to much older and comparatively weaker selection is currently unclear. The bounded estimate of the selection coefficient reported here is at the upper limit of estimates generated to date and of a similar magnitude to estimates generated for resistance alleles in the mosquito Culex pipiens [START_REF] Labbe | Resistance gene replacement in the mosquito Culex pipiens: fitness estimation from long-term cline series[END_REF]. In human populations mutations associated with resistance to malaria infection such as G6PD and sickle cell trait have coefficients of selection of 0.02-0.05 [START_REF] Tishkoff | Genetic analysis of African populations: human evolution and complex disease[END_REF] and 0.05-0.18 [START_REF] Li | The first arrival time and mean age of a deleterious mutant gene in a finite population[END_REF] respectively. In the third actor in the malaria transmission cycle of Plasmodium falciparum a selection coefficient of 0.1 has been obtained for the locus dhfr that confers resistance to the chemotherapeutic agent, pyrimethamine [START_REF] Nair | A selective sweep driven by pyrimethamine treatment in southeast asian malaria parasites[END_REF].

Together with strong and recent positive selection, the major determinant of LD around selected loci will be the rate of recombination. Indications of dramatic variation in the recombination rate across the An. gambiae genome have already been reported [START_REF] Pombi | Variation in recombination rate across the X chromosome of Anopheles gambiae[END_REF][START_REF] Black | The yin and yang of linkage disequilibrium: Mapping of genes and nucleotides conferring insecticide resistance in insect disease vectors[END_REF] and it is possible that, being close to the centromere of chromosome 2L the sodium channel locus is in an area of reduced recombination. However, our Kenyan data are consistent with rates of recombination sufficient to reduce the region hitch-hiked with a selectively advantageous locus in a relatively short period of time. Indeed detection of the signatures of selection for loci with low selection coefficients will be more logistically challenging in An. gambiae than humans because of much lower background levels of LD (Weetman D, Wilding CS, Steen K, Donnelly MJ unpublished). We attempted to amplify microsatellites from around the sodium channel to fully define the extent of the swept region as has been done for drug resistance loci in Plasmodium falciparum [START_REF] Wootton | Genetic diversity and chloroquine selective sweeps in Plasmodium falciparum[END_REF][START_REF] Nair | A selective sweep driven by pyrimethamine treatment in southeast asian malaria parasites[END_REF]. However, the sodium channel is situated in a region with an abundance of repetitive sequences and it was not possible to identify unique locus-specific microsatellite primer pairs. Given the apparently high selection pressure on the L1014F mutation, it is curious that there are no studies, with adequate sample size, that have observed either of the kdr alleles at fixation [START_REF] Santolamazza | Distribution of knock-down resistance mutations in Anopheles gambiae molecular forms in west and west-central Africa[END_REF]). One explanation would be that of over-dominance; however, insecticide bioassays studies suggest this is unlikely to be the case [START_REF] Chandre | Modifications of pyrethroid effects associated with kdr mutation in Anopheles gambiae[END_REF][START_REF] Reimer | Relationship between kdr mutation and resistance to pyrethroid and DDT insecticides in natural populations of Anopheles gambiae[END_REF]) and our estimate of the dominance coefficient shows the kdr L1014F allele to be partially recessive. Therefore it is likely that there is some fitness cost to the L1014F allele, and that this could be attributable to heterogeneity in exposure to pyrethroids in the environment or a consequence of an Hill-Robertson effect where selection One generation per month is assumed following [START_REF] Lehmann | The effective population size of Anopheles gambiae in Kenya: implications for population structure[END_REF]). 95% confidence intervals for each observed data point were calculated according to [START_REF] Newcombe | Two-sided confidence intervals for the single proportion: comparison of seven methods[END_REF].

Expected data generated from simultaneous maximum-likelihood estimates of initial frequency and selection and dominance coefficients (see methods).

Figure legends FIG 1 .

 legends1 Figure legends FIG 1. Bifurcation plots showing patterns of recombination in the centromeric (5' towards the left) and telomeric (3' towards the right) directions. The core is marked by the dark circle and each of the 29 SNPs is represented by a node and a recombination event is represented by a bifurcation. The diameter of the circle at each SNP node is proportional to the numbers of individuals with the same long range haplotype at that position. No bifurcation plot is shown for the L1014F core in Ghanaian M-form populations as only a single haplotype was observed (see results)

Table legends Table 1 .

 legends1 Origin and kdr genotype of specimens used in the study.

Table 2 .

 2 Comparison of median EHH levels between alleles at the kdr loci.

						Number of each kdr genotype
			Total
	Population	Year		Form	L1014S/	L1014F/	L1014F/
			N					L1014S/wt L1014F/wt Wt/wt
		Collected			L1014S	L1014F	L1014S
	Asembo Bay, Kenya, 00°10'S,	2005 1
			48	S	11	-	-	17	-	20
	34°22'E
	Dienga, Gabon, 01°52'S, 12°40'E	'99-2000 2	30	S	-	-	-	4	2	24
	Bakoumba, Gabon, 01º49'S, 13º01'E '99-2000 2	42	S	-	5	8	5	7	17
	Libreville, Gabon, 00°22'N, 09°26'E '99-2000 2	73	S	34	8	31	-	-	-
	Okyereko & Accra area, Ghana, 05°	2002 3
			35	S	-	33	-	-	2	-
	24.9N 00° 36.6W, 05°38'N, 00°15'E
	Okyereko, Ghana, 05° 24.9N 00°	2002 3
			30	M	-	-	-	-	2	28
	36.6W

Acknowledgements.

We thank Nelson Cuamba, Eveline Klinkenberg and Pie Muller for providing wild caught specimens and the associate editor for constructive comments.

Funding.

This work was supported by the National Institutes for Health [U01 AI58542 to Dr Edward Walker, Michigan State University] the Innovative Vector Control Consortium and the Royal Society.

at a kdr locus can interfere with the selection at nearby beneficial mutations [START_REF] Hill | The effect of linkage on limits to artificial selection[END_REF].

The data presented herein show that it is possible to detect genomic signatures of strong positive selection in pest species with large effective population size and generally low levels of linkage disequilibrium. We suggest that such approaches are likely to extremely powerful in many non-model taxa subject to similar selective events, Supplementary Materials.

Supplementary tables are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).

The population name and total numbers of each DNA sample utilized. Molecular form is indicated and the numbers of each kdr genotype shown.

Additional information on the collection sites may be obtained from the publications where the specimens are originally described 1 (Müller et


A.-L Fevre-Gautier

A Beylot

M.-P Boutin

P Michel

L Rouvreau

S Vaxelaire

J Villeneuve

CONSTRUCTION AND DEMOLITION WASTE: A SYSTEMIC ANALYSIS TO IMPROVE RECYCLING

. The use of secondary raw materials, thanks to recycling, has to be strengthened to decrease pressure on ecosystems caused by natural resources extraction and waste elimination. The ASURET project aims at identifying and analyzing barriers to C&D waste recycling, focusing on technical and non-technical factors including waste sorting and characterization, recycling and demolition processes and organizational or human aspects regarding the whole sector and its economic actors. Through a systemic approach, the project draws the basis of a long-term thinking, considering territories as evolutionary deposits of materials. An objective is to provide the C&D actors with a specific approach and decision-aid system, dedicated to resources and waste flows optimization in order to lower environmental impacts of C&D activities. This research builds on a multi-level territorial approach, ranking from very local (the building) to national scale, considering the question of the geographic and temporal scale at which recycling and valorization networks must be thought in order to gather the best economical, technical and social conditions for material flows optimization. The work relies on an analysis of the production, availability and exchanges of C&D waste data, at different territorial scales, at the economic actors' level and throughout the whole responsibility chain of a C&D project. These data are used to conduct a territorial material flow analysis for stock estimation and input / output material flows accounting applied to C&D activities. This methodology is being applied on two French case-studies: Orléans City and the Bouches-du-Rhône department. Expected final results will include the design of a decisionaid system to be tested with the two selected case studies. It will be based on recommendations for a multi-level territorial data collection system and multi criteria analysis, including the results of an economy wide life cycle analysis (LCA) based on the input / output LCA methodology and focused on the use of primary and secondary materials for construction.

The results presented here derive from the work conducted for the ongoing ASURET project, funded by the French National Research Agency from February 2010 to March 2013 and led in partnership with 13 Développement, the French Scientific and Technical Centre for Building (CSTB), INSAVALOR POLDEN and University of Technology of Troyes (UTT).


Xiao Yang
email: xiaoyang@broadinstitute.org

Florian Sikora
email: sikora@univ-mlv.fr

Guillaume Blin
email: gblin@univ-mlv.fr

Sylvie Hamel
email: hamelsyl@iro.umontreal.ca

Romeo Rizzi
email: romeo.rizzi@uniud.it

Srinivas Aluru
email: aluru@iastate.edu

An Algorithmic View on Multi-related-segments: a unifying model for approximate common interval

A set of genes that are proximately located on multiple chromosomes often implies their origin from the same ancestral genomic segment or their involvement in the same biological process. Among the numerous studies devoted to model and infer these gene sets, the recently introduced approximate common interval (ACI) models capture gene loss events in addition to the gene insertion, duplication and inversion events already incorporated by earlier models. However, the computational tractability of the corresponding problems remains open in most of the cases. In this contribution, we propose an algorithmic study of a unifying model for ACI, namely Multi-related-segments, and demonstrate that capturing gene losses induces intractability in many cases.

Introduction

The genetic blueprint of an organism is encoded in a set of DNA sequences, known as chromosomes. During evolution, some subsequences of a chromosome diverged while others, known as genes, were conserved among different organisms. A chromosome is typically represented as a sequence of genes, then evolution is described as a series of discrete events: gene insertion, loss, duplication and inversion. One of the most important goals in comparative genomics is to identify a set of genes that are in proximate locations on multiple chromosomes and their actual chromosomal occurrences. Indeed, preservation of gene co-locality tends to indicate that the corresponding genes either form a functional unit (e.g., operons) or result from speciation or duplication events [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF]. In the literature, the former is termed "gene cluster " [START_REF] Bergeron | [END_REF], whereas the latter is known as "synteny" [START_REF] Yang | [END_REF]. Both were extensively studied during the past decade, and numerous models and algorithms were proposed to define and identify them. From an algorithmic point of view, we present a unified model to capture approximate common intervals and provide tractability results in association with evolutionary events.

Gene Proximity: Properties and Models

Modeling gene proximity based on biological intuition is known to be difficult, but some key properties have been raised by Hoberman and Durand [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF]. We present a formalization of these properties by developing the notion of Multi-related-segments [START_REF] Yang | A Unified Model for Multi-genome Synteny and Gene Cluster Inference[END_REF][START_REF] Yang | An improved model for gene cluster inference[END_REF], meanwhile, show that some of them are inadequately captured by existing models. We consider here related algorithmic aspects.

Key Properties of Gene Proximity

Observing the co-occurrence of a gene set A in different chromosomal segments indicates the common origin of these segments. Genes in A are referred as ancestral genes. Naturally, these segments of interest are subject to evolutionary constraints. The first crucial constraint consists in evidence of any gene of interest as being ancestral. This property is usually related to observing a minimum β occurrences of such a gene among the segments, thereby reducing the possibility of misinterpreting what is in fact a chance occurrence. Secondly, each segment contributes sufficiently to the ancestral gene set. More formally, each segment contains at least m different ancestral genes. Then, consider gene loss and insertion events that may have occurred, such an segment may not necessarily contain all ancestral genes while each may pick up genes independently. To constrain the frequency of these events so that the signal of common origin can still be detected, local and global ancestral gene density constraints apply. The former is captured by allowing at most α interleaving genes between two consecutive ancestral genes, while the latter is captured by allowing a maximum l gene losses in each segment and a maximum t total gene losses among all segments.

Existing Models

Consider k chromosomes, each represented as a permutation over a given gene set A. A conserved segment [START_REF] Nadeau | Lengths of chromosomal segments conserved since divergence of man and mouse[END_REF] consists of a set of genes that occur consecutively in the same order on every input chromosome. Once the constraint of the preserved ordering is removed, it leads to the common interval (CI) model definition [START_REF] Uno | Fast algorithms to enumerate all common intervals of two permutations[END_REF]. If the unordered pair of the first and the last genes of a CI is the same on each chromosome, this CI is moreover called conserved [START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF]. Furthermore, if we relax the constraint that genes in a CI have to be consecutive in each chromosomal occurrence -namely, two genes belonging to a CI can be interleaved by a bounded number of genes not belonging to it -the definition of gene-teams (GT) model [START_REF] Bergeron | The algorithmic of gene teams[END_REF] follows. The GT model is of higher biological relevance since it in addition captures gene insertions. The aforementioned models can be applied to strings to account for gene duplications, but the number of resulting gene sets complying the model may increase exponentially. More recently, approximate common interval (ACI) models were introduced [START_REF] Amir | Improved approximate common interval[END_REF][START_REF] Rahmann | Integer Linear Programs for Discovering Approximate Gene Clusters[END_REF][START_REF] Böcker | Computation of median gene clusters[END_REF][START_REF] Jahn | Efficient computation of approximate gene clusters based on reference occurrences[END_REF], where not all ancestral genes need to occur in every segment. Among these, Median gene cluster (Mgc) model [START_REF] Böcker | Computation of median gene clusters[END_REF] is the most recent formulation, but the complexity of this model remains open.

Multi-related segments model

We now present Multi-related-segments (Mrs) model [START_REF] Yang | A Unified Model for Multi-genome Synteny and Gene Cluster Inference[END_REF][START_REF] Yang | An improved model for gene cluster inference[END_REF], which is defined as consisting of a set of segments of interest, each evolved from an ancestral segment with gene set A via gene insertion, loss, duplication, and inversion events. Formally, a Mrs is defined as follows. To ensure evidence of being ancestral genes, each gene in A occurs in at least β (≥ 2) segments. Each segment of interest has to contain at least m different ancestral genes and is maximal (i.e., not extendable by including surrounding genes) -thus, imposing a constraint on the minimum contribution to A. Similar to the GT model, the local ancestral gene density is constrained allowing at most α non-ancestral genes between any two consecutive ancestral genes. To control global ancestral gene density, we require each segment to induce no more than l gene losses and the total number of gene losses of all segments to be lower than t . Then, given a set of chromosomes and parameters α, β, m , l and t , the general problem is to identify all Mrs.

The formulation of Mrs captures existing models and holds a better biological intuition. Mrs corresponds to a CI when β = k, m = |A| and α = 0, and to a GT when α ≥ 0. Compared with these two models, Mrs further captures gene loss events. Note that this aspect was already considered in the Mgc model [START_REF] Böcker | Computation of median gene clusters[END_REF]. Nevertheless, there are several major differences. Firstly, Mrs captures the same origin of more than two segments in the absence of strong pairwise similarity information, such as differential gene loss [START_REF] Simillion | Recent developments in computational approaches for uncovering genomic homology[END_REF] and uber-operon [START_REF] Che | Detecting uber-operons in prokaryotic genomes[END_REF] -which is not the case for Mgc due to the requirement that segments pairwisely share some common genes. Secondly, the minimum evidence of a gene being ancestral is more flexible in Mrs by requiring β occurrences of any ancestral gene -which has to be at least k 2 in Mgc. Finally, the local ancestral gene density is not required in Mgc -which is, as explained in [START_REF] Hoberman | The incompatible desiderata of gene cluster properties[END_REF], crucial.

From an algorithmic point of view, regarding all above mentioned models, complexity increases when chromosomes are delineated as strings rather than permutations: the problem is still tractable when considering CI [START_REF] Bergeron | Computing Common Intervals of K Permutations, with Applications to Modular Decomposition of Graphs[END_REF][START_REF] Didier | Character sets of strings[END_REF] but folds into the hardness as soon as conserved CI is considered [START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF]. The GT model, which captures gene insertions, duplications and inversions, is polynomial on permutations [START_REF] Bergeron | The algorithmic of gene teams[END_REF] but exponential over strings [START_REF] Pasek | Identification of genomic features using microsyntenies of domains: domain teams[END_REF]. Considering the complexity of ACI models, which further captures gene losses, an algorithm with O(kn 3 + occ) run time over strings was proposed [START_REF] Amir | Improved approximate common interval[END_REF] where occ is bounded by the number of substrings of the genomes. In this paper, we further investigate the complexity of ACI models, by considering from an algorithmic point of view the problem of Mrs inference. Since known algorithmic results are available in previous modeling of gene duplications, insertions and inversions, our focus is on deriving if the problem is tractable when trying to model gene losses.

Complexity Analysis of MRS

One has to note that, in the following, we will set the numerical parameters of the model to specific values. This just consists in an algorithmic trick for ease of proof. We first consider the case when the ancestral gene set A is a priori known. The problem, termed LocateMrs, then corresponds to locate, given k chromosomes S = {S 1 , S 2 , . . . , S k } represented as strings, a feasible Mrs originating from A. We prove that this problem is NP-hard even in the restricted case where |CS(S i )| = |S i | and meanwhile, at most one substring per S i can belong to the resulting Mrs, for every S i ∈ S. It follows that LocateMrs problem is NP-hard. Next, we prove LocateMrs to be fixed-parameter tractable (FPT) [START_REF] Downey | Parameterized Complexity[END_REF], and provide an efficient dynamic programming solution. Then, we prove that the optimization problem to identify all Mrs is hard to approximate. Finally, we show that with the removal of the maximum number of gene loss constraint and the maximum number of substrings per input sequence constraint, a polynomial algorithm can be derived. Due to space constraint, some proofs are deferred to the full version of the paper.

Identify A MRS Given A

Let us consider that no gene insertion is allowed (α = 0), β ≥ 2 and t = l = ∞. Then, by definition, any Mrs consists of substrings involving only genes in A. Thus, each input chromosome can be pre-processed in order to remove any gene not belonging to A, resulting in a sequence of substrings. One may then filter out any substring that does not respect the minimum contribution to A criterion (i.e., using m ). Any remaining substring will be referred as of interest. Finally, since in the Mrs definition, we are looking for maximal substrings (i.e., not extendable by including surrounding genes), any substring of interest will be either kept or fully removed in the solution. string of interest of S j (i.e., j th chromosome), find a subsets S ⊆ S corresponding to a Mrs, such that A = S∈S CS(S) and each character of A appears in at least two elements of S , and ∀S a i , S b j ∈ S , i = j.

Definition 1. LocateMrs: Given a character set A, a string set S = {S 1 1 , S 2 1 , . . . , S 1 2 , S 2 2 , . . . , S 1 k , S 2 k . . .} where S i j corresponds to the i th sub- ST              T1 = u1x1u2x2u3x3 -v1w1v2w2v3w3z3 T2 = v2x2v3x3u4x4 -u2w2u3w3z3v4w4z4 T3 = w3x3v4x4u5x5 -u3v3z3u4w4z4v5w5z5 T4 = w4x4v5x5u6x6 -u4v4z4u5w5z5v6w6 T5 = v1x1w5x5v6x6 -u1w1u5v5z5u6w6 SX                      S1,1 = x1 S2,1 = x2 S6,1 = x6 S1,2 = u1 -w1 S2,2 = u2 -w2 S6,2 = u6 -w6 S1,3 = u1 -v1 S2,3 = u2 -v2 S6,3 = u6 -v6 S1,4 = w1 -v1 S2,4 = w2 -v2 S6,4 = w6 -v6 S3,1 = x3 S4,1 = x4 S5,1 = x5 S3,2 = u3 -v3 S4,2 = u4 -v4 S5,2 = u5 -v5 S3,3 = w3 -v3 S4,3 = w4 -v4 S5,3 = w5 -v5
We will prove that LocateMrs is hard but fixed-parameter tractable. We first consider that |S| = |CS(S)| for any S ∈ S (i.e., S is a permutation). Note that this problem is in NP. Indeed, given a subset S of S, one can check in polynomial time that each character of A appears in at least two elements of S and that no more than one substring S i j of any S j belongs to S . To prove that this problem is moreover NP-hard, we provide a polynomial reduction from the NP-complete problem Exact Cover by 3-sets (X3C) [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-completeness[END_REF]: Given a finite set X = {x 1 , . . . , x |X | } and a family T = {t 1 , . . . , t |T | } of triples over X , is there a subfamily T ⊆ T such that every x i ∈ X is contained in exactly one element of T ?

X3C problem is hard even in the special case where each element of X appears at most three times in T [START_REF] Garey | Computers and Intractability: A guide to the theory of NP-completeness[END_REF]. Then, it is sufficient to consider the case where each element appears either two or three times. Indeed, any triple containing some element that occurs only once has to be part of any solution and can be removed from further consideration. According to the problem definition, a solution corresponds to a selection of one among the at most three occurrences of any element of X . Without loss of generality, we fix the triple order in T .

Let us now provide the construction from any instance (X , T ) of X3C problem (an example is given in Figure 1). For each element x i ∈ X , let x i , u i , v i , w i and z i be some characters. The set S will be built on |T | sequences, which represent the triples of T , and four (resp. three) additional sequences, which represent any element of X occurring twice (resp. three times) in T . Let S T = {T 1 , . . . T |T | } (resp. S X ) be the set of sequences representing the triples (resp. the elements of X ). Moreover, we use the symbol "-" to separate the non-adjacent substrings in a given string, e.g., S = S 1 -S 2 -S 3 . Note that, the order of the characters in these substrings is not important according to the definition of Mrs. Let us first construct S T as follows. First, for each element x i ∈ X occurring twice in T , concatenate u i x i (resp. v i x i ) to T 1 j (initially empty) and v i w i (resp. u i w i ) to T 2 j (initially empty) if the first (resp. second) occurrence of x i appears in the j th triple of T . Second, for each element x i ∈ X occurring three times in T , concatenate u i x i (resp. v i x i and w i x i ) to T 1 j and v i w i z i (resp. u i w i z i and u i v i z i ) to T 2 j if the first (resp. second and third) occurrence of x i appears in the j th triple of T . Let us now construct the set S X . For each element x i ∈ X occurring two times in T , add the following four sequences to S X :

S i,1 = x i , S i,2 = u i -w i , S i,3 = u i -v i , S i,4 = w i -v i .
And, for each element x i ∈ X occurring three times in T , add the following three sequences to S X : S i,1 = x i , S i,2 = u i -v i , S i,3 = w i -v i . We finally define A to be the set of all characters used in the construction. Lemma 1. There exists a solution T ⊆ T to X3C problem over (X , T ) if and only if in the corresponding built instance (A, S) of LocateMrs there exists a subset S ⊆ S corresponding to a Mrs.

Correctness of Lemma 1 implies the following result.

Theorem 1. LocateMrs problem is NP-complete even in the special case where none of the input strings contains duplicated characters and at most one substring S i j of every S j can belong to any solution S . We now prove that LocateMrs belongs to the class of the fixedparameter tractable (FPT) problems [START_REF] Downey | Parameterized Complexity[END_REF]. In other words, it can be solved efficiently by an algorithm exponential only with respect to a fixed parameter -|A| in our case -while polynomial in the size of the input.

Theorem 2. LocateMrs problem is Fixed-Parameter Tractable in |A|

To show this, we provide a dynamic programming solution. According to LocateMrs definition, one has to select exactly one substring of interest among all of them in each sequence S j . A naive algorithm may try all such combinations and check for each if any character appears in at least two substrings. Such an algorithm has an exponential running time. We will prove that by using an efficient dynamic programming strategy, one may hold the exponential factor in the size of the ancestral gene set. Note that one does not need to compute the exact number of times each character occurs but only to ensure that it occurs in at least two substrings in the solution. According to this remark, consider a fixed ordering of characters (a 1 , a 2 , . . . , a |A| ) of A, we compute after adding substring S to the current solution a vector C = (c 1 , c 2 , . . . , c |A| ), where c i ∈ {0, 1, 2} denotes respectively that a i is not contained, contained in one, or contained in at least two substrings. For example, consider A = {1, 2, 3, 4, 5} and current solution S = {124}, one may derive a vector C = (2, 1, 0, 2, 1) after adding substring "1445" to S . The main property of this representation is that, given A, there are only 3 |A| possible vectors. Further, let µ(x) and χ S (x) denote, respectively, the position of x in the fixed ordering of A and the boolean function indicating whether x occurs in S. We define a boolean dynamic table D indexed by the last substring added and the vector C for the current solution. The main recursion1 is defined as follows:

D(S i j , (c 1 , . . . c |A| )) =          1 -if ∃i , j < j s.t. D(S i j , (c 1 , . . . c |A| )) = 1 and ∀1 ≤ l ≤ |A|, χ S i j (x) + c l = min{2, c l } where µ(x) = l 0 -otherwise Algorithm 1 LocateMrs(A, S = {S 1 1 , S 2 1 , . . . , S 1 2 , S 2 2 , . . . , S 1 k , S 2 k , . . .}) 1: Initialize all entries of D to 0 2: for each S i 1 ∈ S do D(S i 1 , (c1, . . . , c |A| )) = 1 where ∀x ∈ S i 1 , c µ(x) = χ S i 1 (x) done 3: for j = 2 to k do 4:
for each S i j ∈ S do Fill out D(S i j , (c1, . . . , c |A| )) done 5: end for 6: for each S i k ∈ S do 7:

if D(S i k , (2, . . . , 2)) = 1 then return True end if 8: end for 9: return False Given this function, one can apply Algorithm 1. The algorithm computes, for each sequence S j the possible character set solution induced by any combination of substrings of interest from sequences S j with j < j. Therefore, any entry D(S i j , (2, . . . , 2)) = 1 corresponds to a Mrs being found. One may, using a simple back-tracking technic, rebuild one optimal solution. Let us now prove the time complexity of this algorithm. In order to fill out D, one has to compute |S| × 3 |A| entries. Indeed, there are at most |S| different substrings and 3 |A| possible character sets. The main recursion needs, for each entry, to browse at most |S| × 3 |A| other entries of D. This leads to an overall O((|S| × 3 |A| ) 2 ) running time algorithm. Hence, the problem is FPT with respect to |A|.

Identify All MRS When A Is Unknown

We will prove that finding all Mrs problem is hard even in the special case where none of the sequences contains duplicated characters and in any solution S , for any sequence S j at most one substring S i j ∈ S (i.e., α = 0, β ≥ 2, t = l = ∞).

First, note that the problem is in NP since given a subset S of S, one can polynomially check that each element of A appears in at least two substrings and no more than one substring of any sequence belongs to S . To prove that this problem is moreover NP-hard, we provide a polynomial reduction from the NP-complete problem X3C [START_REF] Yang | An improved model for gene cluster inference[END_REF] based on a slight modification of the reduction of the previous section. Indeed, if one replaces each of the separations "-" between substrings of interest by a unique character appearing only once in S, then by definition, those added characters will never be part of a Mrs since any character should appear at least twice in a Mrs. Due to the unextandability property of Mrs, one should be able to find neither a smaller nor a bigger substring of interest in each sequence than in the LocateMrs formulation. The rest of the proof still holds, leading to the following theorem.

Theorem 3. Finding a Mrs problem is NP-complete even in the special case where none of the sequences contains duplicated characters and in any solution S , at most one substring from each S j belongs to S .

Let us then consider the optimization version of the problem (Definition 2) where one wants to find a Mrs induced by the maximum unknown ancestral gene set (in other words, one constrains the minimum size of A), and at the same time, at most one substring of each S j can belong to the Mrs. Definition 2. MaxMrs: Given a set of k strings S = {S 1 , . . . S k }, find any possible (A, S ) where S = {S 1 , S 2 , . . . S k : S i is a substring of S i }, A = k i=1 CS(S i ), and |A| is maximum. We will demonstrate that this optimization problem is hard to approximate. Meanwhile, we show that the inapproximability of this problem may stem from forbidding more than one substring per input chromosome, the relaxation of which leads to polynomiality.

In the following, we consider that β ≥ 2, α = 0, m = 1, and t = l = ∞. We prove the inapproximability of MaxMrs below by proposing a reduction from the Minimum Set Cover (MinSC) problem: Given a family F of subsets of a finite universe U, find a set cover F for U -that is a subfamily F ⊆ F whose union is U -of the minimum cardinality. Since any character appearing once in an input string will not be part of a Mrs, we use the symbol "-" to denote any such character. The presence of symbol "-" will induce, in MaxMrs problem, that characters appearing before and after it in any input string cannot be part of the same solution. Given any instance (F, U) of MinSC, where U = {u 1 , . . . u n } and F = {F i : Proof. We first define a polynomial-time subroutine that transforms any solution (A, S ) to an equally good solution where CS(S 0 ) ⊆ A. For any character of S 0 not belonging to A -say u j , add to S one substring S 1 i that was not in S but contains u j , meanwhile, remove from S correspondingly two substrings S 2 i and S m+i . Every such replacement operation will change a given v i by u j in A without decreasing the cardinality of A (i.e. an equally good solution). Once this subroutine has been applied to (A, S ), one can build a cover F = {F i : S 2 i ∈ S }. The subroutine guarantees that all elements of S 0 belong to F -a cover for U. Clearly, |A| -n corresponds to the number of v j s belonging to A. Considering Proof. Consider Minimum 3-SetCover-3 (Min3SC-3), a subproblem of MinSC, where the size of any set in F is bounded by 3 as well as the number of times each character of U occurs in F. We will prove the theorem by contradiction, assuming that MaxMrs admits a Polynomial-Time Approximation Scheme (PTAS), i.e. one would be able to find an approximation algorithm leading to an approximate solution (A AP X , S AP X ), which compared with the optimal solution (A OP T , S OP T ), induces |A AP X | ≥ (1 -) • |A OP T | for a parameter > 0. Accordingly, under the same assumption, we will prove that Min3SC-3 also admits a PTAS, i.e. one would be able to find an approximation algorithm leading to an approximate solution F AP X , which compared with the optimal solution F OP T , induces |F AP X | ≤ (1 + γ) • |F OP T | for a parameter γ > 0 -a contradiction to the fact that Min3SC-3 is APX-hard [START_REF] Papadimitriou | Optimization, approximation and complexity classes[END_REF].

F i = {u 1 i , u 2 i , . . . u n i i }, 1 ≤ i ≤ m}, we define a set of strings S = {S 0 , . . . , S k=2m } with S 0 = u 1 . . . u n , S i = S 1 i -S 2 i = u 1 i u 2 i . . . u n i i - v i and S m+i = v i , for 1 ≤ i ≤ m.
Since each character of U occurs at most three times in F, the size of the ground set used to build F is at most 3n, leading to m ≤ 3n. Moreover, any cover F ⊆ F of U is at least of size n 3 since F is composed of sets of size at most three. Hence, We now prove the following result on restricted instances.

Theorem 5. If one restricts neither the maximum number of gene losses per substring of interest, nor the maximum number of substrings of interest per chromosome, and if every input sequence contains no duplicated characters, finding all the Mrs becomes a polynomial task. Proof. Consider a graph G = (V, E) obtained from S in such a way that a vertex is assigned to every character in each string S i ∈ S and a red (resp. blue) colored edge is created between any two adjacent characters (resp. any two vertices representing identical characters). Given this representation, the notion of character set naturally extends to any subgraph G[V ] of G as the set of represented characters by V . Our method consists in an iterative procedure which stops when none of the following operations can be applied anymore. The results will consist of a set of connected components, each corresponding to a Mrs. The first operation consists in removing from V any vertex which is only incident to red colored edges. This polynomial operation results in the removal of genes not appearing twice in a candidate connected component. The second operation gets rid of candidates not fulfilling the minimum contribution to the ancestral gene set by pruning any red edge-induced subgraph G such that |CS(G )| < m . This operation can be done in linear time by browsing any connected component. Once none of these operations can be done anymore, it is easy to see that each remaining connected component corresponds to a Mrs.

Fig. 1 .

 1 Fig.1. Illustration of the construction on the following instance of X3C: X = {1, 2, 3, 4, 5, 6} and T = {(1, 2, 3), (2, 3, 4),[START_REF] Bergeron | [END_REF][START_REF] Bergeron | The algorithmic of gene teams[END_REF][START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF],[START_REF] Bergeron | The algorithmic of gene teams[END_REF][START_REF] Bergeron | On the similarity of sets of permutations and its applications to genome comparison[END_REF][START_REF] Blin | Comparing genomes with duplications: a computational complexity point of view[END_REF], (1, 5, 6)}. A correspondance between the solutions of the problems is highlighted in bold.

Lemma 2 .

 2 If there exists a cover F ⊆ F for U (i.e. U = F ∈F F ) then there exists a solution (A, S ) (i.e. a Mrs) for the built up instance S of MaxMrs such that |A| = n + m -|F |. Lemma 3. Given a solution (A, S ) for a built up instance S of MaxMrs, we can construct in polynomial-time a cover F ⊆ F for U, such that |F | ≤ m -(|A| -n).

S 1 ,

 1 S 2 ...S m (where S 2 j s appear), there exist at most m -(|A| -n) strings such that S 2 j ∈ S ; inducing that |F | ≤ m -(|A| -n). Theorem 4. MaxMrs is APX-hard even in the special case where, for every input string S i , |CS(S i )| = |S i |.

n 3 ≤

 3 |F OP T | and consequently, m ≤ 9 • |F OP T |.If we have an approximate solution (A AP X , S AP X ), then By Lemma 3,|F AP X | ≤ m -(|A AP X | -n) By assumption, m -(|A AP X | -n) ≤ m -((1 -) • |A OP T | -n) By Lemma 2, |A OP T | = n + m -|F OP T | Which leads to, |F AP X | ≤ • n + • m + (1 -) • |F OP T | (m ≤ 3n ≤ 9|F OP T |) ≤ 12 • |F OP T | + (1 -) • |F OP T | Finally, ≤ (1 + 11 ) • |F OP T |

The base case is made in Algorithm 1

Appendix

Missing proofs for main results

Lemma 1. There exists a solution T ⊆ T to X3C problem over (X , T ) if and only if in the corresponding built instance (A, S) of LocateMrs there exists a subset S ⊆ S corresponding to a Mrs.

Proof. (⇒) Suppose such T ⊆ T exists for the instance (X , T ). Let us prove that we can compute in polynomial time an S ⊆ S corresponding to a Mrs in the corresponding built instance of LocateMrs. For any triple t q ∈ T , if t q ∈ T then add T 1 q to S ; add T 2 q otherwise. Moreover, for 1 ≤ i ≤ |X |, add S i,1 to S . For any element x i ∈ X occuring two times in T , if the first triple containing x i is in T then add S 2 i,2 , S 2 i,3 , S 2 i,4 to S ; add S 1 i,2 , S 1 i,3 , S 1 i,4 to S otherwise. For any element x i ∈ X occuring three times in T , if the first triple containing x i is in T then add S 2 i,2 , S 1 i,3 to S ; if the second triple containing x i is in T then add S 1 i,2 , S 1 i,3 to S ; add S 1 i,2 , S 2 i,3 to S otherwise (see Figure 1). Let us now prove that S is indeed a Mrs, i.e., A = S∈S CS(S) and any element of A appears at least twice in S . First note that, by definition, in T , there is exactly one triple containing any element of X . Thus, by construction, there is exactly one substring per sequence in S . Let us consider any element x i ∈ X which occurs twice in T -say in triples t q and t q for the first and respectively second occurrences of x i . Then, since T is an exact cover, there exists exactly one among t q and t q in T -say t q . Therefore, by construction, {T 1 q , T 2 q , S i,1 , S 2 i,2 , S 2 i,3 , S 2 i,4 } ⊆ S (the reader may, for example, consider x 1 in Figure 1). Consequently, a) two occurrences of x i appear in S (one from S i,1 and one from T 1 q ), b) two occurrences of u i appear in S (one from T 1 q and one from T 2 q ), c) two occurrences of v i appear in S (one from S i,3 and one from S i,4 ), and d) two occurrences of w i appear in S (one from S i,2 and one from T 2 q ). With a similar reasoning, one can check that any of x i , u i , v i , w i appears twice when t q ∈ T . Let us now consider each element x i ∈ X that occurs three times in T -say in triples t q , t q and t q for the first (resp. second and third) occurrence of x i . Then, since T is an exact cover, there exists exactly one among t q , t q and t q in Tsay t q . Therefore, by construction, {T 1 q , T 2 q , T 2 q , S i,1 , S 2 i,2 , S 1 i,3 } ⊆ S (the reader may, for example, consider x 3 in Figure 1). Consequently, a) two occurrences of x i appear in S (one from S i,1 and one from T 1 q ), b) three occurrences of u i appear in S (one from T 1 q , one from T 2 q and one from T 2 q ), c) two occurrences of v i appear in S (one from T 2 q and one from S i,2 ), d) two occurrences of w i appear in S (one from S i,3 and one from T 2 q ), and e) two occurrences of z i appear in S (one from T 2 q and one from T 2 q ). With a similar reasoning, one can check that any of x i , u i , v i , w i , z i appears at least twice when t q ∈ T or t q ∈ T . We have completed the proof that each element of A appears at least twice in S ; inducing that S is indeed a Mrs.

(⇐) Suppose now that such a set S ⊆ S exists for the corresponding built instance of LocateMrs (A, S). We will prove that we can compute in polynomial time a T ⊆ T corresponding to a solution for the instance (X , T ). For any

Let us now prove that T is indeed an exact cover of T . Note that, any solution S is a subset of S since we are looking for unextandable substrings (inducing that any substring of interest is either fully kept or removed). Moreover, recall that we consider here the special case of the problem where in any solution S , for any string S j , at most one substring S i j ∈ S . Let us first prove that given two triples t q and t q both containing one of the two occurrences of x i ∈ X (i.e. the case where x i appears twice in X ), then, exactly one of the substrings T 1 q and T 1 q belongs to S (the reader may, for example, consider x 1 in Figure 1). By contradiction, suppose this is not the case, i.e. {T 1 q , T 1 q } ⊆ S . Then, in S , there are already one occurrence of u i , one occurrence of v i and two occurrences of x i . Note that the set of elements related to variable x i , i.e., {x i , u i , v i , w i }, appears in the following sequences {T q , T q , S i,1 , S i,2 , S i,3 , S i,4 }. Since, by definition, in S , every element should appear at least twice, using exactly one of the substrings of each of these sequences, one should be able to obtain another u i , another v i and two w i . Unfortunately, this is not possible in any combination of the corresponding substrings excluding, by hypothesis, T 2 q and T 2 q ; a contradiction. Now consider the three triples t q , t q and t q , each containing one of the three occurrences of x i ∈ X (i.e. the case where x i appears three times in X ), then, exactly one among the substrings T 1 q , T 1 q and T 1 q belongs to S (the reader may, for example, consider x 3 in Figure 1). By construction, the occurrences of z i belong to T 2 q , T 2 q and T 2 q . In order to obtain at least two occurrences of z i in S , at least two of the substrings among {T 2 q , T 2 q , T 2 q } should be in S . Moreover, since the occurrences of x i belong to {T 1 q , T 1 q , T 1 q , S i,1 }, exactly one of {T 1 q , T 1 q , T 1 q } should be in S . We have proved that for each element x i occuring twice (resp. three times) in T , exactly one of the triples containing x i is kept in T . Thus, T is indeed an exact cover of T .

Lemma 2. If there exists a cover F ⊆ F for U (i.e., U = F ∈F F ) then there exists a solution (A, S ) (i.e., a Mrs) for the built up instance S of MaxMrs such that |A| = n + m -|F |.

Proof. Consider the following solution S = {S 0 } {S 1 i ∈ S : F i ∈ F } {S 2 i , S m+i ∈ S : F i ∈ F }. By definition, S is indeed a Mrs since it cannot be extended (it is made of one unextandable substring of each sequence) and any element appears twice. Moreover, its character set A = U {v i : F i ∈ F } which is of size n + m -|F |.