The invasive spotted-wing drosophila, Drosophila suzukii (Matsumura, 1931), is native to southeast Asia and has become a global pest of fruits (Asplen et al. Reference Asplen, Anfora, Biondi, Choi, Chu and Daane2015; Tait et al. Reference Tait, Mermer, Stockton, Lee, Avosani and Abrieux2021). Over the past 15 years, it has spread rapidly across North America. In Canada, it was first detected in British Columbia in 2009 and in Quebec and Ontario by 2011 (Thistlewood et al. Reference Thistlewood, Gibson, Gillespie and Fitzpatrick2013) and is now widespread across the country. Current management of D. suzukii populations in Canada relies on cultural practices such as short harvest intervals and the removal of unmarketable fruit combined with chemical insecticide applications, applied at regular intervals of 7–14 days through fruit ripening (Ontario Ministry of Agriculture, Food, and Rural Affairs 2022; British Columbia Ministry of Agriculture and Food 2023). In addition, biological control agents Ganaspis brasiliensis Ihering, 1905 and Leptopilina japonica Novković and Kimura, 2011 (both Hymenoptera: Figitidae) have recently established in British Columbia (Abram et al. Reference Abram, Franklin, Hueppelsheuser, Carrillo, Grove and Eraso2022a), and their role in suppression of D. suzukii remains an active area of research.

Understanding the population genetic structure of D. suzukii is potentially useful for refining existing management strategies because genetically distinct populations may differ in phenotypic traits (e.g., hatch rate, generation time, insecticide resistance, susceptibility to parasitoids, and host plant association) that could impact pest management practices (Rota-Stabelli et al. Reference Rota-Stabelli, Ometto, Tait, Ghirotto, Kaur and Drago2020; Olazcuaga et al. Reference Olazcuaga, Foucaud, Deschamps, Loiseau, Claret and Vedovato2022). Given that current management strategies rely heavily on the use of insecticides for population suppression, selection for resistance is a serious threat (e.g., Gress and Zalom Reference Gress and Zalom2019; Ganjisaffar et al. Reference Ganjisaffar, Gress, Demkovich, Nicola, Chiu and Zalom2022). Knowledge of population size, genetic diversity, and rates of gene flow can help build a more resilient management strategy (Green et al. Reference Green, Stenberg and Lankinen2020); for example, the use of untreated host plant refugia in areas where gene flow is known to occur can help to impede insecticide resistance development (Tabashnik et al. Reference Tabashnik, Brévault and Carrière2013). The genetic structure and invasion history of D. suzukii have been the focus of past investigations (Fraimout et al. Reference Fraimout, Debat, Fellous, Hufbauer, Foucaud and Pudlo2017; Lewald et al. Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021) that each found evidence of multiple D. suzukii introduction events from Asia into the United States of America, with a geographic subdivision of eastern and western populations. However, patterns of population structure remain unknown for Canadian D. suzukii populations.

The success of D. suzukii is partially due to its ability to exploit a wide range of host plants that fruit throughout the growing season in different landscapes (e.g., cultivated cherry (Prunus avium Linnaeus), Himalayan blackberry (Rubus armeniacus Focke), and elderberry (Sambucus racemosa Linnaeus); Kenis et al. Reference Kenis, Tonina, Eschen, van der Sluis, Sancassani and Mori2016; Abram et al. Reference Abram, Franklin, Hueppelsheuser, Carrillo, Grove and Eraso2022a). Under laboratory conditions, adaptation of D. suzukii to different host fruits has been observed over multiple generations (Olazcuaga et al. Reference Olazcuaga, Foucaud, Deschamps, Loiseau, Claret and Vedovato2022), but it is unclear whether similar adaptation could occur with the spatial and temporal patterns of fruit use by D. suzukii under natural conditions. If D. suzukii populations are able to associate with or specialise in specific host plants, host fruit–associated genetic structure may occur at a local scale.

Here, we examined the genetic structure of Canadian D. suzukii using double-digest restriction site–associated DNA sequencing (ddRADseq). First, using specimens from four regions, we asked whether there exists large-scale genetic structuring in D. suzukii populations across the country. Second, using D. suzukii that emerged from fruits collected from several host plant species in British Columbia, we assessed whether there is evidence of local population genetic structuring associated with host plant species.

In 2020, we collected fruit samples infested with D. suzukii eggs, larvae, and pupae from four fruit-growing regions across Canada (Fig. 1; Agassiz, British Columbia; Summerland, British Columbia; Vineland, Ontario; and Kentville, Nova Scotia). We collected ripe fruits from six cultivated or wild host plant species between the beginning of June and early October (Table 1). Following the rearing protocol of Abram et al. (Reference Abram, Wang, Hueppelsheuser, Franklin, Daane and Lee2022b), we stored fruits in indoor locations in ventilated containers (12 × 12 × 8 cm, Ziploc^® Medium Square Containers) at ambient or controlled (24 °C) temperatures without wire mesh lining the container bottom. We checked containers for adult D. suzukii emergence two to three times per week. If found, we mouth aspirated adult flies into vials containing 95% ethanol, which we labelled and stored at –20 °C for later molecular analysis. We randomly selected and sexed 10 flies from each host plant + site combination for genetic analysis. We retained representative voucher specimens from each host plant + site combination at –20 °C at the Agassiz Research and Development Centre (Agriculture and Agri-Food Canada, Agassiz, British Columbia) and deposited vouchers at the Royal British Columbia Museum, Victoria, British Columbia.

Figure 1. Principal component analysis and structure barplots depicting spatial genetic clustering among Drosophila suzukii. STRUCTURE results from the first of 10 subsampled replicates of 24 individuals. K values with the greatest statistical support using ΔK and LnP(K) are presented. Photo of male D. suzukii by Yoichiro Watanabe.

Table 1. Summary of collection locations, latitude and longitude coordinates, host plant species, collection dates, and sex ratios for Drosophila suzukii collected for molecular analysis.

BC, British Columbia; NS, Nova Scotia; ON, Ontario.

We extracted genomic DNA from 80 adult D. suzukii with DNeasy Blood & Tissue Kits (QIAGEN, Hilden, Germany) with the manufacturer’s suggested addition of pancreatic ribonuclease A (RnaseA; 4 μL at 100 mg/mL; QIAGEN) to digest RNA. We eluted DNA into 2 × 50 μL of 56 °C Buffer AE to increase DNA yield and concentration, which we kept at −20 °C until library preparation. PstI-MspI restriction enzyme ddRADseq libraries were prepared by staff at the Molecular Biology Service Unit at the University of Alberta (Edmonton, Alberta, Canada) using the protocol of MacDonald et al. (Reference MacDonald, Dupuis, Davis, Acorn, Nielsen and Sperling2020). The two enzymes cut DNA at select restriction sites across the genome to reduce complexity before sequencing. Single-end, 75 base-pair (bp) reads were generated from 200 ng of DNA per individual on one high-output flow cell of a NextSeq 500 (Illumina, California, United States of America).

We demultiplexed DNA sequence data using the process_radtags programme in Stacks 2, version 2.55 (Rochette et al. Reference Rochette, Rivera-Colón and Catchen2019), discarding reads that failed the Illumina purity filter, had uncalled bases, or had Phred scores below 20 over a 15% sliding window of the read length. Resulting reads were 67 bp long after the index sequences were removed. Due to sequencing error associated with the adaptor sequence in the PstI restriction enzyme site, we used Cutadapt, version 3.4 (Martin Reference Martin2011), to trim an additional 5 bp from the 5′ ends, yielding final read lengths of 62 bp. We aligned reads to the near-chromosome level assembly of the D. suzukii genome (Paris et al. Reference Paris, Boyer, Jaenichen, Wolf, Karageorgi and Green2020) with the Burrows–Wheeler Aligner, version 0.7.17 (Li and Durbin Reference Li and Durbin2009), using the BWA–MEM algorithm. We converted SAM files to BAM format and sorted them with SAMtools, version 1.9 (Li et al. Reference Li, Handsaker, Wysoker, Fennell, Ruan and Homer2009), and then used the ref_map programme in Stacks 2 to call single-nucleotide polymorphisms (SNPs), specifying one population containing all individuals. We retained loci that were present in at least 80% of individuals (Paris et al. Reference Paris, Stevens and Catchen2017), with a minor allele frequency cut-off of 0.05. We completed further filtering in VCFtools, version 0.1.16 (Danecek et al. Reference Danecek, Auton, Abecasis, Albers, Banks and DePristo2011), removing individuals with more than 10% missing data and loci with more than 5% missing data. We thinned SNPs to ensure that no two were within 10 000 bp of one another, minimising physical linkage (MacDonald et al. Reference MacDonald, Dupuis, Davis, Acorn, Nielsen and Sperling2020). Our demultiplexed ddRADseq data are published in fastq format in the National Center for Biotechnology Information Sequence Read Archive under accession number PRJNA995753.

We used three independent analyses to quantify population structure. First, we performed principal component analysis in adegenet, version 2.1.5 (Jombart and Ahmed Reference Jombart and Ahmed2011), and visualised output with ggplot2, version 3.3.5 (Wickham Reference Wickham2016). Then we used Bayesian model–based clustering analysis implemented in STRUCTURE, version 2.3.4 (Pritchard et al. Reference Pritchard, Stephens and Donnelly2000), to assess potential genetic clusters (K). To ensure sample sizes were equal for the global STRUCTURE analysis (Puechmaille Reference Puechmaille2016), we used VCFtools to create 10 subsampled data sets from the global data set, each containing all six individuals from Nova Scotia and six randomly subsampled individuals from Agassiz, Summerland, and Ontario (24 of 74 potential individuals per subsampled data set). Sampling three to six individuals per population permits robust assessment of genetic structuring when using SNP data (Shi et al. Reference Shi, Ayub, Vermeulen, Shao, Zuniga and van der Gaag2010; Nazareno et al. Reference Nazareno, Bemmels, Dick and Lohmann2017; Qu et al. Reference Qu, Liang, Wu, Zhao and Chu2020). For our third step, we calculated pairwise F _ST values (Weir and Cockerham Reference Weir and Cockerham1984) using weir-fst-pop in VCFtools and calculated the number of sites with private alleles among groups using private_alleles in poppr, version 2.9.4 (Kamvar et al. Reference Kamvar, Brooks and Grünwald2015).

For each STRUCTURE analysis, we ran 10 replicates for each K value 1–10, stipulating a burn-in of 100 000 followed by 1 000 000 Markov chain–Monte Carlo repetitions. Using the admixture model, we set the alpha prior to 0.5 based on the greatest expected support for K = 2 in principal component analysis (Wang Reference Wang2017), and we used the locpriors option to better resolve spatial structure across the four collection locations (Porras-Hurtado et al. Reference Porras-Hurtado, Ruiz, Santos, Phillips, Carracedo and Lareu2013). We used StructureSelector (Li and Liu Reference Li and Liu2018) to assess the K values with the greatest statistical support, considering both LnP(K) (Pritchard et al. Reference Pritchard, Stephens and Donnelly2000) and ΔK (Evanno et al. Reference Evanno, Regnaut and Goudet2005), and CLUMPAK, version 1.1, to average the 10 replications of each K value and generate Q-matrices (Kopelman et al. Reference Kopelman, Mayzel, Jakobsson, Rosenberg and Mayrose2015). To assess hierarchical genetic structuring (Vähä et al. Reference Vähä, Erkinaro, Niemela and Primmer2007), we re-ran the individuals that constituted separate genetic clusters in the global analyses through Stacks 2 using the above filtering options and then through principal component analysis and STRUCTURE. We set the alpha value to 1 for each hierarchical data set, and we did not use locpriors for data sets that were derived from a single collection location.

After filtering in VCFtools, the global 74-individual data set consisted of 6766 biallelic SNPs with a mean read depth of 31.38× across 15 502 225 filtered reads and 1.56% total missing data. We found two major genetic clusters among the 74 retained D. suzukii individuals in principal component analysis. The first of the 10 subselected, 24-individual STRUCTURE data sets consisted of 5058 biallelic SNPs with a mean read depth of 32.43× across 3 883 504 filtered reads and 1.44% total missing data. K = 2 was best supported using the Evanno et al. (Reference Evanno, Regnaut and Goudet2005) method, and K = 4 was best supported using the Pritchard et al. (Reference Pritchard, Stephens and Donnelly2000) method, although K = 2 was also highly supported (Supplementary material, Figs. S2 and S3). We found similar STRUCTURE results in all subselected data sets (Supplementary material, Fig. S1). Overall, we found that D. suzukii from Agassiz and Summerland comprised one genetic cluster and that D. suzukii from Ontario and Nova Scotia comprised the other cluster (Fig. 1; Supplementary material, Table S1). Mean pairwise F _ST between the western (“BC”) and eastern (“ON + NS”) clusters was 0.071, and the clusters had 1493 and 217 sites with private alleles, respectively. The “BC” cluster had a mean observed heterozygosity of 0.22 (range 0.15–0.30) with a mean expected heterozygosity of 0.25 (range 0.24–0.25), and the “ON + NS” cluster had a mean observed heterozygosity of 0.24 (range 0.19–0.36) with a mean expected heterozygosity of 0.27 (range 0.27–0.28).

Our finding of two geographic genetic clusters is consistent with the strong population structure identified between western and eastern populations of D. suzukii in the United States of America by Lewald et al. (Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021). Although the differences in population genetic methods used in the current study prevent direct comparison of our SNPs to the genotype likelihoods of Lewald et al. (Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021), it is likely that the western and eastern Canadian populations we identified are congruent with those found in the United States of America. Neither our study nor Lewald et al. (Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021) investigated the genetic structure of D. suzukii populations from central Canada and central United States of America. It is possible that “central” individuals are genetically intermediate to the eastern and western clusters on a hypothetical cline that spreads across North America. However, we hypothesise that they belong to the eastern or western genetic cluster. We did not detect significant genetic differences between D. suzukii from Ontario and those from Nova Scotia despite over 1000 km of overland separation, which is good evidence that individual clusters can remain genetically similar over large parts of their range (e.g., along the coasts of North America; Lewald et al. Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021) without forming a genomic cline.

The final hierarchical data set of D. suzukii from British Columbia contained 7063 biallelic SNPs with a mean read depth of 30.50× across 12 361 221 filtered reads and 1.29% total missing data. We found little evidence of principal component analysis clustering within this data set (Fig. 2; n = 58). In the STRUCTURE data set, three clusters (K = 3) had the greatest statistical support using the Evanno et al. (Reference Evanno, Regnaut and Goudet2005) method, whereas 10 clusters (K = 10) had the greatest support using the Pritchard et al. (Reference Pritchard, Stephens and Donnelly2000) method (Supplementary material, Figs. S6 and S7). Such incongruence between principal component analysis and each K method is expected when there is little genetic clustering in hierarchical analysis (e.g., Nelson et al. Reference Nelson, MacDonald and Sperling2022). Mean overall F _ST between host plant groups was very low, varying from −0.010 to +0.016 (Supplementary material, Table S2). Females had 34 sites with private alleles, whereas males had three. Hierarchical principal component analysis, STRUCTURE, and F _ST analysis found little evidence of genetic clustering within D. suzukii collected in British Columbia, indicating that there is little host-associated, sex-associated, or phenological genetic clustering (Fig. 2; Supplementary material, Figs. S4 and S5), although several individuals that emerged from elderberry or raspberry in Agassiz formed two small outlier groups.

Figure 2. Principal component analysis and structure barplots of 58 Drosophila suzukii reared from fruit host species in Agassiz and Summerland, British Columbia. K values with the greatest statistical support using ΔK and LnP(K) are presented.

Hierarchical principal component analysis and STRUCTURE analysis found no evidence of genetic clustering among the individuals collected in Ontario and Nova Scotia (Supplementary material, Figs. S8–S10; n = 16) or in Summerland and Agassiz when analysed separately (not shown due to similarity with Fig. 2). Under laboratory conditions, Olazcuaga et al. (Reference Olazcuaga, Foucaud, Deschamps, Loiseau, Claret and Vedovato2022) found that there is potential for host plant–associated local selection and adaptation in D. suzukii, but our results show no evidence of this occurring under natural conditions in coastal and interior regions of British Columbia. The likelihood of host plant–associated local selection may be relatively low under natural conditions, perhaps due to (1) low genetic diversity among the founding members of the British Columbia cluster of D. suzukii preventing rapid evolution, (2) frequent spatial co-occurrence of different host plants at small scales (Abram et al. Reference Abram, Franklin, Hueppelsheuser, Carrillo, Grove and Eraso2022a), coupled with dispersal ability of D. suzukii between host plants (Tait et al. Reference Tait, Mermer, Stockton, Lee, Avosani and Abrieux2021), and (3) fruiting phenology of single host plants not extending long enough throughout the season to be the only host supporting D. suzukii during its reproductive period (Thistlewood et al. Reference Thistlewood, Rozema and Acheampong2019). As such, D. suzukii would need to reproduce on at least two different host plants to complete its seasonal cycle. Despite our relatively limited sample size at a few sites during a “snapshot” in time, our results showing a lack of host plant–associated genetic structuring suggest that current pest management decisions are best implemented for British Columbia at the landscape level and provide a baseline against which future potential adaptation and selection of invasive D. suzukii populations in western Canada can be measured.

Our finding of a distinct subdivision of western and eastern Canadian D. suzukii populations is consistent with that of Lewald et al. (Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021), who found similar subdivision of eastern and western D. suzukii populations in the United States of America. These findings could have implications for the management of D. suzukii if the two populations have different phenotypic parameters or different levels of susceptibility to integrated pest management tactics. For example, there may be geographic variation in the susceptibility of these D. suzukii populations to native, adventive, or intentionally introduced natural arthropod enemies (reviewed in Wang et al. Reference Wang, Lee, Daane, Buffington and Hoelmer2020; Abram et al. Reference Abram, Wang, Hueppelsheuser, Franklin, Daane and Lee2022b), as has been found in other Drosophila–parasitoid associations (e.g., Kraaijeveld and Godfray Reference Kraaijeveld and Godfray1999). It is also possible that these two D. suzukii populations respond differently to insecticide treatments or develop resistance (e.g., Gress and Zalom Reference Gress and Zalom2019; Ganjisaffar et al. Reference Ganjisaffar, Gress, Demkovich, Nicola, Chiu and Zalom2022) at different rates. Whether there are phenotypic differences relevant to pest management in western versus eastern D. suzukii populations associated with the genetic structure identified by ourselves and others (Fraimout et al. Reference Fraimout, Debat, Fellous, Hufbauer, Foucaud and Pudlo2017; Lewald et al. Reference Lewald, Abrieux, Wilson, Lee, Conner and Andreazza2021) remains an open question that should be addressed in the coming years.