Next Article in Journal
Evaluation of the Insecticidal Activities of α-Pinene and 3-Carene on Sitophilus zeamais Motschulsky (Coleoptera: Curculionidae)
Next Article in Special Issue
Population Genetic Structure and Demography of the Critically Endangered Chequered Blue Butterfly (Scolitantides orion) in a Highly Isolated Part of Its Distribution Range
Previous Article in Journal
Synergism between Hydramethylnon and Metarhizium anisopliae and Their Influence on the Gut Microbiome of Blattella germanica (L.)
Previous Article in Special Issue
DNA Barcoding: A Reliable Method for the Identification of Thrips Species (Thysanoptera, Thripidae) Collected on Sticky Traps in Onion Fields
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Genetic Differentiation of a New World Screwworm Fly Population from Uruguay Detected by SNPs, Mitochondrial DNA and Microsatellites in Two Consecutive Years

by
Luana Walravens Bergamo
1,2,3,*,
Karina Lucas Silva-Brandão
3,4,
Renato Vicentini
1,
Pablo Fresia
5,* and
Ana Maria Lima Azeredo-Espin
1,3,*
1
Departamento de Genética, Evolução, Microbiologia e Imunologia, Instituto de Biologia, Universidade Estadual de Campinas (UNICAMP), Campinas SP 13083-970, Brazil
2
Programa de Pós-Graduação em Genética e Biologia Molecular, Universidade Estadual de Campinas (UNICAMP), Campinas SP 13083-862, Brazil
3
Centro de Biologia Molecular e Engenharia Genética, Universidade Estadual de Campinas (CBMEG-UNICAMP), Campinas SP 13083-875, Brazil
4
Centro de Ciências Naturais e Humanas, Universidade Federal do ABC (CCNH-UFABC), Santo André SP 09210-580, Brazil
5
Unidad Mixta Pasteur + INIA (UMPI), Institut Pasteur de Montevideo, Montevideo 11400, Uruguay
*
Authors to whom correspondence should be addressed.
Insects 2020, 11(8), 539; https://doi.org/10.3390/insects11080539
Submission received: 9 July 2020 / Revised: 2 August 2020 / Accepted: 9 August 2020 / Published: 16 August 2020
(This article belongs to the Special Issue Population Genetics of Insects)

Abstract

:

Simple Summary

The New World screwworm (NWS) fly is a pest related to economic impacts in livestock breeding, mainly because females lay their eggs in wounds and natural cavities of warm-blooded vertebrates. A direct way to control this pest is the sterile insect technique (SIT), which is based on the release of sexually sterilized males to mate with wild females, resulting in the absence of new flies. The first step to the implementation of such program is to know how flies’ populations are structured, that is, if they can freely mate in the field or, otherwise, if they are impeded of mate due to spatial or temporal barriers. Here we investigated genomic polymorphisms to infer temporal structure in a population of the NWS fly from Uruguay in two consecutive years (approximately 18 generations). We found that this time was enough for this population to change its genetic composition, but it was not sufficient to prevent free mate among individuals from one year to the next. We also found signal that some of the polymorphisms are related to the response of the population to insecticide chemical control. This approach can be used in the future to estimate spatial barriers to geographically isolated populations.

Abstract

The New World screwworm (NWS) fly, Cochliomyia hominivorax (Diptera: Calliphoridae), is an economically important ectoparasite currently distributed in South America and in the Caribbean basin. The successful eradication of this species in USA, Mexico and continental Central America was achieved by a control program based on the sterile insect technique (SIT). In order to implement a genetic control strategy over the NWS fly’s current area of occurrence, first, it is necessary to understand the species dynamics and population structure. In order to address this objective, the spatial genetic structure of the NWS fly was previously reported in South America based on different genetic markers; however, to date, no study has investigated temporal changes in the genetic composition of its populations. In the current study, the temporal genetic structure of a NWS fly population from Uruguay was investigated through two consecutive samplings from the same locality over an interval of approximately 18 generations. The genetic structure was accessed with neutral and under selection SNPs obtained with genotyping-by-sequencing. The results gathered with these data were compared to estimates achieved with mitochondrial DNA sequences and eight microsatellite markers. Temporal changes in the genetic composition were revealed by all three molecular markers, which may be attributed to seasonal changes in the NWS fly’s southern distribution. SNPs were employed for the first time for estimating the genetic structure in a NWS fly population; these results provide new clues and perspectives on its population genetic structure. This approach could have significant implications for the planning and implementation of management programs.

1. Introduction

Population genomics is an area of considerable growth, accompanying the development of next-generation sequencing (NGS) technologies, with the potential to improve population studies [1]. Many distinct approaches to obtain hundreds to thousands of molecular markers throughout the genome were developed in recent times [2,3]. One of them is genotyping-by-sequencing (GBS) [4,5], which is based on restriction enzymes to reduce genome representation and is capable of obtaining a huge number of single-nucleotide polymorphisms (SNPs) for several individuals at the same time. In the last years, several studies have been conducted using the GBS technique or similar methodologies for non-model insects [6,7,8,9,10,11,12].
This approach, which combines the sampling of many individuals and the characterization of a huge number of markers, allows the development of temporal studies, where samples from the same population are taken at different times with the aim to observe temporal changes in allelic frequencies, offering insights about the species population dynamics and demography. Effective population size (Ne), an important demographic parameter usually difficult to accurately estimate, can be more robustly assessed with temporal sampling and population genomic data [13,14]. Most studies that aimed to infer demographic dynamics based on genetic data approached that based on historical timescales (i.e., thousands of generations) and not on an ecological timescale (i.e., tens of generations) [15]. In fact, the role of genetic variation on ecological differences must be considered when leading with both temporal and spatial allopatric populations [16]. Understanding demographic changes at an ecological timescale is particularly important for the New World screwworm (NWS) fly populations in its southern distribution (e.g., Uruguay), where the species is probably impacted by seasonal changes throughout the year.
The NWS fly, Cochliomyia hominivorax (Coquerel, 1858) (Diptera: Calliphoridae), is an ectoparasite that causes myiasis in warm-blooded vertebrates, including species of economic importance, consequently leading to negative economic impacts in livestock breeding [17]. In this species, females lay their eggs in wounds and natural host cavities. The emerged larvae feed on host living tissues and after three instars they fall to the ground to pupate, giving rise to adults after some days [17]. The species is intolerant to low temperatures and screwworms eggs fail to develop at temperatures lower than 12.3 °C [18]. The species is currently distributed over the Neotropical region, but its distribution originally extended from the Southern USA to Argentina [19]. This territory reduction involved systematic area-wide management programs based on the sterile insect technique (SIT) that eradicated the species from North America and continental Central America [20].
Several population genetic studies have been conducted in recent years in order to better understand the dynamics and distribution of the NWS fly (reviewed in [21]) with the aim to give basis for future implementation of similar control programs in its current area of occurrence. In South America, two genetic groups were identified based on different molecular markers, the North Amazon group (NAG) and the South Amazon group (SAG) [22,23]. Within the SAG, the results of some of the studies diverged about the finer scale in which genetic differentiation can be detected [24,25,26,27], which can be due to the evolutionary rate of the molecular markers investigated (mitochondrial DNA, microsatellites or LcαE7 sequences). For this reason, the use of more informative molecular markers, such as SNPs obtained with the GBS methodology, could provide higher resolution in population genetic studies, allowing establishing boundaries (i.e., neighborhood size) and management units.
In the current study, we standardized the two-enzyme GBS protocol for the NWS fly and used SNPs markers to investigate temporal genetic changes in a Uruguayan NWS fly population at an ecological timescale. For a better understanding of the information provided by these markers we compared the obtained results with estimates based on mtDNA sequences and microsatellites for the same samples. Our study is the first to implement the GBS technique for this pest species, allowing its future implementation in spatial NWS fly population studies.

2. Material and Methods

2.1. Sampling Information

NWS fly samples were collected from wounds and natural cavities of sheep and cattle from the locality of Cañas (32°21′21.96″ S, 53°49′45.1″ W—Cerro Largo, Uruguay) in two distinct periods: February/March 2015 and March 2016, hereafter named as 2015-Sample and 2016-Sample, respectively. Genetic property was registered under SISGEN #A61EDD2. Larvae were collected from only one wound per animal to lower bias and were kept in absolute ethanol and freezer at −20 °C.

2.2. DNA Extraction

GBS libraries demand high quality DNA and the best quality total genomic DNA was obtained by extracting one to three larvae per wound using the CTAB method [28], with an additional step of RNAse treatment. Total DNA was eluted in 50 µL of AE buffer and stored at −20 °C. Extractions were quantified by fluorescence using the Qubit DNA quantification system (Invitrogen, Carlsbad, CA, USA), and their quality and purity were verified using a NanoDrop UV spectrophotometer (Techno Scientific, Wilmington, DE, USA). High quality DNA (260:280 ratio > 1.8, 230:260 ratio > 2.0) was obtained for 32 and 31 samples from 2015-Sample and 2016-Sample, respectively (a 45% extraction success). The DNA concentration per sample was normalized to 20 ng/µL for GBS libraries construction.

2.3. Single-Nucleotide Polymorphisms

SNPs were obtained using the genotyping-by-sequencing (GBS) approach. Two separated libraries were constructed with the individuals from 2015-Sample and 2016-Sample, based on standard protocols [4,5] with modifications described below. Total DNA was double-digested with two high fidelity restriction enzymes, PstI–MspI (New England Biolabs—NEB, Ipswich, MA, USA). Barcoded adapters were ligated on the rare cut site (i.e., PstI cut site) for each individual separately, and the common adapter was subsequently ligated on the common cut site (i.e., MspI cut site) for all individuals, using T4 DNA ligase (NEB). The samples were separated by year and pooled in two sets, which were amplified by multiplexed PCRs using standard forward primer A and modified reverse primer C with the Q5®® high-fidelity DNA polymerase (NEB). PCR products were purified with the Agencourt Ampure XP beads (Beckman-Coulter, Inc., Brea, CA, USA), and the DNA amount was estimated with a NanoDrop UV spectrophotometer (Techno Scientific, Wilmington, DE, USA). Library quality was accessed with the evaluation of fragment sizes and the presence of adapter dimers in an Agilent BioAnalyzer 2100 at Life Sciences Core Facility (LaCTAD) from University of Campinas (UNICAMP). The two libraries were sequenced separately in an Illumina NextSeq 500 lane (Illumina, Inc. San Diego, CA, USA), using 150-bp single-end reads, at the IBTEC—UNESP (Botucatu, Brazil).
Stacks v1.48 de novo pipeline was used to demultiplex and perform SNP calling without a reference genome [29,30]. First, the process_radtags was used for raw data demultiplexing, removing low quality reads and with uncalled bases reads. Adaptors and restriction enzyme recognition sites were trimmed, and reads were filtered to consider only the initial 70 bp. We excluded 11 individuals after this process (five from 2015-Sample and six from 2016-Sample), due to the low number of reads (less than 1% from the total number of reads in the sequenced lane). This cleaned dataset was the input of ustacks function to assemble de novo loci present in each individual and to SNPs calling (applied parameters include minimum stack depth of 3 and allowed distance between stacks of 2 mismatches). Then, the data from each individual from both samples were merged into a catalog of loci using the cstacks function (allowed distance between catalog loci = 2), which contains all loci and alleles in the considered populations. Subsequently, each individual was matched against the catalog with sstacks. The function populations was executed for additional filtering considering (i) locus that must be present in both sampled years, (ii) only the first SNP of each locus, (iii) locus in at least 75% of individuals and (iv) a minimum minor allele frequency of 0.10.
Lastly, additional data filtering included pruning loci under Hardy–Weinberg disequilibrium and loci under selection from the two samples. Hardy–Weinberg tests were conducted in the web version of Genepop [31,32] considering a p-value of 0.05. Loci putatively under selection were detected with Lositan [33], which calculates the ratio between FST and expected heterozygosity (He) values from each locus and compares them with the neutral distribution to detect outliers (i.e., loci under putative selective pressure) [34]. The program was run three times to lower bias [33]: (i) the first run considered all loci to calculate a neutral mean FST, using 50,000 simulations, 99% confidence interval, infinite alleles mutation model and false discovery rate of 0.1%; (ii) the second run was conducted with the same parameters, except that only neutral loci from the first run were considered to recalculate the mean neutral FST; (iii) the third run considered all loci again and the neutral FST calculated in the second run. Outlier loci in this last run were inferred as under selection and consequently excluded from the dataset used in population genetic analyses. These outlier loci sequences detected by Lositan were blasted, mapped and annotated against the National Center for Biotechnology Information Search database (NCBI) using Blast2Go [35] for gene ontology annotation. File conversions were conducted using PGDSpider v2.0.5.1 [36] for further population genetic analyses.

2.4. mtDNA

Three mitochondrial DNA fragments were amplified according to the procedures described in [23]: the B domain of the control region (CR) (ca. 541 bp) [37], the cytochrome c oxidase subunit I gene (COI) (ca. 892 bp) and the cytochrome c oxidase subunit II gene (COII) (ca. 686 bp). Purified PCR products were sequenced using the BigDye terminator cycle sequencing kit v3.0 (Applied Biosystems, Foster City, CA, USA) in an ABI 3730 automated DNA sequencer (Applied Biosystems, Foster City, CA, USA) at the Life Sciences Core Facility (LaCTAD) from State University of Campinas (UNICAMP).
CR, COI and COII sequences were independently aligned using ClustalX [38] and manually inspected in Mega v6 [39]. COI and COII fragments were translated to protein to confirm the open reading frame. Each insertion/deletion (indel) in the CR fragment was considered a single mutational step and recoded as a single position in the final alignment in Mega. Individual sequences of the three fragments were collapsed into haplotypes with Fabox (available at http://users-birc.au.dk/palle/php/fabox/index.php). Sequences of the three fragments for each individual were concatenated in a unique haplotype for further analyses.

2.5. Microsatellites

Eight polymorphic microsatellite loci were amplified and genotyped according to the procedures previously described [40,41]. PCR products were verified with 6% denaturing polyacrylamide gels stained with silver nitrate and the positive amplicons were multiplexed (4-plex, using different fluorophores (FAM, NED, PET and VIC)) for subsequent sequencing in an ABI 3500 genetic analyzer (Applied Biosystems). Each reaction received 0.4 μL of the molecular weight marker (Liz GeneScan 600 v2, Thermo Fisher Scientific), 9 μL of HI–DI formamide (Thermo Fisher Scientific) and up to 2 μL of the PCR products diluted in multiplex.
Microsatellite alleles were genotyped by size using GeneMarker v2.4.2 [42]. The presence and frequency of null alleles were estimated with Microchecker [43], which was also used to adjust allelic and genotypic frequencies.

2.6. Genetic Diversity and Population Differentiation

Genetic diversity and population structure analyses were conducted for each molecular marker in order to detect temporal changes in the NWS fly population from Cañas.
Filtered neutral SNPs were used to investigate nucleotide diversity. Observed heterozygosity (Ho), expected heterozygosity (He) and Wright’s F-statistic (FIS) were estimated with the function populations available in the Stacks pipeline. Genetic differentiation between 2015-Sample and 2016-Sample was estimated with a nonhierarchical AMOVA with Arlequin v3.5 [44], with the Bayesian approach implemented in Structure v2.3.3 [45] and with a discriminant analysis of principal components (DAPC) [46]. AMOVA significance was obtained by 10,000 permutations. The parameters for the Structure analysis included admixture model with correlated allele frequencies, 100,000 MCMC replications, 10,000 burn-in, K varying from 1 to 5 and 10 iterations for each K value. Evanno et al. method [47] was used to determine the best K value using the application Structure Harvester v0.6.94 [48]. Results from this best K value were summarized with Clumpp v1.1.2 [49] and represented with Distruct v1.1 [50]. The DAPC analysis was conducted with the R package adegenet [51]. The function find.clusters was used to determine genetic clusters, based on principal component analysis (PCA) retaining the maximum number of PCs and the Bayesian information criterion (BIC) to generate a graph of optimum K. This procedure produced a scatterplot and a smaller number of PCs were retained then, following recommendations in the package. Subsequently, the contributions of the alleles to the best clustering scenario was also computed, considering a threshold of 0.005.
Mitochondrial DNA haplotype diversity (Ĥ), nucleotide diversity (π) and overall genetic structure (nonhierarchical AMOVA) were estimated using Arlequin v3.5 [44] with statistical significance obtained with 10,000 permutations. A parsimony haplotype network was constructed with PopART v1.7 [52] using the TCS method with 95% connection limit [53] to recover relationships among concatenated sequences.
The genetic variability of microsatellite loci was accessed through allelic richness (AR) with MSA v4.05 [54]; observed (Ho) and expected (He) heterozygosities with GenAlex v6.5 [55]; and estimates of FIS with the web version of Genepop [31,32]. Each locus was tested for Hardy–Weinberg equilibrium with the web version of Genepop [31,32], considering the Markov chain model from [56] and using the following parameters: 1000 series, 1000 iterations and 1000 steps of “dememorization”. Linkage disequilibrium between pairs of loci within each sample was investigated using the web version of Genepop [31,32]. Genetic structure was estimated by a nonhierarchical AMOVA with Arlequin v3.5 [44], using 10,000 permutations to access statistical significance and by a Bayesian clustering analyses with Structure 2.3.3 [45], using the same parameters described above.

2.7. Demographic Inferences

For both SNPs and microsatellite data, two-sample effective population size (Ne) estimates were calculated using Jorde and Ryman method [57] in NeEstimator v2.01 [58]. This analysis considered a minimum allele frequency (MAF) cutoff of 0.01, 0.02 and 0.05 and 18 generations/year (considering that the NWS fly generation time is approximately 21 days under laboratory conditions at 25 °C [17]).
The demographic history of the NWS fly samples was inferred from mitochondrial data with Tajima’s D [59] and Fu’s Fs [60] neutrality tests and mismatch distribution analysis using Arlequin v3.5 [44]. For both summary statistics, values near zero are indicative of population size stability, negative values are indicative of recent population expansion and positive values are indicative of population bottlenecks [59,60]. Mismatch distribution analysis was performed with 10,000 permutations for statistical significance. The Raggedness index (rg) was used to measure the adjustment to a population expansion model.

3. Results

3.1. Single-Nucleotide Polymorphisms

Sequencing of both 2015 and 2016 GBS libraries yielded more than 135 Gb of raw data (Table 1). After the initial filtering steps, 89% and 79% of these raw sequences were retained for Stacks pipeline downstream steps for 2015-Sample and 2016-Sample, respectively (Table 1).
SNP calling after successive filtering recovered 1137 loci in 52 individuals (11 individuals were filtered from the total 63 due to the low number of reads, resulting in a final dataset of 27 and 25 individuals in 2015 and 2016 Samples, respectively). From the 1137 loci 52 were in Hardy–Weinberg disequilibrium and 257 were detected as outliers (i.e., putatively under selection). From these 257 outliers’ loci, only 50 were successfully annotated against NCBI database and were related to 171 Gene Ontologies (GOs); no information was found for the remaining 207 outlier loci. Most of the annotated sequences were associated with metabolic process (72.22%) when considering the category of biologic process, such as proteolysis (12%), regulation of transcription by RNA polymerase II (10%), positive regulation of gene expression (8%) and negative regulation of transcription (7%) (Figure S1). A small percentage (1%) of sequences has GOs associated with phosphate-containing compound metabolic process. Specifically, the locus 20,670 can be related to organophosphate metabolic process via adenylate cyclase type 2 (Table S1), showing high similarity to sequences from Musca domestica (XP_011296241.1) and Stomoxys calcitrans (XP_013113360.1) (identities 21/22, 95%). The locus 10,178 blasted with the highest similarity to Lucilia cuprina cadherin-related tumor suppressor (KNC26434.1) (identity 22/23, 96%), which is known to be related to Bt-resistance in insects [61,62]. The locus 226,347 was annotated as glutathione S-transferase, with high similarity to sequences from different Diptera species, as Sarcophaga dux (QGJ03700.1) (identity 22/23, 96%), Bactrocera oleae (XP_014086721.1) (identity 23/23, 100%), Ceratitis capitata (XP_004523110.1) (identity 22/23, 96%), Musca domestica (5ZWP_A) (identity 22/23, 96%) and Lucilia cuprina (P42860.2) (identity 22/23, 96%). Glutathione transferase is highly related to insecticide resistance in insects [63,64].
After all the filtering steps described above, the data matrix for downstream population analyses was represented by 828 neutral loci. Mean expected heterozygosities for both 2015- and 2016-Samples were similar, being somewhat lower in 2016-Sample (0.1317 for 2015-Sample and 0.1141 for 2016-Sample) and were slightly higher than the mean observed heterozygosities (0.1113 for 2015-Sample and 0.1044 for 2016-Sample), culminating in small positive FIS estimates (0.0981 for 2015-Sample and 0.0498 for 2016-Sample). Besides the 52 loci with deviations from Hardy–Weinberg equilibrium in both samples (previously excluded), 116 and 94 loci showed deviations in 2015- and 2016-Samples, respectively (data not shown).
AMOVA based on neutral SNPs indicated that populations were only subtly genetically differentiated (FST = 0.01444, p < 0.05). Bayesian inferences of population structure based on neutral SNPs also indicated some population differentiation over time (Figure 1A–C). Delta K from Evanno et al. method [47] indicated K = 3 as the best clustering scenario (Figure 1A), but one of the clusters (presented in yellow) is apparently misleading, since no individuals have a high admixture proportion assigned to it (Figure 1B). Thus, we considered that K = 2 seems to be more plausible since we are dealing with two temporal samplings (Figure 1C). Similarly, the BIC analysis from DAPC resulted in the smallest value with K = 3 (Figure 2A). However, the choice of K is not always clear, considering that clustering programs assign individuals to groups in order to create a caricature of a complex reality [65]. Therefore, although obtaining the smallest value of BIC with K = 3 (results from DAPC analysis considering K = 3 were presented in Figure S2), we considered K = 2 as the best option to explain the distribution of variation in temporal samples of NWS flies (Figure 2B), since the reduction in BIC values between K = 3 and K = 2 was very small. The two clusters when considering K = 2 were mainly determined by 13 loci (Figure S3). Cluster 1 mainly comprised individuals from 2015-Sample, while cluster 2 exclusively comprised individuals from 2016-Sample (Figure 2C). Similar results for DAPC analysis were found when considering the 257 loci putatively under selection combined with the 828 neutral loci (Figure S4), indicating that population genetic structure is not due to selection.

3.2. mtDNA

Concatenated mtDNA sequences resulted in 21 haplotypes for the 63 individuals, from which 13 were unique. The most frequent haplotypes H1 (32/63) and H2 (5/63) are shared between 2015- and 2016-Samples (Table 2). The 2015-Sample is the most diverse in number of haplotypes and both 2015- and 2016-Samples showed high and moderate haplotype diversity estimated by Ĥ, respectively and low nucleotide diversity (π) (Table 2).
Haplotype network based on mtDNA concatenated sequences indicated that 2015- and 2016-Samples differ in haplotypes composition and frequency, although they share the two most frequent haplotypes (Figure 3). Additionally, the 2016-Sample show a star-shaped pattern of haplotype variation which is considered indicative of population expansion.
AMOVA based on mtDNA data revealed genetic structure between 2015- and 2016-Samples (FST = 0.27537, p < 0.001), which indicates that the NWS fly population from Cañas changed considerably from one year to another when the mitochondrial information is considered.

3.3. Microsatellites

Microsatellite loci showed a moderate to high degree of variability, with Ho ranging from 0.355 to 0.833 and He from 0.361 to 0.852 (Table 3). No locus showed significant deviations from Hardy–Weinberg equilibrium after Bonferroni correction (p < 0.006) in the 2015-Sample. Otherwise, four out of eight loci significantly deviated from Hardy–Weinberg equilibrium in 2016-Sample after correction (p < 0.006), with CH15 and CH26 exhibiting heterozygote deficit (i.e., positive FIS) and CH05 and CH14 exhibiting excess of heterozygotes (i.e., negative FIS).
Tests for detecting genotypic linkage disequilibrium (LD) found significant association between pairs of loci only in the 2016-Sample, with 12 cases for p < 0.05 and 8 cases after Bonferroni correction (p < 0.006) (Table S2). AMOVA for microsatellites indicated that populations were subtly differentiated (FST = 0.08385, p < 0.001). Bayesian inferences of population structure for microsatellite data also indicate population differentiation over time, revealing the presence of two genetic clusters (Figure 1D,E).

3.4. Demographic Inferences

Estimates of effective population size (Ne) based on SNPs and microsatellite data diverged about 3.5-fold: estimates based on SNPs indicated around 250 individuals, while those based on microsatellites indicated a Ne of approximately 70 individuals (Table 4).
Neutrality tests based on mtDNA were statistically significant only in one case: Tajima’s D for 2016-Sample (Table 5). The negative value of Tajima’s D is indicative of population expansion. Since only 2016-Sample revealed some signal of population size change, the mismatch distribution analysis was conducted with this population alone (Figure S5). This analysis resulted in a statistically insignificant Raggedness index (rg = 0.2978, p = 1.00), indicating that we cannot reject the null hypothesis of population expansion.

4. Discussion

Single-nucleotide polymorphisms were hereby obtained for the NWS fly (C. hominivorax) using the two-enzyme based genotyping-by-sequencing (GBS) libraries for the first time. The GBS technique [4,5] radically changed the state of the art of population genomics, allowing the simultaneous obtainment of many informative genomic markers of any species for several individuals, at a relatively low cost [2]. This population genomics approach allowed a refinement in population and demographic parameters inference, which otherwise would not be obtained using most regularly employed markers, such as mtDNA sequencing and microsatellites. SNPs are believed to be superior to microsatellites for elucidating historical demography, since they represent an unbiased sampling of genomic variation [66]. Therefore, we can consider that Ne estimates based on SNP data (250 individuals ranging from 211.2 to 287.8) can be more accurate for our NWS fly sampling than microsatellite estimates (70 individuals ranging from 39.5 to 100.5). In spite of that, these results should be considered with caution, since accurate estimation of Ne and population demographic changes in a contemporary time scale (20 generations) is dependent on the number of samples and the targeted number of RADseq (or GBS) loci (i.e., very large SNP data sets, with more than 25,000 SNPs, are necessary to detect population declines) [14], as well as dependent on the proportion of missing data and minor allele frequency (MAF) [67]. Our Ne estimates of approximately 70 and 250 individuals for microsatellite and SNP loci, respectively, is low compared to other insects. Estimates of Ne for Aedes aegypti across 17 geographic localities and time points resulted in values slightly different from the ones we found, considering the same method of Jorde and Ryman [57]: for microsatellite loci, estimated Ne averaged 303.3, ranging from 25.0 to 1181.0, while for SNPs the average was 166.0, ranging from 22.9 to 549.2 [68].
Wounded-host availability, extensive insecticide use by farmers, seasonality and climate are important factors affecting yearly population size of the NWS fly in its southernmost distribution. We consider two likely hypothesis for NWS fly populations dynamics in Uruguay: (i) during winter unfavorable conditions, characterized by cold temperatures and low precipitation, the species migrates to northern latitudes and recolonize the southern region when climate became favorable; alternatively (ii) during the winter the NWS fly populations are reduced, surviving in local refugees like riverine forests and when climate became favorable the populations expand. This cyclic phenomenon of population size changes based on annual seasons is well known for many insect species [69]. However, we were unable to formally test both hypothesis (i.e., migration–recolonization, maintenance–recolonization) due to the absence of samples for additional years and specially samples taken at the beginning and the end of the favorable season.
The three investigated molecular markers, SNPs, mtDNA and microsatellites, indicated temporal genetic structure between two consecutive years for the NWS fly populations from Cañas, Uruguay. However, SNP data variability showed no drastic change from one year to the next: a similar number of loci showed deviations from Hardy–Weinberg equilibrium and mean observed and expected heterozygosity were similar. Besides that, population structure analyses indicated differentiation between years. SNPs are bi-allelic and expected to be less variable than microsatellites [13]. Despite of this, many studies reported sufficient power for SNPs to detect population structure, even in cases where divergence between populations is very low (e.g., [70]). We indeed found small genetic structure between 2015- and 2016-Sample based on SNPs and microsatellites, whereas mtDNA indicates higher structure. That difference suggests differential permanence of males and females from one year to the next in our population, in contrast with data on dispersion of NWS flies [71]. Additionally, random variation due to small population size may result in genetic structure between populations. In fact, populations of NWS flies are constantly controlled at field conditions to diminish their populations. As result, the next year population will be composed by individuals that survived during the anterior year with the addition of individuals that may migrate from neighbor populations.
Although we have not found any important contribution of loci under selection for our NWS fly samples genetic structure, we were still able to annotate some of them as associated with insecticide resistance, which is of interest since it represents a challenge to species control. One of these genes is related to organophosphate metabolic process, considered the main insecticide resistance mechanism for C. hominivorax [72,73]. Another locus is related to Bt-resistance in insects and, despite not being common for C. hominivorax control, it was demonstrated that Bacillus thuringiensis (Bt) isolates have applicability for Calliphoridae species control, including screwworms (Patent US6056953A—Hickle et al. 2000). Possibly, our NWS samples were exposed to Bt in the field, since it is widely used for controlling sympatric herbivorous species, both as transgenic part of the genome of host plants and as sprays [74,75]. We also found a glutathione transferase sequence, a diverse gene family highly associated with insecticide resistance in insects [63,64]. Annotation is fundamental to identify genomic regions under selection since it adds the adaptive process layer to analyses in population genomics studies, leading to a better understand of population structure patterns [76]. Accordingly, loci putatively under selection, many of them annotated as resistance to insecticides functions, are suggested to be the main factors responsible for the genetic structure found in Spodoptera frugiperda (Lepidoptera, Noctuidae) Brazilian populations [10]. Neither mtDNA nor microsatellites could provide such refined information, and future studies with SNPs in NWS fly populations should confirm the potential role of loci under selection in the spatial population genetic structure.
Despite limitations due to the low number of samples, SNPs obtained by GBS from C. hominivorax samples demonstrate to be a potential powerful marker for both spatial and temporal population genomic studies for this species. Additionally, these markers allow the detection of loci putatively under selection that can be associated with insecticide resistance and other important functions for insect pest control.

Supplementary Materials

The following are available online at https://www.mdpi.com/2075-4450/11/8/539/s1, Figure S1: Pie charts of multilevel gene ontology annotation (GO) of loci putatively under selection in each category (cellular component, molecular function and biologic process), Figure S2: DAPC results based on SNP data, Figure S3: A) Contributions of the alleles to the clustering scenario of K = 2 and B) allele frequencies of the 13 loci that are above the threshold of 0.005 in each of the two clusters, Figure S4: DAPC results based on 1085 SNP data (828 neutral loci + 257 outlier loci detected by Lositan), Figure S5: Mismatch distribution analysis for mtDNA sequences from 2016-Sample, Table S1: Gene ontology annotation (GO) of loci putatively under selection, Table S2: Pairwise linkage disequilibrium between microsatellite loci.

Author Contributions

Conceptualization, L.W.B., K.L.S.-B., P.F. and A.M.L.A.-E.; methodology, L.W.B., K.L.S.-B., R.V. and P.F.; formal analysis, L.W.B., K.L.S.-B., R.V. and P.F.; writing—original draft preparation, L.W.B.; writing—review & editing, L.W.B., K.L.S.-B., R.V., P.F. and A.M.L.A.-E.; supervision, A.M.L.A.-E.; project administration, L.W.B.; funding acquisition, K.L.S.-B., R.V. and A.M.L.A.-E. All authors have read and agreed on the published version of the manuscript.

Funding

L.W.B. was supported by a CNPq fellowship (Process # 141777/2014-1). K.L.S.-B. was supported by FAPESP (Process # 2012/16266-7). A.M.L.A.E. was supported by FAPESP (Process # 2015/02079-9). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brazil (CAPES)—Finance Code 001.

Acknowledgments

The authors would like to thank Daniel Lara for assistance in sample collection; Rosangela Rodrigues for technical support with mitochondrial and microsatellite data; Aline da Costa Lima Mores for technical support with GBS.

Conflicts of Interest

The authors declare no conflict of interest.

Data Availability Statement

GBS raw data can be accessed at NCBI SRA Bioproject PRJNA657064 and Biosamples SAMN15814514, SAMN15814515. Mitochondrial DNA sequences were deposited in the GenBank database under the accession numbers: MT875913 - MT875975 for COI, MT875976 - MT876038 for COII and MT876039 - MT876101 for CR.

Additional Information

The manuscript is part of an unpublished dissertation from L.W.B. at the University of Campinas (Brazil).

References

  1. Luikart, G.; England, P.R.; Tallmon, D.; Jordan, S.; Taberlet, P. The power and promise of population genomics: From genotyping to genome typing. Nat. Rev. Genet. 2003, 4, 981–994. [Google Scholar] [CrossRef] [PubMed]
  2. Davey, J.W.; Hohenlohe, P.A.; Etter, P.D.; Boone, J.Q.; Catchen, J.M.; Blaxter, M. Genome-Wide genetic marker discovery and genotyping using next-generation sequencing. Nat. Rev. Genet. 2011, 12, 499–510. [Google Scholar] [CrossRef] [PubMed]
  3. Andrews, K.R.; Good, J.M.; Miller, M.R.; Luikart, G.; Hohenlohe, P.A. Harnessing the power of RADseq for ecological and evolutionary genomics. Nat. Rev. Genet. 2016, 17, 81–92. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Elshire, R.J.; Glaubitz, J.C.; Sun, Q.; Poland, J.; Kawamoto, K.; Buckler, E.S.; Mitchell, S.E. A Robust, Simple Genotyping-by-Sequencing (GBS) Approach for High Diversity Species. PLoS ONE 2011, 6, e19379. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Poland, J.; Brown, P.J.; Sorrells, M.E.; Jannink, J.-L. Development of High-Density Genetic Maps for Barley and Wheat Using a Novel Two-Enzyme Genotyping-by-Sequencing Approach. PLoS ONE 2012, 7, e32253. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Nosil, P.; Gompert, Z.; Farkas, T.E.; Comeault, A.A.; Feder, J.L.; Buerkle, C.A.; Parchman, T.L. Genomic consequences of multiple speciation processes in a stick insect. Proc. R. Soc. B Biol. Sci. 2012, 279, 5058–5065. [Google Scholar] [CrossRef] [Green Version]
  7. Gompert, Z.; Comeault, A.A.; Farkas, T.E.; Feder, J.L.; Parchman, T.L.; Buerkle, C.A.; Nosil, P. Experimental evidence for ecological selection on genome variation in the wild. Ecol. Lett. 2013, 17, 369–379. [Google Scholar] [CrossRef]
  8. Lozier, J.D. Revisiting comparisons of genetic diversity in stable and declining species: Assessing genome-wide polymorphism in North American bumble bees using RAD sequencing. Mol. Ecol. 2014, 23, 788–801. [Google Scholar] [CrossRef]
  9. Silva-Brandão, K.L.; E Silva, O.A.B.N.; Brandão, M.M.; Omoto, C.; Sperling, F.A. Genotyping-By-Sequencing approach indicates geographic distance as the main factor affecting genetic structure and gene flow in Brazilian populations of Grapholita molesta (Lepidoptera, Tortricidae). Evol. Appl. 2015, 8, 476–485. [Google Scholar] [CrossRef]
  10. Silva-Brandão, K.L.; Peruchi, A.; Seraphim, N.; Murad, N.F.; Carvalho, R.A.; Farias, J.R.; Omoto, C.; Cônsoli, F.L.; Figueira, A.; Brandão, M.M. Loci under selection and markers associated with host plant and host-related strains shape the genetic structure of Brazilian populations of Spodoptera frugiperda (Lepidoptera, Noctuidae). PLoS ONE 2018, 13, e0197378. [Google Scholar] [CrossRef]
  11. Anderson, C.; Tay, W.T.; McGaughran, A.; Gordon, K.H.; Walsh, T. Population structure and gene flow in the global pest, Helicoverpa armigera. Mol. Ecol. 2016, 25, 5296–5311. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Wachi, N.; Matsubayashi, K.W.; Maeto, K. Application of next-generation sequencing to the study of non-model insects. Entomol. Sci. 2017, 21, 3–11. [Google Scholar] [CrossRef] [Green Version]
  13. Habel, J.C.; Husemann, M.; Finger, A.; Danley, P.D.; Zachos, F.E. The relevance of time series in molecular ecology and conservation biology. Biol. Rev. 2013, 89, 484–492. [Google Scholar] [CrossRef] [PubMed]
  14. Nunziata, S.O.; Weisrock, D.W. Estimation of contemporary effective population size and population declines using RAD sequence data. Heredity 2017, 120, 196–207. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Nunziata, S.O.; Lance, S.L.; Scott, D.; Lemmon, E.M.; Weisrock, D.W. Genomic data detect corresponding signatures of population size change on an ecological time scale in two salamander species. Mol. Ecol. 2017, 26, 1060–1074. [Google Scholar] [CrossRef] [PubMed]
  16. Owings, C.G.; Spiegelman, C.; Tarone, A.M.; Tomberlin, J.K. Developmental variation among Cochliomyia macellaria Fabricius (Diptera: Calliphoridae) populations from three ecoregions of Texas, USA. Int. J. Leg. Med. 2014, 128, 709–717. [Google Scholar] [CrossRef] [PubMed]
  17. Guimarães, J.H.; Papavero, N.; Prado, Â.P.D. As miíases na região neotropical (identificação, biologia, bibliografia). Rev. Bras. Zool. 1982, 1, 239–416. [Google Scholar] [CrossRef] [Green Version]
  18. Adams, T.S. The Reproductive Physiology of the Screwworm, Cochliomyia Hominivorax (Diptera: Calliphoridae) II. Effect of constant temperatures on oogenesis. J. Med. Entomol. 1979, 15, 484–487. [Google Scholar] [CrossRef]
  19. Wyss, J.H. Screwworm Eradication in the Americas. Ann. N. Y. Acad. Sci. 2000, 916, 186–193. [Google Scholar] [CrossRef]
  20. Vargas-Terán, M.; Hofmann, H.C.; Tweddle, N.E. Impact of Screwworm Eradication Programmes Using the Sterile Insect Technique. In Sterile Insect Technique: Principles and Practice in Area-Wide Integrated Pest Management; Springer: Dordrecht, The Netherlands, 2005; pp. 629–650. [Google Scholar]
  21. Bergamo, L.W.; Fresia, P.; Azeredo-Espin, A.M.L. Phylogeography and Insecticide Resistance in the NWS in South America and Caribbean Regions. In Area-Wide Integrated Pest Management: Development and Field Application; Hendrichs, J., Pereira, R., Vreysen, M.J.B., Eds.; CRC Press: Boca Raton, FL, USA, in press.
  22. Lyra, M.L.; Klaczko, L.B.; Azeredo-Espin, A.M.L. Complex patterns of genetic variability in populations of the New World screwworm fly revealed by mitochondrial DNA markers. Med. Veter. Entomol. 2009, 23, 32–42. [Google Scholar] [CrossRef]
  23. Fresia, P.; Lyra, M.L.; Coronado, A.; De Azeredo-Espin, A.M.L. Genetic structure and demographic history of new world screwworm across its current geographic range. J. Med. Entomol. 2011, 48, 280–290. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Lyra, M.L.; Fresia, P.; Gama, S.; Cristina, J.; Klaczko, L.B.; de Azeredo-Espin, A.M.L. Analysis of Mitochondrial DNA Variability and Genetic Structure in Populations of New World Screwworm Flies (Diptera: Calliphoridae) from Uruguay. J. Med. Entomol. 2005, 42, 589–595. [Google Scholar] [CrossRef] [PubMed]
  25. Torres, T.T.; Lyra, M.L.; Fresia, P.; Azeredo-Espin, A.M.L. Assessing Genetic Variation in New World Screwworm Cochliomyia Hominivorax Populations from Uruguay. In Area-Wide Control of Insect Pests; Vreysen, M.J.B., Robinson, A.S., Hendrichs, J., Eds.; Springer: Dordrecht, The Netherlands, 2007; pp. 183–191. [Google Scholar]
  26. Bergamo, L.W.; Fresia, P.; Azeredo-Espin, A.M.L. 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. PLoS ONE 2015, 10, e0128441. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Bergamo, L.W.; Fresia, P.; Lyra, M.L.; Azeredo-Espin, A.M.L. High Genetic Diversity and No Population Structure of the New World Screwworm Fly Cochliomyia hominivorax (Diptera: Calliphoridae) on a Microgeographic Scale: Implications for Management Units. J. Econ. Entomol. 2018, 111, 2476–2482. [Google Scholar] [CrossRef] [PubMed]
  28. Doyle, J.J.; Doyle, J.L. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 1987, 19, 11–15. [Google Scholar]
  29. Catchen, J.M.; Amores, A.; Hohenlohe, P.; Cresko, W.A.; Postlethwait, J.H. Stacks: Building and Genotyping Loci De Novo From Short-Read Sequences. G3 Genes Genom. Genet. 2011, 1, 171–182. [Google Scholar] [CrossRef] [Green Version]
  30. Catchen, J.; Hohenlohe, P.A.; Bassham, S.; Amores, A.; Cresko, W.A. Stacks: An analysis tool set for population genomics. Mol. Ecol. 2013, 22, 3124–3140. [Google Scholar] [CrossRef] [Green Version]
  31. Raymond, M.; Rousset, F. GENEPOP (Version 1.2): Population Genetics Software for Exact Tests and Ecumenicism. J. Hered. 1995, 86, 248–249. [Google Scholar] [CrossRef]
  32. Rousset, F. genepop’007: A complete re-implementation of the genepop software for Windows and Linux. Mol. Ecol. Resour. 2008, 8, 103–106. [Google Scholar] [CrossRef]
  33. Antao, T.; Lopes, A.; Lopes, R.J.; Beja-Pereira, A.; Luikart, G. LOSITAN: A workbench to detect molecular adaptation based on a Fst-outlier method. BMC Bioinform. 2008, 9, 323. [Google Scholar] [CrossRef] [Green Version]
  34. Beaumont, M.A.; Nichols, R.A. Evaluating loci for use in the genetic analysis of population structure. Proc. R. Soc. B Biol. Sci. 1996, 263, 1619–1626. [Google Scholar] [CrossRef]
  35. Goetz, S.; García-Gómez, J.M.; Terol, J.; Williams, T.D.; Nagaraj, S.H.; Nueda, M.J.; Robles, M.; Talón, M.; Dopazo, J.; Conesa, A. High-Throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008, 36, 3420–3435. [Google Scholar] [CrossRef] [PubMed]
  36. Lischer, H.E.L.; Excoffier, L. PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 2011, 28, 298–299. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Lessinger, A.C.; Azeredo-Espin, A.M.L. Evolution and structural organisation of mitochondrial DNA control region of myiasis-causing flies. Med. Veter. Entomol. 2000, 14, 71–80. [Google Scholar] [CrossRef] [PubMed]
  38. Thompson, J.D.; Gibson, T.J.; Plewniak, F.; Jeanmougin, F.; Higgins, D.G. The CLUSTAL_X windows interface: Flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 1997, 25, 4876–4882. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Tamura, K.; Stecher, G.; Peterson, D.; Filipski, A.; Kumar, S. MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Mol. Biol. Evol. 2013, 30, 2725–2729. [Google Scholar] [CrossRef] [Green Version]
  40. Torres, T.T.; Brondani, R.P.V.; Garcia, J.E.; Azeredo-Espin, A.M.L. Isolation and characterization of microsatellite markers in the new world screw-worm Cochliomyia hominivorax (Diptera: Calliphoridae). Mol. Ecol. Notes 2004, 4, 182–184. [Google Scholar] [CrossRef]
  41. Torres, T.T.; De Azeredo-Espin, A.M.L. Development of new polymorphic microsatellite markers for the New World screw-worm Cochliomyia hominivorax (Diptera: Calliphoridae). Mol. Ecol. Notes 2005, 5, 815–817. [Google Scholar] [CrossRef]
  42. Johathan Liu, C.S.; Hulce, D.; Li, X.; Snyder-Leiby, T. GeneMarker® genotyping software: Tools to increase the statistical power of DNA fragment analysis. J. Biomol. Tech. 2011, 22, S35–S36. [Google Scholar]
  43. Van Oosterhout, C.; Hutchinson, W.F.; Wills, D.P.M.; Shipley, P. Microchecker: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 2004, 4, 535–538. [Google Scholar] [CrossRef]
  44. Excoffier, L.; Lischer, H.E.L. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 2010, 10, 564–567. [Google Scholar] [CrossRef]
  45. Pritchard, J.K.; Stephens, M.; Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 2000, 155, 945–959. [Google Scholar] [PubMed]
  46. Jombart, T.; Devillard, S.; Balloux, F. Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genet. 2010, 11, 94. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Evanno, G.; Regnaut, S.; Goudet, J. Detecting the number of clusters of individuals using the software structure: A simulation study. Mol. Ecol. 2005, 14, 2611–2620. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Earl, D.A.; Vonholdt, B.M. Structure Harvester: A website and program for visualizing Structure output and implementing the Evanno method. Conserv. Genet. Resour. 2011, 4, 359–361. [Google Scholar] [CrossRef]
  49. Jakobsson, M.; Rosenberg, N.A. CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 2007, 23, 1801–1806. [Google Scholar] [CrossRef] [Green Version]
  50. Rosenberg, N.A. Distruct: A program for the graphical display of population structure. Mol. Ecol. Notes 2003, 4, 137–138. [Google Scholar] [CrossRef]
  51. Jombart, T. Adegenet: A R package for the multivariate analysis of genetic markers. Bioinformatics 2008, 24, 1403–1405. [Google Scholar] [CrossRef] [Green Version]
  52. Leigh, J.W.; Bryant, D. Popart: Full-Feature software for haplotype network construction. Method. Ecol. Evol. 2015, 6, 1110–1116. [Google Scholar] [CrossRef]
  53. Clement, M.; Snell, Q.; Walker, P.; Posada, D.; Crandall, K. TCS: Estimating Gene Genealogies. In Parallel and Distributed Processing Symposium, Proceedings International; IEEE: Fort Lauderdale, FL, USA, 2002; p. 184. [Google Scholar]
  54. Dieringer, D.; Schlötterer, C. Microsatellite analyser(MSA): A platform independent analysis tool for large microsatellite data sets. Mol. Ecol. Notes 2003, 3, 167–169. [Google Scholar] [CrossRef]
  55. Peakall, R.; Smouse, P.E. GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research—An update. Bioinformatics 2012, 28, 2537–2539. [Google Scholar] [CrossRef] [Green Version]
  56. Guo, S.W.; Thompson, E.A. Performing the Exact Test of Hardy-Weinberg Proportion for Multiple Alleles. Biometrics 1992, 48, 361. [Google Scholar] [CrossRef] [PubMed]
  57. Jorde, P.E.; Ryman, N. Unbiased Estimator for Genetic Drift and Effective Population Size. Genetics 2007, 177, 927–935. [Google Scholar] [CrossRef] [PubMed]
  58. Do, C.; Waples, R.S.; Peel, D.; Macbeth, G.M.; Tillett, B.J.; Ovenden, J.R. NeEstimatorv2: Re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol. Ecol. Resour. 2013, 14, 209–214. [Google Scholar] [CrossRef] [PubMed]
  59. Tajima, F. The Effect of Change in Population Size on DNA Polymorphism. Genetics 1989, 123, 597–601. [Google Scholar] [PubMed]
  60. Fu, Y.X. Statistical Tests of Neutrality of Mutations Against Population Growth, Hitchhiking and Background Selection. Genetics 1997, 147, 915–925. [Google Scholar]
  61. Bravo, A.; Gill, S.S.; Soberón, M. Mode of action of Bacillus thuringiensis Cry and Cyt toxins and their potential for insect control. Toxicon 2007, 49, 423–435. [Google Scholar] [CrossRef] [Green Version]
  62. Pardo-López, L.; Soberón, M.; Bravo, A. Bacillus thuringiensisinsecticidal three-domain Cry toxins: Mode of action, insect resistance and consequences for crop protection. FEMS Microbiol. Rev. 2013, 37, 3–22. [Google Scholar] [CrossRef] [Green Version]
  63. Kostaropoulos, I.; Papadopoulos, A.I.; Metaxakis, A.; Boukouvala, E.; Papadopoulou-Mourkidou, E. Glutathione S–Transferase in the defence against pyrethroids in insects. Insect Biochem. Mol. Biol. 2001, 31, 313–319. [Google Scholar] [CrossRef]
  64. Fang, S.-M. Insect Glutathione S-Transferase: A Review of Comparative Genomic Studies and Response to Xenobiotics. Bull. Insectol. 2012, 65, 265–271. [Google Scholar]
  65. Jombart, T.; Collins, C. A Tutorial for Discriminant Analysis of Principal Components (DAPC) Using Adegenet 2.0.0; Imperial College London–MRC Centre for Outbreak Analysis and Modelling: London, UK, 2015. [Google Scholar]
  66. Brumfield, R.T.; Beerli, P.; Nickerson, D.A.; Edwards, S.V. The Utility of Single Nucleotide Polymorphisms in Inferences of Population History. Trends Ecol. Evol. 2003, 18, 249–256. [Google Scholar] [CrossRef]
  67. Marandel, F.; Charrier, G.; Lamy, J.B.; Le Cam, S.; Lorance, P.; Trenkel, V.M. Estimating Effective Population Size Using RADseq: Effects of SNP Selection and Sample Size. Ecol. Evol. 2020, 10, 1929–1937. [Google Scholar] [CrossRef] [PubMed]
  68. Saarman, N.P.; Gloria-Soria, A.; Anderson, E.C.; Evans, B.R.; Pless, E.; Cosme, L.V.; Gonzalez-Acosta, C.; Kamgang, B.; Wesson, D.M.; Powell, J.R. Effective Population Sizes of a Major Vector of Human Diseases, Aedes Aegypti. Evol. Appl. 2017, 10, 1031–1039. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  69. Shpak, M.; Wakeley, J.; Garrigan, D.; Lewontin, R.C. A Structured Coalescent Process for Seasonally Fluctuating Populations. Evolut. N. Y. 2009, 64, 1395–1409. [Google Scholar] [CrossRef] [PubMed]
  70. Mesnick, S.L.; Taylor, B.L.; Archer, F.I.; Martien, K.K.; Treviño, S.E.; Hancock-Hanser, B.L.; Medina, S.C.M.; Pease, V.L.; Robertson, K.M.; Straley, J.M.; et al. Sperm whale population structure in the eastern and central North Pacific inferred by the use of single-nucleotide polymorphisms, microsatellites and mitochondrial DNA. Mol. Ecol. Resour. 2011, 11, 278–298. [Google Scholar] [CrossRef] [PubMed]
  71. Mayer, D.G.; Atzeni, M.G. Estimation of Dispersal Distances for Cochliomyia hominivorax (Diptera: Calliphoridae). Environ. Entomol. 1993, 22, 368–374. [Google Scholar] [CrossRef]
  72. De Carvalho, R.A.; Torres, T.T.; De Azeredo-Espin, A.M.L. A survey of mutations in the Cochliomyia hominivorax (Diptera: Calliphoridae) esterase E3 gene associated with organophosphate resistance and the molecular identification of mutant alleles. Veter. Parasitol. 2006, 140, 344–351. [Google Scholar] [CrossRef] [PubMed]
  73. Carvalho, A.R.; De Azeredo-Espin, A.M.; Torres, T.T. Deep sequencing of New World screw-worm transcripts to discover genes involved in insecticide resistance. BMC Genom. 2010, 11, 695. [Google Scholar] [CrossRef] [Green Version]
  74. Ffrench-Constant, R.; Daborn, P.J.; Le Goff, G. The genetics and genomics of insecticide resistance. Trends Genet. 2004, 20, 163–170. [Google Scholar] [CrossRef]
  75. Tabashnik, B.E.; Brévault, T.; Carrière, Y. Insect resistance to Bt crops: Lessons from the first billion acres. Nat. Biotechnol. 2013, 31, 510–521. [Google Scholar] [CrossRef]
  76. Calla, B.; Demkovich, M.; Siegel, J.P.; Gomes Viana, J.P.; Walden, K.K.O.; Robertson, H.M.; Berenbaum, M.R. Selective Sweeps in a Nutshell; The Genomic Footprint of Rapid Insecticide Resistance Evolution in an Insect Introduction. bioRxiv 2019. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Structure results based on single-nucleotide polymorphisms (SNPs) and microsatellite data. (A) DeltaK from Evanno et al. method [47] indicating three clusters as the best K value for SNPs; (B) Structure inference based on SNPs, considering K = 3; (C) Structure inference based on SNPs, considering K = 2; (D) DeltaK from Evanno et al. method [47] indicating two clusters as the best K value for microsatellites; (E) Structure inference based on microsatellites considering K = 2. In (B,C,E), each color represents a genetic group and each vertical line represents an individual (for SNPs there is a total of 52 individuals, while for microsatellites there are 63). Each individual is colored proportionally to its probability of assignment to each genetic group.
Figure 1. Structure results based on single-nucleotide polymorphisms (SNPs) and microsatellite data. (A) DeltaK from Evanno et al. method [47] indicating three clusters as the best K value for SNPs; (B) Structure inference based on SNPs, considering K = 3; (C) Structure inference based on SNPs, considering K = 2; (D) DeltaK from Evanno et al. method [47] indicating two clusters as the best K value for microsatellites; (E) Structure inference based on microsatellites considering K = 2. In (B,C,E), each color represents a genetic group and each vertical line represents an individual (for SNPs there is a total of 52 individuals, while for microsatellites there are 63). Each individual is colored proportionally to its probability of assignment to each genetic group.
Insects 11 00539 g001
Figure 2. Discriminant analysis of principal components (DAPC) results based on SNP data. (A) Graph of Bayesian information criterion (BIC) values versus number of clusters. The lowest BIC value indicates K = 3 as the optimal number of clusters. However, we considered K = 2 as the best option to explain the distribution of variation in temporal samples of New World screwworm (NWS) flies; (B) density plot of the two clusters; (C) graph indicating the population of origin of the individuals from each cluster. Cluster 1 is composed of 27 individuals from 2015-Sample and 13 individuals from 2016-Sample, while cluster 2 has only individuals from 2016-Sample (12 individuals).
Figure 2. Discriminant analysis of principal components (DAPC) results based on SNP data. (A) Graph of Bayesian information criterion (BIC) values versus number of clusters. The lowest BIC value indicates K = 3 as the optimal number of clusters. However, we considered K = 2 as the best option to explain the distribution of variation in temporal samples of New World screwworm (NWS) flies; (B) density plot of the two clusters; (C) graph indicating the population of origin of the individuals from each cluster. Cluster 1 is composed of 27 individuals from 2015-Sample and 13 individuals from 2016-Sample, while cluster 2 has only individuals from 2016-Sample (12 individuals).
Insects 11 00539 g002
Figure 3. Haplotype network from mtDNA concatenated sequences of control region (CR), cytochrome c oxidase subunit I gene (COI) and cytochrome c oxidase subunit II gene (COII). Circles sizes are proportional to the number of sequences representing the haplotype.
Figure 3. Haplotype network from mtDNA concatenated sequences of control region (CR), cytochrome c oxidase subunit I gene (COI) and cytochrome c oxidase subunit II gene (COII). Circles sizes are proportional to the number of sequences representing the haplotype.
Insects 11 00539 g003
Table 1. Sequencing results from 2015-Sample and 2016-Sample GBS libraries.
Table 1. Sequencing results from 2015-Sample and 2016-Sample GBS libraries.
2015-Sample2016-Sample
Total number of sequences145,692,719136,389,952
Reads containing adapter sequence9,343,2606,878,774
Ambiguous barcodes934,98313,656,605
Low quality reads51,524164,258
Ambiguous tags5,134,7077,798,267
Retained reads130,228,245 (89%)107,892,048 (79%)
Table 2. Mitochondrial DNA genetic diversity. N—number of individuals; Nh—number of haplotypes; Ĥ—haplotype diversity; π—nucleotide diversity.
Table 2. Mitochondrial DNA genetic diversity. N—number of individuals; Nh—number of haplotypes; Ĥ—haplotype diversity; π—nucleotide diversity.
SampleNNhHaplotypes (# Individuals)Ĥ (SD)π (SD)
20153219H1 (8), H2 (4), H3 (2), H4 (2), H5 (2), H6 (3), H8 (2), H9 to H17 (1)0.9375 (0.0286)0.004530 (0.002372)
2016318H1 (24), H2 (1), H7 (2), H18 to H21 (1)0.4538 (0.1109)0.001858 (0.001064)
Table 3. Microsatellite genetic diversity for each locus in each NWS fly population. N—sample size; Na—number of alleles per locus; AR—allelic richness; Ho—observed heterozygosity; He—expected heterozygosity; FIS—inbreeding coefficient. Significant departures from Hardy–Weinberg equilibrium were indicated by an asterisk (*) in the value of He (p < 0.006).
Table 3. Microsatellite genetic diversity for each locus in each NWS fly population. N—sample size; Na—number of alleles per locus; AR—allelic richness; Ho—observed heterozygosity; He—expected heterozygosity; FIS—inbreeding coefficient. Significant departures from Hardy–Weinberg equilibrium were indicated by an asterisk (*) in the value of He (p < 0.006).
Sample Loci
CH01CH05CH09CH12CH14CH15CH24CH26
2015-SampleN2332313024212132
(N = 32)Na (AR)7 (7)5 (4.97)5 (5)11 (11)4 (4)8 (7.98)6 (6)11 (10.75)
Ho0.5220.6560.3550.8330.3750.6670.7140.781
He0.6930.6650.3610.8040.5330.8110.7870.815
FIS0.26770.02910.0337−0.01970.31570.20110.11630.0572
2016-SampleN3131313131193130
(N = 31)Na (AR)6 (5.99)6 (6)5 (5)7 (6.97)5 (4.90)5 (5)4 (4)9 (9)
Ho0.5480.6450.6130.7100.6450.5790.6130.833
He0.6770.583 *0.6240.7290.616 *0.778 *0.6700.852 *
FIS0.2056−0.09090.03390.0435−0.03180.28130.10090.0391
Table 4. Ne estimates based on SNP and microsatellite (SSR) data (Jorde and Ryman method [57]).
Table 4. Ne estimates based on SNP and microsatellite (SSR) data (Jorde and Ryman method [57]).
Lowest Allele Frequency Used
0.050.020.010+
Harmonic Mean Sample Size27.12727.227.5
SNPFs0.079690.080430.078110.07739
F′0.036310.037100.034690.03435
Ne
(95% CI, parametric)
247.9
(211.4–287.2)
242.6
(211.2–276.1)
259.5
(233.3–287.0)
262.0
(237.3–287.8)
SSRFs0.177220.171670.170520.17003
F′0.135310.130130.129050.12859
Ne
(95% CI, parametric)
66.5
(39.5–100.5)
69.2
(43.9–100.1)
69.7
(45.1–99.6)
70.0
(46.5–98.2)
Table 5. Neutrality tests for mtDNA and their respective p-values (* indicates significant p-value).
Table 5. Neutrality tests for mtDNA and their respective p-values (* indicates significant p-value).
Tajima’s Dp-ValueFu’s FSp-Value
2015-Sample−0.328330.416−2.770430.16600
2016-Sample−2.38060<0.001 *1.106130.72700

Share and Cite

MDPI and ACS Style

Bergamo, L.W.; Silva-Brandão, K.L.; Vicentini, R.; Fresia, P.; Azeredo-Espin, A.M.L. Genetic Differentiation of a New World Screwworm Fly Population from Uruguay Detected by SNPs, Mitochondrial DNA and Microsatellites in Two Consecutive Years. Insects 2020, 11, 539. https://doi.org/10.3390/insects11080539

AMA Style

Bergamo LW, Silva-Brandão KL, Vicentini R, Fresia P, Azeredo-Espin AML. Genetic Differentiation of a New World Screwworm Fly Population from Uruguay Detected by SNPs, Mitochondrial DNA and Microsatellites in Two Consecutive Years. Insects. 2020; 11(8):539. https://doi.org/10.3390/insects11080539

Chicago/Turabian Style

Bergamo, Luana Walravens, Karina Lucas Silva-Brandão, Renato Vicentini, Pablo Fresia, and Ana Maria Lima Azeredo-Espin. 2020. "Genetic Differentiation of a New World Screwworm Fly Population from Uruguay Detected by SNPs, Mitochondrial DNA and Microsatellites in Two Consecutive Years" Insects 11, no. 8: 539. https://doi.org/10.3390/insects11080539

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop