Combining genetic crosses and pool targeted DNA-seq for untangling genomic variations associated with resistance to multiple insecticides in the dengue vector Aedes aegypti

In addition to combating vector-borne diseases, studying the adaptation of mosquitoes to insecticides provides a remarkable example of evolution-in-action driving the selection of complex phenotypes. Indeed, most resistant mosquito populations show multi-resistance phenotypes as a consequence of the variety of insecticides employed and of the complexity of selected resistance mechanisms. Such complexity makes challenging the identification of alleles conferring resistance to specific insecticides and prevents the development of molecular assays to track them in the field. Here we showed that combining simple genetic crosses with pool targeted DNA-seq can enhance the specificity of resistance allele’s detection while maintaining experimental work and sequencing effort at reasonable levels. A multi-resistant population of the mosquito Aedes aegypti was exposed to three distinct insecticides (deltamethrin, bendiocarb and fenitrothion) and survivors to each insecticide were crossed with a susceptible strain to generate 3 distinct lines. F2 individuals from each line were then segregated with 2 insecticide doses. Bioassays supported the improved segregation of resistance alleles between lines. Hundreds of genes covering all detoxifying enzymes and insecticide targets together with more than 7,000 intergenic regions equally spread over mosquito genome were sequenced from pools of F0 and F2 individuals unexposed or surviving insecticide. Differential coverage analysis identified 39 detoxification enzymes showing an increased gene copy number in association with resistance. Combining an allele frequency filtering approach with a Bayesian FST-based genome scan identified multiple genomic regions showing strong selection signatures together with 50 non-synonymous variations associated with resistance. This study provides a simple and cost-effective approach to improve the segregation of resistant alleles in multi-resistant populations while reducing false positives frequently arising when comparing populations showing divergent genetic backgrounds. The identification of these insecticide resistance markers paves the way for the design of novel DNA-based resistance tracking assays. Author summary In addition to combating vector-borne diseases, understanding how mosquitoes adapt to insecticides provides a remarkable example of evolution-in-action. However, the variety of insecticides used and the complexity of adaptive mechanisms make it difficult to identify the genetic changes conferring resistance to each insecticide. Here we combined simple controlled crosses with mass DNA sequencing for enhancing the specificity of resistance gene detection. A multi-resistant mosquito population was exposed to three distinct insecticides and survivors were crossed with a susceptible strain to generate 3 distinct mosquito lines. Individuals from the second generation of each line were then segregated based on their resistance to each insecticide. Bioassays supported the improved segregation of genetic resistance markers between lines. Hundreds of genes potentially involved in resistance together with thousands non-genic regions equally spread over mosquito genome were sequenced from individuals from each line. Genomic analyses identified detoxification enzymes showing an increased gene copy number in association with resistance and multiple genomic regions showing strong selection signatures and carrying point mutations associated with resistance. Such approach improves the specificity of resistance gene detection in field mosquito populations resisting to multiple insecticides and paves the way for the design of novel DNA-based resistance tracking tools.


Author summary
In addition to combating vector-borne diseases, understanding how mosquitoes adapt to 31 insecticides provides a remarkable example of evolution-in-action. However, the variety of 32 insecticides used and the complexity of adaptive mechanisms make it difficult to identify the 33 genetic changes conferring resistance to each insecticide. Here we combined simple 34 controlled crosses with mass DNA sequencing for enhancing the specificity of resistance 35 gene detection. A multi-resistant mosquito population was exposed to three distinct 36 insecticides and survivors were crossed with a susceptible strain to generate 3 distinct 37 and affects all insecticide families used in public health [11], often leading to reduced vector 73 control efficacy [12][13][14]. Although attempts are made to develop alternative arbovirus control 74 strategies [15] their large scale implementation will require decades. Until this, characterizing 75 molecular mechanisms underlying resistance is crucial for tracking down resistance alleles 76 and improving resistance management strategies [16]. 77 Resistance of mosquitoes to chemical insecticides can be the consequence of various 78 mechanisms, such as non-synonymous mutations affecting the protein targeted by 79 insecticides, a lower insecticide penetration, its sequestration, or its biodegradation often 80 called metabolic resistance [6,17]. In Ae. aegypti, resistance to pyrethroids, the main 81 insecticide family used against mosquitoes, is mainly the consequence of target-site 82 mutations affecting the voltage-gated sodium channel targeted by these insecticides (Knock 83 Down Resistance 'kdr' mutations) and metabolic mechanisms [11,18]. Several kdr mutations 84 have been identified in this species and the causal association between the V410L, S989P, 85 V1016G/I and F1534C mutations and pyrethroid resistance has been confirmed [19][20][21][22][23]. 86 Most of these mutations can be genotyped on individual mosquitoes by PCR-based assays, 87 providing essential allele frequency data for resistance management. Conversely, metabolic 88 resistance is far less understood in Ae. aegypti although this type of resistance is frequent 89 and often accounts for a significant part of the resistance phenotype [6]. Such resistance 90 mechanism is caused by an increased activity of detoxification enzymes. These 91 detoxification enzymes include cytochrome P450 monooxygenases (P450s or CYPs for 92 genes), carboxy/cholinesterases (CCEs), glutathione S-transferases (GSTs) and UDP-93 glycosyl-transferases (UDPGTs) although other families can be involved [17,18,24]. Their 94 high diversity (~ 300 genes in Ae. aegypti) and the complexity of biodegradation pathways 95 make challenging the identification of those conferring resistance to a specific insecticide. 96 Theoretically, metabolic resistance can be the consequence of an increased expression of 97 one or multiple detoxification enzymes metabolizing the insecticide and/or the selection of 98 variants showing a higher insecticide metabolism rate due to conformal modifications. As 99 over expression is frequently associated with over transcription, most candidate genes were 100 identified based on their differential transcription in resistant populations as compared to 101 susceptible counterparts using transcriptomics [11,18,24,25]. Although these approaches 102 identified several detoxification enzymes involved in insecticide biodegradation, they mostly 103 failed to pinpoint the underlying genomic changes, thus impairing the high-throughput  However, this study used natural resistant populations displaying multi-resistance 114 phenotypes, thus not allowing to properly discriminating between alleles specifically 115 associated with resistance to the insecticide in question and those associated with resistance 116 to other insecticides. Furthermore, this approach did not allow breaking up the genetic 117 linkages between the genomic variations identified, thus potentially leading to false positives. 118 In this context, the present study aimed at better understanding the origin of complex 119 insecticide resistance phenotypes in the mosquito Ae. aegypti. More precisely, we combined 120 genetic crosses and targeted DNA-seq in an attempt to identify genomic variations 121 specifically associated with resistance to distinct insecticides in a multi-resistant Ae. aegypti 122 population. After exposure to three insecticides of distinct chemical families (the pyrethroid 123 deltamethrin, the organophosphate fenitrothion and the carbamate bendiocarb), survivors to 124 each insecticide were crossed with a susceptible strain to generate three F2 lines. Each F2 125 line was then phenotyped with two increasing doses of its respective insecticide and 126 survivors were used to identify CNV and polymorphism variations associated with resistance 127 in hundreds of target genes including all detoxification enzymes and insecticide target 128 proteins. In addition, the capture of thousands of intergenic regions regularly distributed over 129 mosquito genome also allowed crossing up these data with a genome-wide screening of 130 selection signatures associated with resistance. Overall, this study contributes to improve our 131 understanding of the complex genomic bases of metabolic resistance to insecticides and 132 paves the way for the design of novel insecticide resistance tracking tools in this major 133 arbovirus vector. Bioassays performed on the initial composite population (F0 Guy-R) confirmed its high 140 resistance to the pyrethroid insecticide deltamethrin with resistance ratio (RR 50 ) over 316-fold 141 as compared to the susceptible strain Bora-Bora (Fig. 1A). This population also showed 142 moderate resistance to the carbamate bendiocarb and the organophosphate fenitrothion with 143 RR 50 of 14-fold and 3-fold respectively. As expected, resistance to each insecticide 144 decreased after crossing F0 survivors to each insecticide with the susceptible strain with F1 145 resistance ratios decreasing to 25-fold, 7-fold and 2-fold for deltamethrin, bendiocarb and 146 fenitrothion respectively. Deltamethrin resistance was even lower in F2 (10-fold) while 147 fenitrothion resistance remains low (1.8-fold) and bendiocarb resistance slightly increased to 148 10-fold. Assessing the cross resistance of each line to all insecticides confirmed the partial 149 segregation of resistance alleles after controlled crosses (Fig. 1B). Although F2 individuals 150 from each line showed a higher survival when exposed to its respective insecticide, this trend 151 was only significant for deltamethrin in link with the lower resistance to other insecticides. 152 153

Target-site resistance mutations
Assessing kdr mutations frequencies from targeted DNA-seq reads data confirmed the high 155 frequency of the three kdr mutations V410L, V1016I and F1534C in F0 Guy-R composite 156 population corroborating its high deltamethrin resistance level (Fig. 2). Exposing F0 Guy-R 157 individuals to the LD 80 of each insecticide did not segregate these mutations in survivors. The 158 segregation of kdr mutations became more evident in F2 individuals with higher allele 159 frequencies observed in F2 individuals of the Delt line surviving to high dose of deltamethrin. 160 As expected because genetically constrained by two successive mutation events in Ae. were observed when the number of genotyped mosquitoes was low (S1 Fig.).  conditions and passed quality and coverage filters (S2 Table). These variations were mostly  Table). Most of them were located in proximity of genes potentially involved in 211 metabolic resistance and included genes carrying non-synonymous variations associated 212 with resistance (see below). Among these regions, 5 were located on chromosome 1, 213 including two GST and 1 P450 clusters. The P450 cluster (3 CYP304 genes at ~287 Mb) 214 using field mosquito populations focused on resistance mechanisms to a given insecticide or 261 a given chemical family. However, these kinds of studies did not fully discriminate alleles 262 associated with resistance to different insecticides, which may lead to false positives. In this 263 context, the present study attempted at demonstrating that the identification of resistance 264 alleles can be improved by combining simple genetic crosses and targeted DNA pool 265 sequencing, while maintaining experimental work and sequencing costs at reasonable levels. As compared to deltamethrin, more genes were affected by CNV associated with bendiocarb 328 and fenitrothion resistance. Indeed, even though resistance levels to these two insecticides 329 were lower, the absence of resistance mutations affecting the target of these insecticides in 330 Ae. aegypti because of genetic constraints [  Although none of these mutations were found associated with resistance in our study, some 391 mutations identified in these CCEs are located near the catalytic triad (e.g. I330M in 392 CCEae3A and D332G in AEL005123) or the active site (e.g. P293A in AAEL019678). 393 Altogether, the present study allowed identifying multiple detoxification genes located in dose-response bioassays with deltamethrin, fenithrothion and bendiocarb were performed on 465 the F0 Guy-R composite population to assess its resistance phenotype and identify the LD 80 466 of each insecticide to be used for the initial F 0 segregation. These bioassays were conducted 467 with at least five doses of deltamethrin (0.05 % to 1%), fenitrothion (0.0125% to 0.4%) and Resistance levels of F1 and F2 individuals from each line were obtained following the same 474

procedure. 475
Cross-resistance profiles of F2 individuals from each line to all insecticides were then 476 evaluated using single-dose bioassays. For each insecticide, the dose was calibrated in 477 order to obtain a mortality ranging from 20% to 40% in the corresponding F2 line. Doses 478 used were as follows: deltamethrin 0.05% for 60 min, bendiocarb 0.5% for 60 min, order to cover > 95% of Ae. aegypti genome using the following criteria: target region size = 504 220 bp; optimal distance between 2 regions = 150 kb ± 10 kb; region distance to any 505 annotated gene > 5 kb; avoid repeated and redundant regions; avoid regions with GC 506 richness > 70% or single nucleotide richness > 50%; avoid regions with undefined 507 nucleotides (N); do not consider supercontigs < 150 kb; avoid regions located within 75.5 kb 508 of supercontig boundaries. All genomic regions targeted by the study are detailed in S4 509 Library' protocol vB.4. Briefly, 3 µg of genomic DNA from each sample were fragmented 512 using a Bioruptor (Diagenode), purified, ligated to adaptors and amplified by PCR using 513 Herculase II DNA polymerase (Agilent Technologies). After QC of library size and quantity, 514 libraries were hybridized to biotinylated baits and purified using Dynal MyOne streptavidin 515 beads (Invitrogene). Captured DNA fragments were amplified, purified and multiplexed 516 before sequencing. Sequencing was performed on an Illumina NextSeq500 and generated 517 more than 300 million 75 bp paired reads with an average of 23.3 million reads per sample. 518 Reads were assigned to each sample (unplexing) and adaptors were removed. Reads 519 quality was checked for each sample using FastQC 520 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and reads were loaded into 521 Strand NGS v3.1.1 (Strand Life Science) for further analyses. 522 Reads mapping and filtering. In order to minimize false positives arising from mapping bias 523 in high-redundancy and low complexity regions, CNV were only identified from coding 524 regions. Reads were mapped against all Aaeg L5 exons using a padding of 35 bp and the 525 following parameters: minimum 90% identity, maximum 5% gap, mean insert size of 167 ± 526 30 bp, mismatch penalty = 4, gap opening penalty = 6, gap extension penalty = 1, clipping 527 penalty = 5, min align read length = 30, ignore reads with more than 5 matches, trim 3' end if 528 base quality < 25. Reads were then filtered on their quality metrics and mapping quality using 529 the following criteria: mean read quality ≥ 28, N allowed ≤ 2, alignment score ≥ 90, Mapping 530 quality ≥ 40, read length ≥ 35, remove non primary multiply mapped reads, remove inter-531 chromosomal split reads. Finally, mate missing, translocated and duplicated reads were 532 removed. For polymorphisms analysis, reads were mapped against the whole Aaeg L5 were expected to increase from F0 to F0 survivors, decrease from F0 survivors to F2 after 546 crossing with the susceptible strain, and increase in F2 survivors. No dose response 547 condition was applied for CNV in order to minimize the confounding effect of target site 548

mutations. 549
Polymorphisms and selection signatures. Variants were call against the whole Aaeg L5 550 genome using the following parameters: coverage > 30 in all conditions, confidence score cut 551 off = 100, ignore loci with homopolymer stretch > 4, ignore loci with average base quality ≤ 552 15, ignore loci with strand bias ≥ 50 and coverage ≥ 50, ignore reads with mapping quality ≤ 553 20, ignore variants with less than 4% supporting reads. Among all variants called only those 554 polymorphic among our conditions (i.e. showing ≥ 5 % variation in at least two conditions) 555 were retained and their genic effects were computed. 556 Associations between polymorphisms and resistance to each insecticide were assessed by 557 combining an allele frequency filtering approach with an F ST -based approach. 558 The frequency filtering approach was based on the expected resistance allele frequency 559 variations across conditions taking into account their initial frequency. Frequency thresholds 560 used are described in detail in Table 1. Basically, the frequency of alleles positively 561 associated with resistance was expected to increase from unexposed F0 individuals to F0 562 survivors, decrease from F0 survivors to unexposed F2 individuals (following crossing with 563 the susceptible strain), and increase again in F2 survivors in association with the insecticide 564 dose. Different initial allele frequency thresholds were used for identifying alleles associated 565 with deltamethrin (≥30% in F0 Guy-R) and those associated with bendiocarb and fenitrothion 566 resistance (≥15% in F0 Guy-R) to account for the higher deltamethrin resistance level of the 567 initial composite population. The frequency of deleterious alleles (i.e. those negatively 568 associated with resistance) was expected to behave reciprocally. 569 The F ST -based approach aimed at assessing departure from neutrality using the Bayesian

21.
Yanola J, Somboon P, Walton C, Nachaiwieng W, Somwang P, Prapanthadara LA. 675 High-throughput assays for detection of the F1534C mutation in the voltage-gated sodium 676 channel gene in permethrin-resistant Aedes aegypti and the distribution of this mutation 677 throughout Thailand.