Standing genetic variation in laboratory populations of insecticide‐susceptible Phlebotomus papatasi and Lutzomyia longipalpis (Diptera: Psychodidae: Phlebotominae) for the evolution of resistance

Abstract Insecticides can exert strong selection on insect pest species, including those that vector diseases, and have led to rapid evolution of resistance. Despite such rapid evolution, relatively little is known about standing genetic variation for resistance in insecticide‐susceptible populations of many species. To help fill this knowledge gap, we generated genotyping‐by‐sequencing data from insecticide‐susceptible Phlebotomus papatasi and Lutzomyia longipalpis sand flies that survived or died from a sub‐diagnostic exposure to either permethrin or malathion using a modified version of the Centers for Disease Control and Prevention bottle bioassay. Multi‐locus genome‐wide association mapping methods were used to quantify standing genetic variation for insecticide resistance in these populations and to identify specific alleles associated with insecticide survival. For each insecticide treatment, we estimated the proportion of the variation in survival explained by the genetic data (i.e., “chip” heritability) and the number and contribution of individual loci with measurable effects. For all treatments, survival to an insecticide exposure was heritable with a polygenic architecture. Both P. papatasi and L. longipalpis had alleles for survival that resided within many genes throughout their genomes. The implications for resistance conferred by many alleles, as well as inferences made about the utility of laboratory insecticide resistance association studies compared to field observations, are discussed.


| INTRODUC TI ON
When populations occupy stressful environments or are faced with novel challenges, extinction can result unless populations rapidly adapt (Bell, & Gonzalez, 2009;Gomulkiewicz, & Holt, 1995;Orr, & Unckless, 2014). Such evolutionary rescue can prevent species loss from global change, but can also hamper attempts to eradicate pests and pathogens. Populations can adapt to such challenges from standing genetic variation or from new mutations, both of which can have various levels of genetic complexity from single genes of large effect to many genes with individually smaller effects (Barrett & Schluter, 2008;Chevinm, Lande, et al., 2010;Chevin, Martin, et al., 2010;Orr, 2005). Standing genetic variation is expected to contribute disproportionately to the early stages of rapid adaptation, especially in multicellular Eukaryotes (e.g Barrett & Schluter, 2008;Burke et al., 2010;Messer, & Petrov, 2013;Rêgo et al., 2019). As such, a better understanding of standing genetic variation for survival and fecundity in novel, stressful environments is critical for predicting a specie's response and persistence under such conditions. Human activities have caused novel selective pressures and resulted in species decline (Bell & Collins, 2008). One of the best examples of human-induced selective pressures is insecticide exposure and the evolution of insecticide resistance. Since the 1940s, synthetic insecticides have been used successfully to reduce disease transmission from arthropod vectors, but today resistance has been documented in the most important vectors including mosquitoes, black flies, triatomines, lice, fleas, and sand flies to the insecticides that are used in vector control campaigns (Alexander & Maroli, 2003;Hemingway & Ranson, 2000;Palumbi, 2001;Rivero et al., 2010). The continued application of insecticides, both in designed application and indiscriminately causing repeated exposures, and subsequent evolution of resistance have dampened enthusiasm for disease eradication in favor of sustained control (Hemingway et al., 2002;Nauen, 2007). Despite this, insecticides remain a reliable tool for controlling vectors, but their utility is eroding as fewer insecticides remain viable for control (Hemingway & Ranson, 2000;Hemingway et al., 2016;Kelly-Hope et al., 2008). New classes of insecticides have only recently been developed and are being evaluated for efficacy and implementation (Agumba et al., 2019). Integrated approaches are needed to possibly move away from total reliance on insecticides (Kelly-Hope et al., 2008).
Phlebotomine sand flies (Diptera: Psychodidae) are vectors that transmit Leishmania protozoans that cause leishmaniasis in humans, a disfiguring, stigmatizing, and lethal disease which causes tens of thousands of deaths each year worldwide (Alvar et al., 2012;Hotez, 2008;World Health Organization [WHO], 2010;WHO, 2013). Only females in the genera Phlebotomus and Lutzomyia are the competent, putative vectors of these parasites (Akhoundi et al., 2016).
Phlebotomus papatasi (Scopoli) and Lutzomyia longipalpis (Lutz and Neiva) are two important vectors of Leishmania protozoans to people in the Old World and New World, respectively (Belo et al., 2013;Maroli et al., 2013). Sand flies, including P. papatasi and L. longipalpis, remain for the most part, susceptible to insecticides (Coleman et al., 2011). There is, though, increasing evidence of insecticide resistance in the Middle East, southern Asia, and South America (Alexander et al., 2009;Dinesh et al., 2010;Faraj et al., 2012;Hassan et al., 2012Hassan et al., , 2015Karakus et al., 2017;Khan et al., 2015;Saeidi et al., 2012;Surendran et al., 2005). Despite the recent findings of resistance in sand fly populations around the world, there is little knowledge about the genetic and molecular mechanisms of resistance in these populations. An understanding of these mechanisms will be crucial for the success of sand fly control programs to reduce the leishmaniasis burden without exacerbating resistance. Vector control programs based on known mechanisms of insecticide resistance in sand fly populations will have a starting point to make informed, effective control decisions about using alternative insecticides or using other integrated control methods (Alexander et al., 2009;Alexander & Maroli, 2003;Faraj et al., 2012;Surendran et al., 2005).
Conventional insecticide resistance testing often focuses primarily on the mechanisms of target-site insensitivity and metabolic detoxification (ffrench-Constant et al., 2004;Hemingway et al., 2004;Nauen, 2007). However, resistance is likely more complicated.
Many genes with different mechanisms can collectively contribute to the resistance phenotype Vontas et al., 2005Vontas et al., , 2007. For example, whole-genome sequencing also has revealed high complexity of copy number variation at insecticide resistance loci in malaria mosquitoes . More robust methods are now needed to scan the sand fly genome for genetic markers associated with insecticide exposure survival.
The goal of this study is to quantify standing genetic variation for survival following insecticide exposure in laboratory populations of insecticide-susceptible P. papatasi and L. longipalpis. To that end, we used genotype-by-sequencing (GBS) and multi-locus genome-wide association methods to quantify standing genetic variation for resistance to two insecticides (malathion and permethrin) and identify genetic loci associated with insecticide resistance (Comeault et al., 2014(Comeault et al., , 2015Romay et al., 2013). While such methods result in only a modest density of genetic markers relative to whole-genome sequencing, they provide a cost-effective approach to sequence a sufficient number of individuals for genetic mapping feasible in nonmodel systems. We discuss the strengths and limitations of such approaches for mapping in more detail in light of our specific results in the discussion.
Following insecticide exposure, all sand flies were captured via mechanical aspiration and released into 1-pint cardboard containers with a mesh top onto which a cotton ball saturated with 30% sugar-water was placed and served as an energy/water source. The containers were kept in the same growth chamber as the insecticidesusceptible colonies. Sand flies were held in these containers for 24-h when mortality was observed as a complete cessation of movement (Denlinger et al., 2015;Perea et al., 2009). Here, phenotypes for insecticide resistance were assigned as a binary score of surviving (1) or dying from exposure (0). For P. papatasi, permethrin exposure resulted in 64.6% survival (n = 128 alive, n = 64 dead) while malathion led to 23.4% survival (n = 45 alive, n = 147 dead; Table 1).

| DNA sequence alignment and variant calling
Custom perl scripts were used to first demultiplex pooled DNA sequences, wherein identifier barcodes served to assign DNA sequences to individual sand flies . Reference We used the "aln" algorithm from "bwa" (version 0.7.5; ) to align the DNA sequences to the L. longipalpis or P. papatasi reference genome (Table 1). We allowed for a maximum of four nucleotide differences, no more than two mismatches in a 20 bp seed, and a quality threshold for read trimming set to 10. Along with these parameters, only reads with a single best match were aligned. A small number of sand flies with few aligned sequences were removed before subsequent analyses. Sequence coverage was 6× for P. papatasi compared to 16× for L. longipalpis, consistent with the larger genome size of P. papatasi than L. longipalpis (Table 1).

TA B L E 1
Summary statistics for the Phlebotomus papatasi and Lutzomyia longipalpis permethrin and malathion treatments: insecticide exposure survival, number of individuals yielding DNA, number of DNA sequences obtained and the percent that aligned to the reference P. papatasi and L. longipalpis genomes from VectorBase, the average genome coverage, the number of variants called for each species (combining both treatments), and the minor allele frequency (MAF) correlation between survivors and dead individuals in each treatment Using "samtools" and "bcftools" (version 0.1.19; , sequence alignments were sorted and indexed for variant calling.
Treatment groups of each sand fly species were combined for this process. As recommended for Illumina HiSeq data, coefficients to cap mapping quality and the number of reads per position were set to 50.
Bases with a quality score below 15 and reads with a mapping quality below 10 were ignored. The prior for θ was set to 0.001, and only single nucleotide variants (SNVs) where the posterior probability of an invariant nucleotide was below 0.01 were retained (Li, 2011). Each SNV set was filtered based on a 128-read minimum for overall coverage, a four-read minimum for the non-reference allele, a minimum mapping quality of 30, and a maximum of 20% of individuals with missing data.
Our filtering criteria were selected to ensure sufficient coverage to estimate allele frequencies while accounting for uncertainty in genotype (Buerkle & Gompert, 2013), and while also avoiding locus drop-in and drop-out (that is, where mutations create or remove DNA sequence motifs cut by the restriction enzymes). Final sets of 38,657 and 18,856 SNVs (~1 SNV per 10 kpb for each species) were retained for P. papatasi and L. longipalpis, respectively (Table 1).

| Estimating genotypes, allele frequencies, and linkage disequilibrium
We estimated allele frequencies for each species and insecticide treatment. Maximum likelihood allele frequency estimates were obtained using an expectation-maximization algorithm that accounts for uncertainty in genotypes Li, 2011). Relative to methods that rely on first calling genotypes, this approach has the advantage of allowing for the inclusion of individuals with a range of sequence coverage and weighting their contributions to the allele frequency estimates by the information carried in their sequence data (Buerkle & Gompert, 2013).
Genotype estimates are required for association mapping. Thus, we next used a Bayesian approach to estimate genotypes for each SNP and individual. Our empirical Bayesian approach uses the allele frequency estimates to define prior probabilities for genotypes, such g denotes the counts of, for example, the non-reference allele (0, 1 or 2 in diploids) and p denotes the corresponding allele frequency.
Posterior probabilities were then obtained according to Bayes rule as Pr(g| D, p) = [Pr(D|g) Pr(g)]/Pr(D), where Pr(D|g) defines the likelihood of the genotype given the sequence data and quality scores as calculated by samtools and bcftools. We then obtained point estimates (posterior means) of genotypes as Pr(g = 0|D,p)*0 + Pr(g = 1| D,p)*1 + Pr(g = 2|D,p)*2. This results in genotype estimates that take on values between 0 and 2 (copies of the non-reference allele) but that are not constrained to be integer valued).
Pairwise linkage disequilibrium (LD) was calculated in each species from our genotype estimates using the "geno-r2" function "vcftools" (version 0.1.15; Danecek et al., 2011). Specifically, we measured LD as the squared correlation between genotypes at pairs of SNPs and computed LD for all pairs of SNPs in 100 kb windows.

| Genome-wide association mapping
Binary Bayesian sparse linear mixed models (BSLMMs) were fit with "gemma" (version 0.98; Zhou et al., 2013) to estimate genetic contributions to variation in insecticide survival and to identify SNVs with this phenotypic variation. Here, survival outcomes were modeled as a function of a polygenic term (denoted u) and a vector of the potential measurable SNV effects (denoted β) (Zhou et al., 2013). A Markov chain Monte Carlo (MCMC) algorithm with variable selection was used to infer posterior inclusion probabilities (PIPs) for SNVs with a non-zero measurable effect on insecticide susceptibility, and modelaverage point estimates (MAPEs) were derived from those PIPs (Zhou et al., 2013). The polygenic term in each BSLMM represents expected deviations from a phenotypic mean based on all SNVs while accounting for phenotypic covariance that arise between sand flies due to relatedness or genetic similarity (i.e., observed kinship; Zhou et al., 2013).
Relatedness was also considered when estimating individual SNV effects (β) and their PIPs with kinship matrices.
Parameters for estimating genetic architecture were derived from the hierarchical structure of the BSLMM (Guan & Stephens, 2011;Lucas et al., 2018;Zhou et al., 2013). Altogether, the parameters indicate the proportion of the phenotypic variance explained (PVE) by additive genetic effects (based on β and the polygenic term), the proportion of PVE explained by measurable-effect SNVs (PGE) or those implicated by LD (β alone), and the number of SNVs with effects that explain phenotypic variance (nγ).
Thirty independent MCMC chains were run for binary BSLMMs, wherein a probit link function was used to connect the binary response (survival outcome) to a latent quantitative risk variable.
MCMC chains included 100,000 burn-in steps, 1 million sampling steps, and a thinning interval of 10. We assessed convergence to the posterior distribution by calculating the Gelman-Rubin potential scale reduction diagnostic for PVE, PGE and nγ in R with the "CODA" package (version 0.19.3; Plummer et al., 2006;R Core Team, 2013); values of this statistic for were generally less than 1.1 consistent with convergence. To reduce bias in estimation, inferences were carried out using the combined values from all iterations across chains (Cowles & Carlin, 1996).

| Insecticide survival predictions
We used five-fold cross-validation to evaluate the predictive power of the genome-wide association mapping models. To do this, we refit the BSLMM model five times for each data set (species and insecticide treatment). In each case, we used a random 80% of the observations as a training set to fit the model and the other 20% to evaluate the model.
We fit the BSLMM models via MCMC with 100,000 steps as a burnin, followed by 1 million sampling steps with a thinning interval of 10.
The fit model was used to predict the survival phenotype of the test individuals, that is to obtain genomic-estimated breeding values for each of the test individuals based on the additive effects of genes were captured by both β and u in the BSLMMs Lucas et al., 2018). We used the full set of predictions across the five-fold cross-validation sets to assess predictive performance. This was done using the R package "ROCR" (version 1.0.7; Sing et al., 2005); receiver operator characteristic (ROC) curves were constructed to interpret the area under the curve (AUC) and determine the predictive power in correctly classifying survival outcomes.

| Variant effect predictions
We used the Ensembl Variation Effect Predictor on VectorBase to characterize the genomic context and possible consequences of each SNV in the data set, that is to classify SNV based on their effect if in exons (e.g., synonymous, missense, etc.) or genomic context if not (e.g., intron, 3′ UTR, 5′ UTR, intergenic, etc.) (Giraldo-Calderón et al., 2015;McLaren et al., 2010McLaren et al., , 2016. We then summarized the annotations for the 100 SNVs most associated with survival in each treatment for each species and used randomization tests (1000 randomizations each) to determine whether any category was overrepresented relative to null expectations.

| Genetic variation
As expected, allele frequencies were highly correlated between surviving and dead sand flies for each species and treatment (Table 1, Figure S1). Average allele frequency differences (i.e., the mean, absolute difference in the frequency of each allele) between surviving and dead flies were 0.042 (malathion) and 0.033 (permethrin) in L. longipalpis and ~0.025 (both treatments) in P. papatasi ( Figure 1).
Nonetheless, change for some SNVs was much higher, with maximum values of 0.23-0.32 across species and insecticide treatments.
Also as expected, greater allele frequency differences between surviving and dead flies was seen for SNVs with higher minor allele frequencies (i.e., more genetic variation; Pearson correlations between 0.36 and 0.49, all p < 0.001).
Linkage disequilibrium decayed with physical genomic distances in both P. papatasi and L. longipalpis (Figure 2). Nonetheless, nontrivial LD persisted at a sufficient distance for the SNV markers to likely exhibit LD with at least a reasonable proportion of causal variants. In particular, with a marker density of ~1 SNV per 10 kb, we would expect most causal variants to be within 5 kb of at least one SNV maker. At the scale of 5 kb, mean LD measured by r 2 was 0.021 in P. papatasi (maximum = 1.0) and 0.047 in L. longipalpis (maximum = 0.80).

| Genome-wide association mapping
Point estimates of the proportion of variation in survival explained by additive genetic effects (PVE) ranged from 14.7% for P. papatasi exposed to malathion to 90.1% for L. longipalpis exposed to malathion (Table 2). However, these estimates were associated with considerable uncertainty (Table 2). Moreover, with the exception of P. papatasi exposed to permethrin, we lacked sufficient data for precise estimates of the proportion of the PVE attributable to individual genetic variants with measurable effects (PGE) versus nearinfinitesimal effects. Estimates of the number of causal variants with measurable effects (nγ) were lower in P. papatasi than L. longipalpis for both insecticides (Table 2).
Consistent with this, several SNVs had modest to large PIPs for P. papatasi exposed to permethrin (six SNVs with PIPs LD was present for SNVs with high effects on permethrin survival ( Figure S2).

| Insecticide survival predictions
Standing genetic variation in P. papatasi was moderately sufficient in predicting permethrin survival (AUC = 0.68, which denotes a 36% increase in predictive power relative null expectations of AUC = 0.5; Figure 5a), but was no better than null expectation in terms of predicting malathion survival (AUC = 0.36; Figure 5b). In L. longipalpis, we observed a small but non-zero increase in predictive power relative to a null model for both permethrin (AUC = 0.53; Figure 5c) or malathion (AUC = 0.59, Figure 5d) exposure.

| Variant effect predictions
In both species, most SNVs occurred outside of genes (e.g., VEP categories intergenic, upstream or downstream of genes) (Figure 6a

| D ISCUSS I ON
We found evidence of standing genetic variation for reduced susceptibility to permethrin and malathion in susceptible lab colonies of P. papatasi and L. longipalpis. This suggests a potential for these species to evolve resistance to these insecticides. We discuss our  to sequence a sufficient number of flies to identify some genetic markers associated with resistance and the overall contribution of genetics to the trait (i.e., the PVE). Thus, we think this approach was useful, but some caution is warranted when interpreting our results.
Second, our results apply specifically to these lab colonies. It is simply unclear at this time whether or to what extent the same genetic variants are segregating in nature. Still, given the fact that standing genetic variation is often shared among populations, we think it is quite likely that at least some of these same causal variants are segregating in nature.

| Genetic architecture
In P. papatasi, survival to a sub-diagnostic dose of permethrin is her-  Perhaps the susceptible population of P. papatasi did not already have the genetic variation to survive malathion's different mode of action from permethrin. However, it is important to note that the effectual LC that flies were exposed to was higher for malathion than for permethrin, and thus the sand flies could harbor additional standing genetic variation for survival to lower concentrations of malathion. Posterior inclusions probabilities for the highest ranking SNVs were also much lower, and not surprisingly, our power to predict survival to exposure to malathion was considerably lower too.
Phenotypic variation for survival ability in permethrin-exposed

| Gene associations
Intergenic variants and variants associated with genes were among the top five highest ranking SNVs in all four treatment groups. The variants associated with genes were found in genes or upstream or downstream of them. Some genes do not yet have an annotated function in the sand fly genomes. The genes that are annotated have a diverse range of metabolic and biochemical functions (Tables S1-S4). We must be cautious, though, in our inferences. Despite being able to analyze tens of thousands of variants, only a small portion of the genome is sequenced with GBS. Some of the variants we found associated with survival to an insecticide exposure may be causal; organophosphates, like malathion. They are up-or downregulated in resistant insects (Chambers & Oppenheimer, 2004;Vontas et al., 2007) and are important for synthesis and conformation of detoxifying enzymes in the presence of organophosphates (Ahmed et al., 1998). Zinc fingers (high MAPE scores in the malathion exposure in both sand fly species) are transcriptional repressors (Kasai & Scott, 2001). In Musca domestica, mixed functional oxidase (MFO), a class of insecticide detoxifying enzymes, promoters bind transcription repressor genes that contain zinc finger moieties. The MFO promoters in pyrethroid-resistant M. domestica bind the repressor genes less than in susceptible individuals because of polymorphisms in the repressor gene. This causes increased transcription of MFOs, which are able to detoxify pyrethroid insecticides (Gao & Scott, 2006;Perera et al., 2008). It is possible that the upstream variant of the zinc finger encoding gene contributes to MFO repression. Decreased MFOs can also confer resistance because they first must enzymatically activate insecticide, which they later detoxify. With fewer MFOs, there are fewer bioactivated insecticides (Scott, 1999). Perhaps variants near or within zinc fingers contribute to increased or decreased MFO expression and either can lead to insecticide resistance.
Several SNVs were found that are associated with proteins in the L. longipalpis malathion-exposure treatment (Table S4). A SNV was found to be associated with a protein containing a disulfide isomerase function. GSTs in insects are known to alter isomerase activity (Sheehan et al., 2001). In the same treatment, microtubule associated protein RP/EB were upregulated found in lambdacyhalothrin resistant Aphis glycines. Microtubule associated proteins interact with postsynaptic proteins in the nervous system. They could help stabilize dendrites to normalize nerve function when malathion disrupts synaptic transmission by inhibiting acetylcholinesterase (Lepicard et al., 2014). Intra-flagellar transport proteins were less abundant in imidacloprid-resistant Myzus persicae (Meng et al., 2014). Glycosyltransferases are detoxification enzymes, and overexpression of some uridine diphosphate-glycosyltransferases has been shown to confer resistance in lepidopteran agricultural pests . Lastly, a SNV was found associated with a gene that transcribes a protein with an alpha/beta hydrolase fold activity. Carboxylesterase and cholinesterase enzymes, such as acetylcholinesterase, evolved from a core alpha/beta hydrolase, and these enzymes frequently confer insecticide resistance (Hotelier et al., 2010).

| Standing genetic variation and adaptation
Despite more genetically homogenous laboratory populations of sand flies (Lanzaro et al., 1998;Mukhopadhyay et al., 1997Mukhopadhyay et al., , 1998Mukhopadhyay et al., , 2001, insecticide exposure survival is a known heritable trait and can lead to resistance (Feyereisen, 1995;Hemingway et al., 2002;Rivero et al., 2010). In theory, alleles for survival will increase in frequency toward fixation with continued selection, disseminate throughout the population, and result in greater population survival over the course of continued exposure (Xu et al., 2012). The rate of evolution in a population depends on multiple factors, including the initial allele frequency (Roush & McKenzie, 1987). The insecticidesusceptible colonies used in this experiment were derived from 30year inbred populations that were most likely homozygous for many loci and maybe during that time emergent pre-adaptive alleles were removed through purifying selection and/or through stabilizing selection. Despite evidence of sufficient standing genetic variation for selection to act upon, this variation could have been very little.
Selection for resistance in a laboratory population falls within the phenotypic distribution of the susceptible population, often below the LC 100 for an insecticide (ffrench-Constant et al., 2004;Oakeshott et al., 2013;Roush & McKenzie, 1987). This selection process is conducted to allow survivors for subsequent generations. In doing so, existing, common variation is selected for, which produces polygenic resistance.
Because of the homogeneity of laboratory populations, very low initial frequency of resistance alleles, the high fitness costs of those resistance alleles, and the weakness of the selection process, the evolution of resistance from major-effect alleles is potentially unlikely (Lande, 1983;McKenzie et al., 1992). Even a LC 90 of an insecticide has the potential to produce polygenic resistance (McKenzie & Batterham, 1994).

F I G U R E 5
Predictive power of the genome-wide association models based on receiver operating characteristic (ROC) curves. ROC curves are shown for Phlebotomus papatasi survival when exposed to malathion (area under the ROC curve [AUC] =0.36) and permethrin (AUC = 0.68) (a) and Lutzomyia longipalpis when exposed to malathion (AUC = 0.59) and permethrin ( Our lineages are being exposed to an approximate LC 50 of permethrin and malathion, so it is certainly expected that we will find evidence of polygenic resistance. Monogenic resistance can be successfully selected for in the laboratory if selection concentration is set above the LC 100 of an insecticide (McKenzie & Batterham, 1998). With diagnostic doses for many insecticides for sand flies recently described (Denlinger, Creswell, et al., 2016), selection for major-effect alleles is possible in the future.
Resistance selection in field populations is much greater (above the LC 100 for an insecticide) and can be outside of the phenotypic range of insecticide tolerance. This can result in the rapid selection of rare, major-effect mutations that can lead to monogenic or oligogenic  resistance that present as target-site insensitivity, metabolic detoxification, or both epistatically (Edi et al., 2014;ffrench-Constant et al., 2004;Hardstone et al., 2009;McKenzie & Batterham, 1998;Saavedra-Rodriguez et al., 2008;Whitten et al., 1980). Here, large sizes of field populations act as a source of rare mutations, whereas the small population sizes of inbred individuals in a laboratory population only lead to an accumulation of small effect-size mutations (ffrench-Constant, 2013;McKenzie et al., 1992). It is the heterogeneity of field populations that allows for rare variants to exist (Groeters & Tabashnik, 2000). Interestingly, rare variants may precede the selection for resistance. For example, In Australia, mutations for organophosphate resistance in Lucilia blow flies predated the use of malathion. Examples of standing genetic variation of resistance alleles in field populations, prior to insecticide use, demonstrate that these alleles are under balancing selection and do not carry a high enough fitness cost (ffrench-Constant, 2007). Alleles already present in populations are known to quickly increase in frequency from human-induced evolution (Messer et al., 2016). This may be why resistance has evolved very rapidly when insecticides are first introduced as a control method (Hemingway & Ranson, 2000).
Laboratory strains initiated from field populations with monogenic resistance may not always evolve monogenic resistance because of the factors associated with polygenic resistance selection (Groeters & Tabashnik, 2000;Kasai et al., 2014;Zhu et al., 2013). This may be why Fawaz et al., (2016) did not find target-site insensitivity mutations in their laboratory colony initiated from Egyptian P. papatasi. Even so, resistance in the field may be more polygenic than initially perceived, and this could be due to fitness costs and pleiotropy from major-effect mutations. Microarrays have found many genes with various functions involved in resistance, more than could be found by simply testing for known resistance mechanisms including target-site insensitivity and metabolic detoxification (Djouaka et al., 2008;Pedra et al., 2004;Vontas et al., 2005Vontas et al., , 2007. These findings demonstrate that insecticide resistance, in both the field and laboratory, is a complx phenotype that combines major-effect changes (target-site insensitivity and metabolic detoxification) and many other alleles that are beginning to be discovered and understood. We found that selecting for insecticide exposure survival in laboratory colonies of sand flies is possible but challenging. There is sufficient standing genetic variation in our colonies for polygenic resistance mechanisms. Polygenic resistance is not frequently found in field populations of insects because of greater selection pressure and larger pools of genetic diversity, but it is possible (Groeters & Tabashnik, 2000;Raymond & Marquine, 1994). Polygenic insecticide resistance found in the field is maintained by low mutation rates and minimal migration, both of which are a source of new alleles for monogenic resistance (Raymond & Marquine, 1994;Zhu et al., 2013).

| Resistance control implications
A question that remains is whether polygenic resistance is likely in field populations of sand flies. Sand flies are weak fliers, distribute poorly, and are vagile, which together can lead to small, genetically structured populations (Belen et al., 2011;Doha et al., 1991;Hamarsheh et al., 2007;Khalid et al., 2012;Morrison et al., 1993;Orshan et al., 2016). The weaker effect of selection in smaller populations, and the stronger effect of drift, could dilute resistant alleles should they arise through mutation (Lanfear et al., 2014). Additionally, smaller populations are less likely to be rescued and more likely to go extinct (Willi et al., 2006), but this is not always true (Ferriere & Legendre, 2012). Compound these factors with little gene flow from poor migration, or with gene flow from susceptible sand flies that were unexposed to insecticide due to inadequate insecticide coverage in the environment, and susceptible alleles could remain commonplace in a population. These maladapted alleles, under insecticide selection pressure, could build up a migration load should there be migration (Bolnik & Nosil, 2007). From a control standpoint, these features could be an exploitable opportunity for a failure of evolutionary rescue that may not be seen in other insect vectors. Rapid evolutionary adaptation may not be realistic in these fragmented populations in nature because of potentially little standing genetic variation, and they would be susceptible to stochastic population decline and extinction with the relative inability for adaptation to save them (Gonzalez et al., 2013). Additionally, our findings that the SNVs associated with survival to permethrin and malathion are mostly independent suggests that cross-resistance in sand flies to multiple insecticide classes may require many SNVs and/or mechanisms. Alternative classes of insecticides would remain viable in the presence of resistance, which would be advantageous for sand fly control programs.
For our laboratory populations, predictions, not assumptions and conclusions, should be made about the mechanisms of insecticide resistance in field populations (Mukhopadhyay et al., 1997).
The results from this experiment should serve as a model, not a standard or representative of sand flies in the field. More research of survival and resistance mechanisms using GBS needs to be investigated in natural populations and incorporated into effective integrated vector management programs. GBS's utility in scanning entire genomes of vectors for markers associated with insecticide exposure survival, in both field and laboratory populations, should be incorporated into studies examining the genetic mechanisms of insecticide resistance. GBS will enhance research that examines insecticide use, refuge populations, and gene flow for when insecticide coverage for vectors is uneven, heritability and dominance levels of resistance, fitness costs, and the dynamics of polygenic resistance becoming monogenic resistance (Mallet, 1989;McKenzie et al., 1992;Neve et al., 2009;Tabashnik et al., 2003).

ACK N OWLED G EM ENTS
We are grateful to the many undergraduate research assistants and volunteers in the Bernhardt lab for their assistance with maintaining and rearing the sand fly colonies.