Incongruent Nuclear and Mitochondrial Genetic Structure of New World Screwworm Fly Populations Due to Positive Selection of Mutations Associated with Dimethyl- and Diethyl-Organophosphates Resistance

Livestock production is an important economic activity in Brazil, which has been suffering significant losses due to the impact of parasites. The New World screwworm (NWS) fly, Cochliomyia hominivorax, is an ectoparasite and one of the most important myiasis-causing flies endemic to the Americas. The geographic distribution of NWS has been reduced after the implementation of the Sterile Insect Technique (SIT), being eradicated in North America and part of Central America. In South America, C. hominivorax is controlled by chemical insecticides, although indiscriminate use can cause selection of resistant individuals. Previous studies have associated the Gly137Asp and Trp251Leu mutations in the active site of carboxylesterase E3 to resistance of diethyl and dimethyl-organophosphates insecticides, respectively. Here, we have sequenced a fragment of the carboxylesterase E3 gene (ChαE7), comprising part of intron iII, exon eIII, intron iIII and part of exon eIV, and three mitochondrial gene sequences (CR, COI and COII), of NWS flies from 21 locations in South America. These markers were used for population structure analyses and the ChαE7 gene was also investigated to gain insight into the selective pressures that have shaped its evolution. Analysis of molecular variance (AMOVA) and pairwise FST analysis indicated an increased genetic structure between locations in the ChαE7 compared to the concatenated mitochondrial genes. Discriminant analysis of principal components (DAPC) and spatial analysis of molecular variance (SAMOVA) indicated different degrees of genetic structure for all markers, in agreement with the AMOVA results, but with low correlation to geographic data. The NWS fly is considered a panmitic species based on mitochondrial data, while it is structured into three groups considering the ChαE7 gene. A negative association between the two mutations related to organophosphate resistance and Fay & Wu’s H significant negative values for the exons, suggest that these mutations evolved under positive selection.

Introduction lipids metabolism) and phosphorylated (associated to OP resistance) forms, enabling a better understanding of the effects of both mutations on the wild-type function of the enzyme [25].
The two mutations (Gly137Asp and Trp251Leu) were identified in the NWS fly and it is hypothesized that these confer resistance through a similar mechanism based on homology with mutations in other species for which resistance has been proven through bioassays [18,22,[26][27][28]. Because of the absence of a susceptible lineage, bioassays are difficult to conduct in NWS fly and the detection of positive selection signals in field populations represents an alternative approach to indirectly demonstrate the association of both mutations with insecticide resistance [29,30].
Another important factor to successfully implement a control program based on SIT is the adequate delimitation of the geographic areas to be targeted [31]. Many population genetic studies, based on several molecular markers, were conducted with this aim in the NWS fly, which resulted in different population structure scenarios and interpretations [32][33][34][35][36][37][38][39]. At least four genetic regional groups were described for the NWS fly throughout its current geographical range, including Cuba (CG), the Dominican Republic (DRG), and North and South of the Amazon Region (NAG and SAG, respectively), with island populations being derived from mainland ones [38]. In the mainland, a split between populations from North/Central America and South America, preceded by an expansion process initiated in North America, occurred after the Last Glacial Maximum. A subsequent split between the Pleistocene and Holocene epochs resulted in the NAG and SAG regions in South America [39]. Within these two areas, the NWS fly did not share mitochondrial haplotypes [38] and the genetic structure was associated with a barrier in the north of the Amazon basin, with NWS fly populations from the Amazon only sharing haplotypes with those in the SAG but not in the NAG region [40]. Despite the wide area sampled, a low population differentiation was observed without a geographic correlation in the SAG group [38], probably due to the population expansion process [39].
In this study, we characterized a fragment of carboxylesterase E3 gene (ChαE7) and investigated the geographic distribution of its genetic variability in NWS fly populations from SAG group locations. We also carried out comparative analyses of ChαE7 variability by also investigating the same parameters of mitochondrial genes. Neutrality tests based on frequency spectrum and linkage disequilibrium tests were also used to detect positive selection signals in the mutations associated with insecticide resistance.

Ethics Statement
Samples were obtained in private farms and no specific permissions were required for the present studies. Additionally, we confirm that the field studies did not involve endangered or protected species.

Carboxylesterase E3 (ChαE7) gene characterization
The total genomic DNA of five NWS fly individuals from Brazil (2) and Uruguay (3) was extracted by the phenol:chloroform method with modifications [36]. Polymerase chain reaction (PCR) amplification of a ChαE7 fragment comprising intron iII, exon eIII, intron iIII and part of exon eIV (Fig 1) was carried out initially using the specific primers 7F0a (5´-GGTATACCA TACGCCCAAC -3´) and RN2 (5´-AACAGTAATCCCTCGTACG-3´). Both primers were designed based on NWS fly cDNA [41] using Primer 3 software [42]. PCR amplifications were conducted in 15μL containing 10X PCR buffer, 200 μM of each dNTP, 1 unit of bovine serum albumin (BSA, New England Biolabs, Ipswich, MA, USA), 2.5 μM MgCl 2 , 1 μM of each primer, 1 unit of Taq DNA polymerase (Fermentas International Inc., Canada) and~15-20 ng of DNA. Amplifications were performed with an initial denaturing step of 3 min at 96°C, followed by 35 cycles of 1 minute at 95°C, 1 minute at 62°C and 2 minutes at 72°C, and ending with a final elongation step of 10 minutes at 72°C. The reaction products were sequenced and used for nested PCR with new primers (7FIn1: 5´-ATTGTGTCTCCCTGCAAGTG-3´and 7R1aN: 5 -CGTTTAGTTTCTGGAGCC-3´) based on the obtained sequences. The nested PCRs were conducted as above but with a reduced MgCl 2 concentration (1.5 μM) and an annealing temperature of 52°C. PCR products were purified with the Illustra GFX PCR DNA and Gel Band purification Kit according to manufacturer's instructions (GE Healthcare, Little Chalfont, UK) and cloned into the pGEM-T vector (Promega Co., Madison, WI, USA). Clones were sequenced with forward and reverse M13 primers using the Big Dye Terminator Cycle Sequencing ABI Kit 3.1 in an automatic sequencer ABI 3700 (Applied Biosystems, Inc.; Foster City, CA, USA). Chromatograms were inspected with FinchTV 1.4.0 software (Geospiza Inc., Seattle, WA, USA) taking into account Phred values [43,44]. Contiguous, overlapping bidirectional sequences (contigs) of each individual were obtained using CAP3 software [45], aligned with the algorithms of ClustalX2 [46] and Muscle [47], as implemented in Mega 5.2 software [48], and manually edited. Sequence length and nucleotide composition (A%, T%, C%, G%, A+T%, C+G%, AT-skew = [A-T]/[A+T] and GC-skew = [G-C]/[G+C]) of forward strands were compared with Lucilia cuprina (GenBank accession numbers: U56636 and AY691508), Drosophila melanogaster (GenBank accession number: NM079537) and Musca domestica (GenBank accession number: AF133341).

Population samples and sequencing
Stored genomic DNA (-80°C) from 201 NWS fly larvae collected in 21 sampling locations from South America, previously analyzed in [36] and [38], were used in this study ( Table 1).
Fragments of three mitochondrial genes, the cytochrome c oxidase subunit I (COI) and subunit II (COII) and the B domain of the control region (CR), were amplified and sequenced using the same protocols and conditions for 54 flies that were not previously analyzed by [38]. Part of the ChαE7 fragment previously characterized, including 157 bp from intron iII, the exon eIII (150 bp), the intron iIII (63 bp) and 282 bp from exon eIV, was amplified from all the individuals using the forward primer 7FIn2 (5´-ACCATCGGTGAGTTGAGAG-3´) and reverse primers RN2 (5´-AACAGTAATCCCTCGTACG-3´) or 7R3a (5´-ATCCTTATCATTATTTT CACCC -3´) (Fig 1), leading to a predicted fragment size of 652 bp. The forward primer Exons and introns were named according to L. cuprina (e = exon; I = intron) [27]. The position of primers used for gene region characterization (7F0a, RN2, 7FIn1 and 7R1aN) and for posterior populational analyses (7FIn2, RN2 and 7R3a) are indicated by blue and red arrows, respectively. RN2 is indicated by a green arrow (used for both characterization and population analyses). (7Fln2) was designed based on sequences previously obtained in this study, with Primer 3 software [42]. The reverse primer 7R3a was designed as described previously [28]. PCR amplifications were performed in a final volume of 25μL containing 10X PCR buffer, 100 μM of each dNTP, 1 unit of bovine serum albumin (BSA, New England Biolabs, Ipswich, MA, USA), 2.5 mM MgCl 2 , 0.2 μM of each primer, 1 unit of Taq DNA polymerase (Fermentas International Inc., Canada) and~15-20 ng of DNA. An initial denaturing step of 3 minutes was performed at 96°C, followed by 35 cycles of 1 minute at 95°C, 1 minute at 56°C and 2 minutes at 72°C, and a final elongation step of 10 minutes at 72°C. The PCR products were purified and sequenced as before. Sequences were trimmed, evaluated for quality and assembled in contigs using the Geneious 6.0.6 software (Biomatters Ltd.; Auckland, NZ). When heterozygotes with alleles containing insertions/deletions (indels) of different lengths within intron iII were detected, PCR amplifications were repeated and the products were purified and cloned into the pGEM-T vector (Promega Corp). At least 6 clones from each individual were sequenced, in order to obtain sequences from both alleles.

Sequences analyses
Sequences from all the 201 individuals were aligned with the algorithms of ClustalX2 [46] and Muscle [47], manually edited in Mega 5.2 [48] and inspected for quality with the Gblock Server 0.91b [49]. Mitochondrial sequences were aligned considering each fragment separately (CR, COI and COII) and subsequently concatenated for analyses. COI, COII and ChαE7 exons (eIII and eIV) were checked for open reading frames using Mega 5.2 [48]. For the CR fragment and ChαE7 intron iII sequences, each indel was considered as a single mutation and all indels were recoded as single positions in the final alignment. Individual sequences were collapsed into unambiguous haplotypes using FaBox [50], and these were compared and named according to sequences reported previously in the case of the mitochondrial data [38][39][40].

Population analyses
General diversity indices, such as number of haplotypes, haplotype diversity or expected heterozygosity (Ĥ) and nucleotide diversity (π) [51], were calculated for each dataset (mitochondrial CR, COI, COII and nuclear ChαE7) independently using the DnaSP 5.0 [52] and Arlequin 3.5 [53] softwares. The genetic structure of NWS fly populations was first investigated using the Analysis of Molecular Variance (AMOVA) and pairwise F ST . This was performed using the Arlequin 3.5 software [53] separately for each dataset and considered each sampled location as a population and all as belonging to a unique population. The statistical significance values were accessed using 10,000 permutations. We also carried out a discriminant analysis of principal components (DAPC) [54] implemented in adegenet 1.3-8 from the R package [55]. This is a multivariate analysis without an a priori population model that identifies clusters (K) of individuals independent from geographic origin. The analysis was performed considering K = 2 to K = 10.
In order to describe genetic and geographically concordant groups of sampling locations we used spatial analysis of molecular variance (SAMOVA) [56]. This was performed based on genetic pairwise differences between sequences and geographic coordinates of sampling locations, through 100 annealing processes and 1,000 permutations using the SAMOVA software [56]. Different numbers of groups were tested (K = 2 to K = 10) and the F SC , F ST and F CT indices were plotted to investigate the number of groups (K) that best represent the spatial distribution of genetic variability.
Mantel tests [57] were performed using the mitochondrial and nuclear datasets to detect signals of isolation by distance (IBD). These tests were carried out considering a genetic distance matrix, constructed based on Slatkin's linearized F ST 's [58], and a geographic distance matrix, that contains linear distances in kilometers between each pair of locations. Tests were performed using the Isolation By Distance Web Service (IBDWS) 3.23 [59] through 10,000 randomizations.

Demographic analyses
Tajima's D [60], Fu's Fs [61], Fu and Li's D Ã and F Ã [62] summary statistics were estimated after defining the best grouping configuration of sampling locations to investigate deviations from neutrality of NWS fly samples. For all tests, values near zero indicate that population sizes are stable, significant negative values indicate an excess of low frequency variants which are associated with positive selection or a demographic scenario of population expansion, and positive values indicate that a population has an excess of intermediary frequency variants resulting from the action of a balancing selection or population bottleneck [60][61][62]. Tests were performed for the concatenated mitochondrial sequences (CR, COI, COII) and each part of the ChαE7 gene separately (introns iII and iIII, exons eIII and eIV) using DnaSP 5.0 [52]. The significance of Fu's Fs was obtained by coalescence, with 1,000 replicates and a 95% confidence interval, and two tailed tests were used to determine the statistical significances of the other three tests.
Mismatch distribution analysis based on the distribution of pairwise differences between sequences was performed to better understand the demographic scenario of NWS fly samples. A unimodal distribution indicates recent expansion, while a multimodal distribution indicates population contraction, structure, influence of migration or subdivision [63][64][65][66][67][68][69]. This analysis was performed using the Arlequin 3.5 software [53], with the significance obtained by 10,000 permutations. The sum of the square deviations (SSD) and the Raggedness index (rg) were used to measure the adjustment to a population expansion model.

Positive selection
The investigation of ChαE7 positive selection signals was done based on the population structure of the sampling locations.
The frequencies of polymorphic nucleotides in the two codons associated to insecticide resistance (Gly137Asp and Trp251Leu) and the Hardy-Weinberg equilibrium between these were obtained with the Arlequin 3.5 software [53].
Fay and Wu's H neutrality test [70] was performed using the DnaSP 5.0 software [52]. The orthologous sequence of D. melanogaster (GenBank accession number: NM079537) was used as the outgroup and statistical significance was accessed by coalescence through 1000 replicates, considering a confidence interval of 95%.

Carboxylesterase E3 (ChαE7) Gene Characterization
ChαE7 sequences from five individuals were obtained (GenBank accession numbers: KR067547-KR067551), although two of these were partial due to the absence of the initial sequence of intron iII. Sequences were approximately 3 Kb, comprising 111 bp of exon eII, around 2000 bp of intron iII (variable sizes due to presence of indels), 150 bp of exon eIII, 63 bp of intron iIII and 318 bp of exon eIV (Fig 1). Although the NWS fly eII, eIII, iIII and eIV regions have been described previously [41], this is the first report of the intron iII sequence. The two mutations associated with insecticide resistance, Gly137Asp and Trp251Leu, are located in exons eIII and eIV, respectively.
Nucleotide composition of exons (eIII and eIV) and introns (iII and iIII) of the five NWS fly individuals and the other Muscomorpha species (L. cuprina, D. melanogaster and M. domestica) is presented in Table 2. For introns, M. domestica was excluded due to absence of sequence.
A high A+T content is observed in both exons and introns of NWS fly individuals, with minimum and maximum values of 60.06% and 80.95%. This was related to a higher proportion of T than A, indicated by the negative AT-skew. The GC-skew shows positive values for all individuals when exons were considered (indicating more Gs than Cs), while this was negative for most of the individuals in the case of the introns. However, both AT-skew and GC-skew were near to zero for intron iII (0.01-0.05), indicating that A and T nucleotides, as well as G and C, occurred in similar proportions.
L. cuprina had a more similar pattern to the NWS fly, consistent with their closer phylogenetic relationship, with an A+T content of 64.67% and 56.29% for eIII and eIV exons, respectively. The AT-skew was negative and the GC-skew positive for both exons of L. cuprina. However, the values were near to zero for exon eIV (-0.04 and 0.01).
D. melanogaster and M. domestica showed a more balanced distribution of nucleotide frequencies. However, D. melanogaster also had a high A+T content for intron iIII and iII sequences with values of 64.91 and 76.19%, respectively.
Haplotype distribution in sampling locations, nucleotide diversity (π) and haplotype diversity/expected heterozygosity (Ĥ) estimated from concatenated mitochondrial and ChαE7 gene sequences are shown in Table 3. Low nucleotide diversity was observed for both mitochondrial and nuclear genes, ranging from 0.0015 to 0.0055 in mitochondrial and 0.0052 to 0.0227 in the ChαE7 datasets. However, high haplotype diversity was found, ranging from 0.7556 to 1 and 0.1895 to 0.8929 in the mitochondrial and ChαE7 datasets, respectively.  (2)

Structure analyses
The non-hierarchical AMOVA, considering all sampled locations as a unique population, showed higher structure when considering the ChαE7 gene, with fixation indexes (F ST ) of 0.235 and 0.094 for the ChαE7 and mitochondrial datasets respectively (Table 4). Pairwise F ST estimates with the ChαE7 gene showed significant values between approximately 65% of all locations pairs (137/210) (S1 Table), decreasing to 33% (70/210) when considering mitochondrial data (S2 Table). In the latter case, BCR, BSA and APL were the locations which had the most significant values (18/20, 14/20 and 12/20, respectively).
Examination of the mitochondrial data showed that the NWS fly populations were genetically and geographically homogeneous, according to DAPC and SAMOVA analyses. In addition, application of the Mantel test showed no correlation between geographic and genetic distances (r = 0.0445, p > 0.05; S1 Fig), indicating an absence of isolation-by-distance and no limited gene flow among NWS fly sampled locations. Therefore, all locations analyzed were considered as belonging to a unique group (K = 1) based on the mitochondrial data.
In contrast, an increased degree of genetic structure was detectable based on the carboxylesterase E3 (ChαE7) gene according to DAPC and SAMOVA analyses, with K = 3 as the most probable geographic arrangement (Fig 2). SAMOVA indicated that K = 3 was the most probable scenario due to high F CT and low F SC values (Fig 2A). A North-South pattern was observed based on the DAPC results, with a higher frequency of individuals belonging to group Red and group Green in higher and lower latitudes, respectively (Fig 2B). Group Yellow was more frequent in some southern locations, such as BSA and UCO. The three genetically and geographically distinct groups, based on both analyses, were: I) BTO; II) BGN, BGO, BCA, BSS, BES, BFV and III) BCR, BAQ, BCI, PYB, BSA, UST, UPM, UDA, BPM, UBM, UCC, UCO, UJS, APL (Fig 2B). Mantel testing indicated no correlation between geographic and genetic distances (r = 0.2458, p = 0.05) (S1 Fig), with absence of isolation-by-distance among NWS flysampled locations.

Demographic inference
The demographic history of NWS fly was inferred based on the groups defined by DAPC and SAMOVA, considering one group for mitochondrial data and three for the ChαE7 gene. Tajima's D, Fu's Fs and Fu & Li's D Ã and F Ã summary statistics yielded significant negative values for the mitochondrial data, while no significant results were found for the ChαE7 gene, except the Tajima's D and Fu & Li's F Ã for the intron iII in group I (Table 5). Mismatch distribution analysis resulted in a unimodal distribution of pairwise differences between mitochondrial sequences, adjusted to the theoretical distribution of a demographic scenario of recent population expansion (P [SSD] = 0.4442, P [raggedness] = 0.6298). This analysis resulted in multimodal distributions for the three groups of populations with significant values of SSD and raggedness index for the ChαE7 gene, rejecting the null hypothesis of population expansion (Fig 3).

Positive selection
From a total of 49 alleles, 31 had a glycine and 18 had an aspartate in position 137, while 43 had a tryptophan, 5 a serine and 1 a leucine in position 251 (S3 Table). The presence of glycine (coded by GGG or GGC) and aspartate (coded by GAC or GAT) in position 137 are associated, respectively, with susceptibility and resistance to diethyl-organophosphates in dipterous species. The presence of tryptophan in position 251 (coded by TGG) defines the susceptible genotype, while serine and leucine (coded by TCG and TTG, respectively) are associated with dimethyl-organophosphates resistance. The substitution of the amino acid tryptophan by a leucine in position 251 was associated previously with dimethyl-organophosphates resistance in L. cuprina [24]. However, another study found that dimethyl-organophosphate resistant lineages of M. domestica have a serine substitution in place of a tryptophan [71], suggesting that substitution with small residues, such as leucine, glycine or serine could be the factor that leads to the increase of insecticide hydrolysis by carboxylesterase E3. No alleles with residues predicted to confer resistance in both sites positions were observed (i.e. aspartate in position 137; serine or leucine in position 251).
Residue combinations in these two positions for 25 susceptible alleles and 24 resistant alleles did not reflect the real frequency of susceptibility and resistance in field populations. To explore this issue further, we specifically analyzed the two residues associated with insecticide resistance at the nucleotide level and considering allele frequencies taking into account the three groups for the ChαE7 gene.
Considering codon 137, the nucleotide present in the second position defines the susceptibility (guanine; named Gly137-G) or resistance (adenine; named Asp137-A). Similarly, the nucleotide present in the second position of codon 251 defines susceptibility (guanine; named Trp251-G) and the two resistant types (thymine and cytosine; named Leu251-T and Ser251-C, respectively). The frequencies of each nucleotide in these two codon positions for each group of sampling locations are showed in Fig 4. Group I (consisting only of BTO) showed the same frequencies of Gly137-G and Ser251-C (56%), and Asp137-A and Trp251-G (44%). Group II showed the same pattern of complementary frequencies for the Gly137-G and Ser251-C (75%) and for Asp137-A and Trp251-G (25%). In group III, the frequency of Asp137-A was 64% and that of Trp251-G was 82%, while the frequencies of Gly137-G and Ser251-C were 36% and 16%, respectively. This showed that there was a higher frequency of individuals in groups I and II with a genotype conferring resistance to dimethyl-organophosphates, while there was a higher frequency of diethyl-organophosphates resistant flies in group III. Another difference between these groups is the low frequency (less than 3%) of Leu251-T that appeared only in group III.
These two polymorphic codon positions (second codon positions from codons 137 and 251) were not in Hardy-Weinberg equilibrium in the three groups of sampled locations, as the observed heterozygosity was smaller than expected (Table 6).  Discussion Partial characterization of the carboxylesterase E3 ChαE7 gene allowed its use as a nuclear molecular marker for comparisons with previously studied mitochondrial markers [38,39] and for investigation of the genetic and geographic structure, demographic history and selection pressure of NWS fly populations inhabiting the southern region of the Amazon river basin. Based on the homologous ChαE7 sequence from L. cuprina [27], we expected to generate a fragment of around 921 bp following PCR amplification, instead of the approximately 2600 bp obtained for the NWS fly. This difference was due to the increased length of intron iII in the NWS fly. L. cuprina has a length of only approximately 280 bp [27], which is almost seven times smaller than the observed length for NWS fly (around 2000 bp). As a further comparison, the length of intron iII in D. melanogaster is 85 bp [72] and that in M. domestica is 63 bp [26]. This indicates that the length of intron iII has high variability within Muscomorpha. In contrast, intron iIII showed a smaller variation in length in the same species, as this was 63 bp in the NWS fly, 56 bp in L. cuprina [27], 58 bp in D. melanogaster [72] and 62 bp in M. domestica [26]. Unlike introns, the sequences of exons eIII and eIV are conserved with respect to the length, with only the D. melanogaster exon eIV presenting three more base pairs in the region analyzed, totalizing 321bp.
The NWS fly ChαE7 nucleotide composition was similar to that of L. cuprina, consistent with their closer phylogenetic relationship. Both species showed a bias toward A+T content, which was suggested previously to be a remarkable feature of the mitochondrial genome of Calliphoridae species [73,74]. In general, a negative AT-skew and a positive GC-skew was observed, indicating different nucleotides usage. The skewness is a measure of the relative number of As to Ts, in the case of AT-skew, and Gs to Cs, in the case of GC-skew [75].
The low genetic distance among haplotypes in both datasets reflects the low nucleotide diversity (π), although haplotype diversity (Ĥ) was high in both cases. This pattern of diversity is typically observed in a demographic scenario of recent populations expansion [76].
All population analyses pointed to an increased genetic structure of the carboxylesterase E3 ChαE7 gene in relation to the mitochondrial gene sequences. This was indicated by the AMOVA F ST value which was 2.5 times greater and a higher percentage of pairwise F ST significant values for the ChαE7 gene. In a similar manner, the DAPC and SAMOVA data showed the ChαE7 genetic structure of NWS fly populations, revealing three geographically non-concordant groups (K = 3). There was also an apparent North-South pattern of variability distribution, which did not appear to result from isolation-by-distance (IBD), as indicated by the non-significant Mantel tests between genetic and geographic distances. A previous investigation suggested that IBD can lead to false positives in tests of population structure and detection of loci evolving under selection [77], although the present finding of no correlation between geographic and genetic distances matrices indicates that the population structure detected is not biased or confounded by IBD. In this scenario, the three ChαE7 groups could result from selection because of regional differences of insecticide usage. For mitochondrial data, the AMOVA F ST value was low and the BCR, BSA and APL locations showed the highest percentage of pairwise F ST significant values. BCR is located in Mato Grosso do Sul (MS), the central region of Brazil, while BSA is located in Rio Grande do Sul (RS), the southern region of Brazil. Both states are noted for high Brazilian cattle production, being the second and tenth in amount of slaughter in 2013, respectively [78]. Not all locations sampled in these two states showed an increased percentage of significant pairwise F ST values, indicating an unusual pattern of both locations which could be linked to the central points for cattle interchange between places. APL is the unique population sampled in Argentina (located near southern Brazil), which is also a region of high cattle production. The results of the DAPC and SAMOVA analyses supported a panmixia scenario for the mitochondrial gene data, consistent with a previous study [38], and therefore did not reveal a genetic/geographic structure of the sampled locations.
These findings are similar to results from other insect pests. The codling moth Cydia pomonella exhibits a low level of genetic structure in four countries (France, Italy, Armenia and Chile) based on seven microsatellite loci, probably due to its high capacity of migration, while a higher population differentiation was detected with insecticide resistance markers [79]. In addition, population differences in resistance are larger than differences in microsatellite variation for the Australian diamondback moth (Plutella xylostella), without evidence of spatial structure for resistance through IBD [80].
Tajima's D, Fu's Fs and Fu & Li's D Ã and F Ã significant negative values and unimodal distribution of pairwise differences indicate a history of population expansion according to mitochondrial gene data. However, the finding of no significant results for these neutrality tests and multimodal distribution of pairwise differences for the ChαE7 gene indicates a different history which, in addition to the population analysis, supports the potential role of selection.
The results of Fay & Wu's H [70] neutrality test indicated positive selection on both exons (eIII and eIV) from groups II and III but not for group I. The absence of a positive selection signal in the case of group I could be due to sampling size since it was composed of only nine individuals from one location (BTO). Resistant genotypes had a high frequency in all sampled areas (Fig 4), with dimethyl-organophosphate resistance represented more in lower latitudes and diethyl-organophosphate resistance was greater in higher latitudes. These differences of ChαE7 allele frequencies could reflect the usage of different insecticide compounds in each geographic location. Alternatively, they could reflect differences in alleles frequencies of a Modifier locus, which has been related to a fitness compensation of the negative effects of both resistant mutations Gly137Asp and Trp251Leu (i.e. developmental instability) in L. cuprina [81][82][83].
High levels of insecticide resistance are also observed for Aedes aegypti in three Brazilian states (Rio de Janeiro, Espírito Santo and Ceará) [84,85] and for the house fly Musca domestica in three Argentinean populations [86]. The horn fly, Haematobia irritans, showed high frequencies of the kdr mutation, associated with high levels of pyrethroid resistance, in southern Chile [87]. Our results, in addition to many other studies, indicate that insecticide resistance is widespread throughout South America. It is likely that use of different compounds for control of distinct insect species has led to resistance through multiple mechanisms, which should be cause for great concern.
A negative association between the mutations in the carboxylesterase E3 ChαE7 gene that confer resistance was detected, indicating that they are not present in the same allele. This was reported previously for L. cuprina [88] and NWS fly in four of the locations under investigation in the present study (UST, UDA, UBM and UCC) [89].
Taken together, the observed low heterozygosity in the two polymorphic codon positions associated with susceptibility and resistance, the negative association between nucleotides in the same codons and the significant negative Fay & Wu's H values for exons suggest that the ChαE7 gene may be undergoing evolution under positive selection. The present findings indicate that this could be due to the use of organophosphate compounds on NWS fly populations in South America. This is consistent with previous studies on the demographic history of this species [38,39]. Many recent population studies have investigated effects on mitochondrial and nuclear markers together, in attempts to understand the history of the species through genetic variation comparisons (for example, see [90,91]).
Livestock production began in the Americas around 500 years ago, resulting in the likely passive distribution and unclear geographic structure of this species and of some concordant insect pests. For example, approximately half of tick species have a genetic structure based mainly on the movement capacity of their respective hosts. Less mobile hosts typically support a low gene flow with highly structured tick populations. In contrast, highly mobile hosts allow a higher gene flow between ticks, leading to a low population structure [92]. As in the case of tick populations, the NWS fly has a notorious relation with its host during the larval stage and can be dispersed through host movement. However, the geographic structure of NWS fly populations has not been completely resolved and no clear boundaries between groups of populations have been identified. Genotyping-By-Sequencing (GBS) is a recent promising method that could be used to help achieve this aim, as this could lead to identification of hundreds or thousands of molecular markers (single nucleotide polymorphisms, SNPs) without the requirement of a completely sequenced genome [93].
When selecting a target region for a SIT based program, it is important to analyze the insecticide susceptibility of populations. Here, we have confirmed that the two mutations in the carboxylesterase E3 ChαE7 gene associated with organophosphate resistance in many species, Gly137Asp and Trp251Leu, are also related to this mechanism in NWS fly through Darwinian selection. Additionally, we showed that NWS fly populations from SAG have high frequencies of organophosphate resistant genotypes, which can represent a barrier for the success of a control program based on SIT.

Conclusions
The genetic data presented here showed that at least three groups of NWS fly populations can be found in SAG, but with no clear geographic correlation, and resistance to organophosphate insecticides is widespread and found in high frequency in most of the analyzed locations. These aspects may represent barriers to implement a successful control program based on SIT. New studies based on dense SNP genotyping will be necessary in order to refine previous results of population structure and to identify additional molecular markers that could have evolved under positive selection. This will allow the detection and monitoring of other genes involved in insecticide resistance and their frequencies in natural populations.