Geographic Population Genetic Structure of the New World Screwworm, Cochliomyia hominivorax (Diptera: Calliphoridae), Using SNPs

Abstract

The New World screwworm, Cochliomyia hominivorax, (Coquerel 1858) (Diptera: Calliphoridae) is a major parasite of livestock and other warm-blooded animals, including humans, in the Western Hemisphere. It historically ranged from the southern US, through Central America, South America (except Chile and Southern Argentina), and much of the Caribbean. Female C. hominivorax are attracted to open wounds where they oviposit up to 450 eggs per batch next to the wound or periodically body orifices. The larvae hatch and crawl into the wound to feed on host tissue whereby they cause myiasis (infestation with fly larvae). While feeding using mouth hooks on the anterior end, they keep the wound open in order to breathe through posterior spiracles and cause the wound to increase in size. These wounds then become more attractive for additional flies, secondary infections can occur, and without treatment this frequently leads to death. In 1975, the cost to the cattle industry in the U.S. alone was estimated at $60–100 million dollars in annual losses (Bushland 1975), corrected for inflation this is approximately $300–500 million dollars in 2021. C. hominivorax has since been eradicated from the U.S., Mexico, and Central America to the Panama-Colombian border where a permanent barrier zone is maintained by releasing millions of sterile flies weekly (Scott et al. 2017). Eradication was completed through the use of sterile insect technique (SIT) which was pioneered for C. hominivorax by scientists at the U.S. Department of Agriculture (USDA) (Knipling 1957). SIT for screwworms involves treating pupae with gamma radiation which damages the germline cells and thereby sterilizes the fly. The irradiated flies are then released. Sterilized males seek out and mate with wild females which results in no viable offspring being produced. The wild population is then suppressed and after multiple releases the population can be completely eradicated. Currently eradicated areas benefit from an estimated $1.3 billion annually (Vargas-Terán et al. 2005).

In 2016–2017 there was an outbreak of C. hominivorax in the Florida Keys (Skoda et al. 2018). At least 10% of the entire population of endangered Key deer (Odocoileus virginianus clavium) were discovered to have myiasis. These Key deer died or were euthanized as a result (Skoda et al. 2018) and it was estimated that an additional 10% of the population died undetected. The Florida Department of Agriculture and Consumer Services were the first to report the potential outbreak and sent samples to the USDA Animal and Plant Health Inspection Service (APHIS) National Veterinary Services Laboratories (NVSL) which were then identified as C. hominivorax. Rapid response by the USDA Agricultural Research Service (ARS) and USDA-APHIS International Service using quarantine measures and SIT was successful in re-eradicating the area according to response protocols (USDA-APHIS 2018). In only six months after the first report of cases of C. hominivorax in Florida, it was declared re-eradicated March 23, 2017. However, the geographic source of the Florida Keys outbreak remains unknown.

Identifying the source of outbreaks and reconstructing invasions can lead to knowledge of introduction pathways and the development of preventive control methods (Estoup and Guillemaud 2010). Population genetics has aided in outbreak source identification in other species such as the spongy moth Lymantria dispar (Linnaeaus) (Lepidoptera: Erebidae) (Wu et al. 2020), red clover casebearer moth, Coleophora deauratella (Lienig & Zeller) (Lepidoptera: Coleophoridae) (Mori et al. 2016), and the red fire ant Solenopsis invicta (Buren) (Hymenoptera: Formicidae) (Ascunce et al. 2011). Previous attempts to identify the source of the C. hominivorax Florida outbreak using the COI, COII, and EF1α genetic loci were unsuccessful (Dupuis et al. 2018). Nonetheless, it was determined that there was likely one single introduction pathway because between Florida samples there was little genetic variation and low to zero haplotype nucleotide diversity compared to other countries included in their analysis (Dupuis et al. 2018). There have been multiple population genetic studies of C. hominivorax to examine geographic structuring using nuclear and mitochondrial genes (McDonagh et al. 2009, Fresia et al. 2011, Fresia et al. 2013) as well as allozymes (Taylor and Peterson 1994), RAPD-PCR (Skoda et al. 2013), RFLPs (Roehrdanz and Johnson 1988, Taylor et al. 1996, Lyra et al. 2009), AFLPs (Alamalakala et al. 2009) and microsatellites (Griffiths et al. 2009, Torres and Azeredo-Espin 2009, Mastrangelo et al. 2014). However, in some cases, these studies found low geographic resolution or incongruence which points to a complex population structure and the need for molecular markers with higher resolution.

To overcome the limitations of previously used techniques and to increase resolution by using thousands of markers instead of only a few, a screening of the entire genome for single nucleotide polymorphisms (SNPs) using a genotyping by sequencing (GBS) approach was employed to assess the geographic population structuring and develop a reference database. GBS has been successfully used to explore temporal genetic variation in C. hominivorax from Uruguay (Bergamo et al. 2020), but this is the first study to incorporate genome-wide SNPs for C. hominivorax at a large geographic scale.

Methods

Sample Collection

Samples from Trinidad and Tobago were obtained in July of 2019. The Dominican Republic was sampled in March–November 2019. Sites were selected using remote sensing and further investigated upon arrival and local farms were also sampled. Larvae were collected from myiasis cases in cattle, horses, sheep, pigs, and dogs using only 1 larva per host. Adult flies were collected by baiting with liver traps and aerial sweep netting of arriving flies or sticky-traps baited with Swormlure-4. Larvae and adults were placed in individual vials with 95% ethanol and stored at 4°C. Argentina flies were collected as larvae in 2018. Florida outbreak samples were collected between September 2016 and January 2017. Larval, pupal, and adult samples from Jamaica, Brazil, Costa Rica, Panama, Peru, Colombia, and Mexico were collected between 1984–2017 and were kept in colony at the Commission for the Eradication and Prevention of Screwworms (COPEG) in Panama. Additionally, flies from the SIT production strain, J06, were also included. Samples originated from a total of 11 geographic areas (Table 1, Fig. 1).

Fig. 1.

Geographic origins of samples used in sequencing.

Open in new tabDownload slide

DNA Extraction and Sequencing

DNA was extracted from individual larvae, pupae, and adults. For samples from Trinidad and Tobago, Dominican Republic, and Florida, 3^rd instar larvae were dissected to remove the digestive tract to avoid contamination from host tissues. Dissections were completed using a #12 surgical blade scalpel (Feather, NY) and Swiss #5 forceps (Miltex, PA) under a stereomicroscope. Each larva was placed with dorsal side up and an incision was made near the anterior spiracles longitudinally toward the posterior spiracles. Internal organs were carefully separated from body wall using forceps and nondigestive tract organs, as well as integument, were removed and placed in a microfuge tube.

After completing a set of dissections, tissues were immediately homogenized using liquid nitrogen with a pestle in a microtube. Whole insects were also homogenized using liquid nitrogen. Samples were then incubated with lysis buffer, RNase A, proteinase K, and dithiothreitol overnight at 55°C. Next, DNA was isolated using a modified phenol-chloroform method. A mixture of phenol, chloroform, and isoamyl alcohol (25:24:1) was added to the sample and were slowly inverted to mix 25 times, instead of vortexing to avoid shearing. The resulting aqueous layer was removed using wide-bore tips and added to a new tube, after which the process was repeated for a second round and two rounds of chloroform isoamyl alcohol (24:1). After obtaining the final aqueous layer, sodium acetate was used to pellet DNA overnight at –20°C. Resulting pellet was washed with ethanol and dried. Afterwards, DNA was resuspended in Tris HCl buffer and quantified by Qubit (Invitrogen, MA). Samples were divided into preliminary (n = 34) and final (n = 82) batches for a total of 116. Preliminary samples did not have digestive tracts removed and tissue was homogenized using dry ice. Sequenced samples are detailed in Table 1.

All samples were sequenced at the University of Minnesota Genomics Center (UMGC) using genotyping by sequencing (GBS) reduced representation library approach to obtain small nucleotide polymorphisms (SNPs). GBS libraries were created by UMGC using Btgl and TagI restriction enzymes. Illumina primers and individual barcodes were ligated to DNA fragments and fragments were then amplified using PCR. Fragments were size selected and then sequenced across two lanes on Illumina NextSeq 550 platform. Reads were generated as single end fragments of 150bp in length.

Bioinformatic Analyses

Demultiplexing was performed by UMGC including prefiltering, base calling, and uncertainty assessments. After which, 535,444,823 raw reads were obtained from UMGC. All analyses after obtaining raw reads were completed through the Linux command line on the SCINet Ceres cluster. Reads were first trimmed to remove padding sequences (i.e., barcodes) using UMGC script gbstrim. Reads were then filtered to obtain only those with a Phred score of 33 or higher, which is the standard quality control cut-off, using the process radtags function of STACKS v2.55 (Catchen et al. 2013). After which, reads were indexed and aligned to the C. hominivorax genome (Scott et al. 2020) using the Burrows-Wheeler Aligner BWA v0.7.17 (Li and Durbin 2010). The identification and filtering of SNPs was then continued using STACKS. Methods detailed in Rochette and Catchen (2017) were used to optimize parameters. These parameters included using a significance level of 0.05 to call variant sites, 0.05 to call genotypes, and only including loci that were present in 80% of individuals. Hardy-Weinberg equilibrium was tested using PLINK v1.90 and SNPs in disequilibrium were removed. In order to obtain independent SNPs and avoid linkage disequilibrium only one single SNP per locus was included using the write_single_snp function of STACKS.

Pairwise F_ST values were calculated between all populations using the populations function of STACKS. To examine the population structuring two analyses were completed using the Bayesian clustering method STRUCTURE v2.3.4 (Pritchard et al. 2000) and a Discriminant Analysis of Principal Components (DAPC) (Jombart et al. 2010). STRUCTURE is a model-based approach that assigns individuals (indicated by vertical bars) to a population (indicated by color) or multiple populations if there is admixture, and gives the percent membership probability for that individual assuming markers are in Hardy-Weinberg equilibrium. DAPC identifies clusters using a nonmodel-based multivariate approach that does not make any assumptions about marker neutrality. Therefore, STRUCTURE was run on the dataset of only markers in Hardy-Weinberg equilibrium while DAPC was completed using all markers. STRUCTURE was completed with a burn-in of 10,000 followed by 10,000 generations. The value of K (number of clusters) was selected using the posterior probability method according to Pritchard (2009). In order to provide a more detailed view of the genomic variation, the results for a range of K values are shown (Fig. 2). Results were visualized using DISTRUCT v1.1 (Rosenberg 2004). DAPC was completed using the R package adegenet v2.1.5. The optimal K was selected based on the BIC criteria from the manual using the find.clusters function and DAPC was used to plot clusters (Jombart et al. 2010). Genetic diversity was also evaluated after removing non-neutral markers using STACKS by calculating the observed and expected heterozygosity (H_O H_E), nucleotide diversity (π), and the inbreeding coefficient (F_IS). A phylogenetic analysis was completed using the maximum likelihood program Randomized Axelerated Maximum Likelihood, RAxML v8.2.12 (Stamatakis 2014). Using the rapid bootstrap algorithm, bootstrap branch support percentage was calculated using 1,000 bootstraps returning the best-scoring tree. The unrooted tree was visualized using FigTree v1.4.4.

Fig. 2.

Population STRUCTURE graphs showing the assignment probability for each individual fly. The selected value of K = 12 is shown as well as a range of K = 10–12. The y axis of each graph is the membership probability from 0.0–1.0 and shade corresponds to population (color online). Sample origin is shown along the top JAM (Jamaica), PAN (Panama), BRA (Brazil), COS (Costa Rica), COL (Colombia), TTO (Trinidad & Tobago), PER (Peru), FLO (Florida Outbreak), DOM (Dominican Republic), ARG (Argentina), J06 (production strain), and MEX (Mexico).

Open in new tabDownload slide

Results

The total resulting number of neutral and non-neutral SNPs was 17,905 and after removing markers in Hardy-Weinberg disequilibrium, the total resulting number of neutral SNPs was 9,877 (Supp. Table S1). These SNPs were used in the different population genetic analyses and it was found that there is geographic structuring occurring in C. hominivorax as shown by pairwise F_ST values (Table 2), STRUCTURE (Fig. 2), DAPC (Figs. 3 and 4), and RaxML (Fig. 5). Pairwise F_ST values ranged from 0.105 to 0.711 which indicates “moderate” to “very great” genetic differentiation according to Hartl and Clark (1997). The highest pairwise F_ST value was obtained for comparisons of samples from Peru and Brazil. Comparisons of Florida, Trinidad and Tobago, and Dominican Republic all had F_ST values indicating only moderate genetic differentiation. The genetic diversity statistics are shown in Table 3. The inbreeding coefficients were generally higher for colony reared samples compared to wild-caught samples, and the nucleotide diversity was generally lower for colony reared samples compared to wild-caught samples that matches expectations when comparing colony to wild samples.

Fig. 3.

Discriminant Analysis of Principal Components (DAPC) for all samples. Letters represent individual samples corresponding to the locations in the legend.

Open in new tabDownload slide

Fig. 4.

Discriminant Analysis of Principal Components (DAPC) with Argentina and production strain samples excluded for closer examination. Letters represent individual samples corresponding to the locations in the legend.

Open in new tabDownload slide

Fig. 5.

Phylogenetic unrooted maximum likelihood tree via RAxML. Bootstrap support percentages are shown for each branch.

Open in new tabDownload slide

Results of the STRUCTURE analysis showed a K = 12 as having the highest posterior probability with the marginal likelihood increasing from K = 1–11 and decreasing at K = 12 (Pritchard 2009). But in order to show a more complete picture of the genomic differentiation values of K = 10–12 are shown (Fig. 2). At a K of 10 there appears to be a Central/North America cluster consisting of Mexico, Costa Rica, and Panama but these three populations can be differentiated at a K of 11. While samples from Jamaica, Brazil, Panama, Peru, and Colombia are consistently assigned as their own population clusters throughout. The STRUCTURE graph shows the Florida outbreak samples as their own genotype but somewhat similar to Trinidad & Tobago, and the Dominican Republic.

The first DAPC (Fig. 3) shows Argentina and the Production Strain as the most differentiated and the rest of the populations as one cluster. In order to better examine the other samples, a second DAPC was performed excluding Argentina and the Production Strain (Fig. 4). In the second DAPC, Florida, Dominican Republic, and Trinidad & Tobago are clustered together. Also, Jamaica, and some of the Panama and Colombia samples are clustered together. All Brazil samples cluster independently as well as the Peru and Costa Rica samples. Lastly, the RaxML shows consistent geographic clusters (Fig. 5). Most branches have bootstrap support percentages of 100% according to 1,000 rapid bootstraps, with two branches having a support value of 60% and 80%.

Discussion

In this study, we found twelve populations structured according to geography for C. hominivorax. Our results suggest a cluster we term as the “outer Caribbean population” including Dominican Republic and Trinidad and Tobago bordering the Atlantic Ocean rather than inside the Caribbean Sea. Samples from Jamaica, inside the Caribbean Sea, were genetically distinct from the outer Caribbean and all other populations with pairwise F_ST values ranging between 0.276–0.685 (Table 1). Jamaica also consistently clustered as its own population in the STRUCTURE graphs K = 10–12 (Fig. 2). In the DAPC Jamaica clustered with Panama samples. The Panama samples consist of two colonies originating from Panama (Table 1). These two colonies can be seen differentiated in the DAPC as well as the STRUCTURE graph. Samples from Mexico, Costa Rica, and Panama appear as a Central/North America population in the STRUCTURE analysis at a population size of 10. These three locations also cluster together on the right side of the phylogeny (Fig. 5). These samples can be differentiated independently at a cluster size of 11.

Previous population genetic studies in C. hominivorax have shown that there is high genetic variability in this species. But if, and how, this variability is partitioned according to geography has been debated due to conflicting results obtained using different molecular markers. Discordance between studies using different molecular markers has been observed previously. For example, the American lobster (Homarus americanus) (Milne-Edwards) (Decapoda: Nephropidae) showed virtually no population structuring when it was examined using allozymes and RAPD markers, but using thousands of SNPs allowed for the identification of fine-scale geographic structuring (Benestan et al. 2015). The results of our study are consistent with previous results using RFLPs in mtDNA (Taylor et al. 1996) wherein three assemblages of mitochondrial haplotypes were found, North and Central America, Jamaica, and South America. We also found concordant results with a study using microsatellites (Torres and Azeredo-Espin 2009) that showed population structuring for Caribbean islands but, similar to our study, found that Dominican Republic samples had the lowest pairwise F_ST values (i.e., most similar) when comparing to the Port of Spain, Trinidad and Tobago at 0.1455 and this comparison for our study was at 0.117.

Our results differed compared to Lyra et al. (2009) who found that flies from the Caribbean islands were their own independent entities but South America had comparatively low geographic structure, whereas the results of the current study show geographic structuring in South America. When comparing Colombia and Brazil we found “very great” differentiation (F_ST = 0.285). This supports findings from Fresia et al. (2011) using COI and COII wherein it was found that there were four main groups Cuba, Dominican Republic, and North of Amazon region, and South of Amazon region. This was further explored in Fresia et al. (2013) wherein a north to south colonization was proposed for the continental Americas with two splits occurring. First, at the end of the Last Glacial Maximum, between North/Central America and South American populations and then secondly a split between the Pleistocene and Holocene eras between the North and South Amazonian populations. Mastrangelo et al. (2014) found that the barrier must exist North of the Amazon basin as their samples in northern Brazil matched the South of the Amazon group. Our results also support this grouping as our Brazilian samples came from Goias, Brazil which would belong to the South of Amazon, and they were genetically distinct from Colombia which can be seen in all STRUCTURE models in Fig. 2.

A limitation of this study is the inclusion of flies from laboratory colonies (Peru, Brazil, Jamaica, and Colombia) in addition to wild-caught specimen (Argentina, Florida Outbreak, Trinidad and Tobago, and Dominican Republic). However, this provided the unique opportunity to examine the genetic diversity of historic eradicated populations (Mexico, Costa Rica, and Panama). Although these flies have undergone multiple generations in the laboratory, they were maintained as isolated and distinct origin-specific colonies. Thus, the purpose of this study was a proof-of-concept to determine if the use of GBS would allow for the identification of more detailed geographic structuring with higher resolution and indeed, this was obtained. This might help explain the exceptionally high F_ST values (e.g., 0.71 comparing two laboratory strains) due to changes in allele frequency related to the adaptation and genetic drift associated with laboratory colonization. However, high F_ST values were also obtained for wild flies in Fresia et al. (2011) comparing Cuba and Brazil with an F_ST of 0.85. Additionally, clustering was not only obtained for laboratory flies but significant clustering was also evident in the wild flies which indicates that GBS methods used was valid for both sample types. These results warrant the future inclusion of more wild-caught specimens to build this geographic reference database. This would also allow for the observation of the degree of differentiation that occurs in this species for laboratory colony strains when compared to wild populations from those locations. This could also have applications for the control of this species to know the level of genetic drift that is occurring in the production facility.

The patterns of population structure seen in this study could potentially be explained by one of two mechanisms, the movement of livestock, or flight dispersal in combination with wind patterns. Livestock was introduced into the Americas about 500 years ago and the anthropogenic migration of livestock has been proposed as a driver of current C. hominivorax distribution patterns (Fresia et al. 2011). The phylogeographic patterns for this fly show a population expansion during its dispersal towards its current range, which matched the general pattern of human Native American growth and expansion (Fresia et al. 2013). Therefore, it is plausible that the exportation of livestock or companion animal with myiasis from the outer Caribbean could have brought the flies to the Florida Keys. A similar introduction occurred in France in 1982 when a dog that had recently arrived from Brazil was found to be infested with C. hominivorax (Chermette 1989). Additionally, an outbreak of C. hominivorax occurred in Libya in 1988. Although the origin is unknown, it was suggested the flies probably came from a shipment of contaminated livestock (Lindquist et al. 1992).

The second potential mechanism, this fly’s flight capabilities and global wind patterns, could also explain the pattern of gene flow found in this study for the outer Caribbean population. C. hominivorax commonly flies a distance of 50 km (Mayer and Atzeni 1993) but has been documented to disperse as much as 290 km in its lifetime (Hightower et al. 1965). In addition, the prevailing Northeast trade winds travel from east to west (Gillespie et al. 2012) and ocean currents are affected by wind patterns. Jamaica is part of the Gulf stream while the outer Caribbean is part of the Antilles current. This might explain the stronger genetic differentiation found for Jamaica compared to the increased gene flow found between the outer Caribbean population. Unfortunately, we were not able to obtain flies from Haiti or Cuba which according to this hypothesis, might also cluster as the outer Caribbean population. Thus, we suggest future studies should include samples from these locations.

Here we show a successful proof-of-concept study in the usefulness of SNPs to elucidate fine patterns of geographic structuring in C. hominivorax. The advantage of characterizing population structure using SNPs is a finer-scale resolution with a broader distribution across the genome. The current study demonstrates the ability to identify geographic population structuring comprising at least twelve populations. However, with the inclusion of colony-reared samples, conclusions as to the origin of the Florida outbreak cannot be made at this time. The Florida outbreak samples clustered along with samples from the outer Caribbean population which consists of the Dominican Republic and Trinidad and Tobago. Next, we intend to complete a study using only wild-caught flies and the inclusion of more geographic locations. These results allowed for the initial creation of a geographic database that can be improved upon in order to reference past outbreaks as well as mitigate future outbreaks.

Acknowledgments

The authors would like to thank those that aided in sampling from APHIS-IS Paula Morales, Mario Ambrosino, Renita Sewsaran, and Sallyann Henry. Also, those that helped sampling in the Dominican Republic: Lissette Gomez (COV), Carlos Sanchez, Joel Sanchez, Robinson González, Luis Almánzar, Élida Ramos, Gavino García, Yeliana del Orbe, Teodoro Vásquez, Tania Galarza, and Corporino Nova (DIGEGA-Ministerio de Agricultura). We are also grateful to those that helped with sampling in in Trinidad and Tobago including Jevan Sieusankar, Nicola De Leon, Niger Lowhar, Annmarie Hosein, Cherlyn Ann Wharwood and Victoria Lashley, Michael Downes, and Kern Johnson. Additionally, Keith George, Akeishia Inniss-Bacchus, Paul Crooks, Annika Gordon-Dillon, Sandra Rodriguez, Kate Yeates, the Ministry Personnel, the Veterinary staff of the Ministry of Agriculture, Land and Fisheries, Animal and Production and Health Division, Port of Spain and the Veterinary staff of the Division of Agriculture, and The Tobago House of Assembly, Trinidad and Tobago. Lastly, we would like to thank John B. Welch (APHIS) and Paul V. Hickner (USDA-ARS) for insight on this study and comments on manuscript. This research used resources provided by the SCINet project of the USDA Agricultural Research Service, ARS project 3094-32000-041-00D. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the USDA. This research was supported in part by the U.S. Department of Agriculture, Agricultural Research Service. The U.S. Department of Agriculture prohibits discrimination in all its programs and activities on the basis of race, color, national origin, age, disability, and where applicable, sex, marital status, familial status, parental status, religion, sexual orientation, genetic information, political beliefs, reprisal, or because all or part of an individual’s income is derived from any public assistance program. (Not all prohibited bases apply to all programs.) Persons with disabilities who require alternative means for communication of program information (Braille, large print, audiotape, etc.) should contact USDA’s TARGET Center at (202) 720-2600 (voice and TDD). To file a complaint of discrimination, write to USDA, Director, Office of Civil Rights, 1400 Independence Avenue, S.W., Washington, D.C.20250-9410, or call (800) 795-3272 (voice) or (202) 720-6382 (TDD). USDA is an equal opportunity provider and employer.

Data Availability

Sequence Read Archive (SRA)/ BioProject accession number: PRJNA775625.

Disclaimer: This article reports the results of research only. Mention of a proprietary product does not constitute an endorsement or a recommendation by the USDA for its use. The USDA is an equal opportunity provider and employer.

References Cited