Spotted-wing drosophila, Drosophila suzukii (Diptera: Drosophilidae), exhibits large-scale spatial genetic structure across Canada but not fruit host–associated genetic structure

Abstract

Spotted-wing drosophila, Drosophila suzukii, is a global pest of soft fruits that is capable of reproducing on a wide range of cultivated and wild plant species. In Canada, D. suzukii was first reported in British Columbia in 2009 and is now widespread across the country. Understanding the genetic structure of D. suzukii populations could be important for pest management if there are phenotypic differences between genetically distinct populations. For example, insect pest populations could respond differently to directional selection imposed by insecticides, differ in their host plant preferences, and vary in their susceptibility to biological control agents. Here, we used double-digest restriction site–associated DNA sequencing to examine large- and fine-scale patterns of the genetic structure of D. suzukii reared from fruit hosts in Canada. We found that this species has a large-scale spatial genetic structure; the flies collected formed two distinct genetic clusters, one of which was distinct to western Canada and the other to eastern Canada. At the local scale, D. suzukii populations showed no evidence of host-associated structuring in British Columbia, suggesting that pest management tactics may be best applied at the landscape level. Our results highlight the need to investigate phenotypic differences between western and eastern D. suzukii populations in Canada.

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. 2015; Tait et al. 2021). 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. 2013) 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. 2022a), 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. 2020; Olazcuaga et al. 2022). 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 2019; Ganjisaffar et al. 2022). Knowledge of population size, genetic diversity, and rates of gene flow can help build a more resilient management strategy (Green et al. 2020); 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. 2013). The genetic structure and invasion history of D. suzukii have been the focus of past investigations (Fraimout et al. 2017; Lewald et al. 2021) 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. 2016; Abram et al. 2022a). Under laboratory conditions, adaptation of D. suzukii to different host fruits has been observed over multiple generations (Olazcuaga et al. 2022), 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. (2022b), 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. (2020). 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. 2019), 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 2011), 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. 2020) with the Burrows–Wheeler Aligner, version 0.7.17 (Li and Durbin 2009), using the BWA–MEM algorithm. We converted SAM files to BAM format and sorted them with SAMtools, version 1.9 (Li et al. 2009), 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. 2017), with a minor allele frequency cut-off of 0.05. We completed further filtering in VCFtools, version 0.1.16 (Danecek et al. 2011), 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. 2020). 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 2011), and visualised output with ggplot2, version 3.3.5 (Wickham 2016). Then we used Bayesian model–based clustering analysis implemented in STRUCTURE, version 2.3.4 (Pritchard et al. 2000), to assess potential genetic clusters (K). To ensure sample sizes were equal for the global STRUCTURE analysis (Puechmaille 2016), 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. 2010; Nazareno et al. 2017; Qu et al. 2020). For our third step, we calculated pairwise F _ST values (Weir and Cockerham 1984) 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. 2015).

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 2017), and we used the locpriors option to better resolve spatial structure across the four collection locations (Porras-Hurtado et al. 2013). We used StructureSelector (Li and Liu 2018) to assess the K values with the greatest statistical support, considering both LnP(K) (Pritchard et al. 2000) and ΔK (Evanno et al. 2005), and CLUMPAK, version 1.1, to average the 10 replications of each K value and generate Q-matrices (Kopelman et al. 2015). To assess hierarchical genetic structuring (Vähä et al. 2007), 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. (2005) method, and K = 4 was best supported using the Pritchard et al. (2000) 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. (2021). 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. (2021), 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. (2021) 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. 2021) 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. (2005) method, whereas 10 clusters (K = 10) had the greatest support using the Pritchard et al. (2000) 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. 2022). 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. (2022) 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. 2022a), coupled with dispersal ability of D. suzukii between host plants (Tait et al. 2021), 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. 2019). 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. (2021), 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. 2020; Abram et al. 2022b), as has been found in other Drosophila–parasitoid associations (e.g., Kraaijeveld and Godfray 1999). It is also possible that these two D. suzukii populations respond differently to insecticide treatments or develop resistance (e.g., Gress and Zalom 2019; Ganjisaffar et al. 2022) 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. 2017; Lewald et al. 2021) remains an open question that should be addressed in the coming years.