Between 2017 and 2020, we established 19 laboratory lines of T. japonicus in a containment facility certified by the Canadian Food Inspection Agency at Agriculture and Agri-Food Canada’s Agassiz Research and Development Centre (Agassiz, British Columbia, Canada). We originally collected progenitors of these lines from wild populations in Switzerland, China, the USA, and Canada between 2009 and 2021 (Fig. 1; Suppl. material 1: table S1). We established 16 iso-female lines by taking a single female from a field-collected parasitized pentatomid egg mass and rearing at least 10 generations of its offspring. We established three additional lines using multiple individuals and reared each as mixed laboratory lines (Suppl. material 1: table S1). For each country of collection, we established two to eight lines. We reared all lines on H. halys egg masses following Wong et al. (2021). Elijah Talamas (Florida Department of Agriculture and Consumer Services) or Francesco Tortorici (Department of Agricultural, Forest and Food Sciences, University of Torino) identified vouchers of each line to species-level (Talamas et al. 2015b; Talamas et al. 2017). We deposited representative vouchers in the Florida State Collection of Arthropods; the Canadian National Collection of Insects, Arachnids, and Nematodes; the Summerland Research and Development Centre arthropod collection; and the Royal British Columbia Museum.

Figure 1.

Collection locations for 19 laboratory lines of Trissolcus japonicus. Inset maps depict sampling regions a the Pacific Northwest of North America b Switzerland, and c Beijing and Hebei provinces, China. Symbol shape depicts geographic area of collection: diamond = Switzerland, triangle = Pacific Northwest of North America, circle = China, square = Ontario, Canada. Square symbols for the two locations in London, Ontario, Canada have been jittered for visibility.

We extracted genomic DNA using DNeasy Blood and Tissue DNA kits (QIAGEN, Hilden, Germany) by following the manufacturer’s protocol, but we added a bovine ribonuclease A treatment (RNaseA, 4 uL at 100 mg/mL, QIAGEN) to digest RNA. We eluted DNA into 2× 50 uL of 56 °C Buffer AE to increase DNA concentration and yield, then we stored DNA at -20 °C until ddRAD library preparation. PstI-MspI library preparation and sequencing were performed by sequencing facility staff as outlined in MacDonald et al. (2020) on an Illumina NextSeq 500 at the University of Alberta (Edmonton, Alberta, Canada). We generated single-end 75 base-pair reads in two separate sequencing runs; we sequenced two samples twice to assess run effects.

Obtaining high-concentration, high-quality DNA from small-bodied organisms is an inherent challenge when preparing DNA libraries (e.g., Andersen et al. 2016; Paspati et al. 2019; Campbell et al. 2020), particularly for hymenopteran parasitoids (Cruaud et al. 2018; Gebiola et al. 2019; Ferguson et al. 2020). However, pooling individuals from female iso-lines has become standard practice in RRS and whole genome analyses of hymenopterans (Leung et al. 2020). We determined that extractions from individual T. japonicus did not yield enough DNA for successful ddRADseq (<<200 ng, TDN pers. obs.), and as such, we pooled 10 diploid female T. japonicus from individual lines for each DNA extraction (hereafter referred to as ‘pools’) to meet the minimum required concentration of DNA. Comparison of pooled close relatives may artificially inflate the perceived ‘true’ genetic distance between wild T. japonicus populations (e.g., Rodríguez-Ramilo and Wang 2012), however, our primary objective was to determine the genetic relatedness among laboratory lines. Because Bayesian model-based clustering can be more easily influenced by close relatives than other methods (Waples and Anderson 2017; O’Connell et al. 2019), we employed both Bayesian and non-Bayesian analysis to assess genetic structuring. We prepared a minimum of five pools per line (Table 1), however, we could only prepare four pools from the Ontario, Canada 2021 line due to low specimen availability when sampling. Each pool acted as one unit in analyses.

We used Stacks 2 version 2.55 (Rochette et al. 2019) to demultiplex and process raw DNA reads. We removed reads if they 1) contained Phred scores below 20 over 15% of their length, 2) failed the Illumina chastity filter, or 3) had uncalled bases. We used the ‘barcode rescue’ option to retain reads with one mismatched base in its 8 base adaptor sequence, then we removed all adaptor sequences. Due to some sequencing error in the PstI restriction site, we removed an additional 5 bases from the 5’ end of each read using Cutadapt version 3.4 (Martin 2011), resulting in final lengths of 62 bases. We used the denovo_map pipeline in Stacks 2 to call single nucleotide polymorphisms (SNPs), specifying one ‘population’ in the popmap and a minor allele frequency of 0.05. In accordance with Paris et al. (2017), we retained SNPs that were present across at least 80% of the pools. We completed final filtering using VCFtools version 0.1.16 (Danecek et al. 2011), retaining SNPs with a minimum read depth of five and discarding SNPs and pools that had more than 10% missing data, and we used the thin option to keep only one SNP per stack.

We assessed population genetic structuring for three datasets: one containing all 19 laboratory lines (‘full dataset’), one with seven lines chosen for proportional representation of potential geographic clusters (‘geographic dataset’), and one with all six lines from the Pacific Northwest of North America (‘Pacific Northwest dataset’) (Table 1). We used principal component analysis (PCA) and Bayesian model-based clustering for assessment. We performed PCA in adegenet version 2.1.5 (Jombart and Ahmed 2011) and we visualised results in ggplot version 3.3.5 (Wickham 2016) using R version 4.1.2 (R Core Team, 2021). We implemented Bayesian model-based clustering in structure version 2.3.4 (Pritchard et al. 2000) using the admixture model. For structure analysis, we randomly sub-selected pools from lines to ensure equal sample size from each (Puechmaille 2016), then ran a burn-in period of 100,000 sweeps followed by 1,000,000 sweeps for each of 10 replications for each potential subpopulation (K) between 1–20 (full dataset) or 1–15 (geographic and Pacific Northwest datasets). We defined each laboratory line in a given dataset as a ‘location’ prior to better resolving the genetic structure (Porras-Hurtado et al. 2013). We set the alpha prior of the full and geographic datasets to 1/7, and that of the Pacific Northwest dataset to 1/3, reflecting the expected results of K=7 or K=3 (Wang 2017). We assessed statistical support of each K value using LnP(K) (Pritchard et al. 2000) and ΔK (Evanno et al. 2005) in StructureSelector (Li and Liu 2018). Finally, we generated Q-matrices from the 10 replicates of each K for each dataset in CLUMPAK version 1.1 (Kopelman et al. 2015).

DNA sequences are available as fastq files in the National Center for Biotechnology Information Sequence Read Archive (NCBI SRA) as BioProject PRJNA933214.

In total, we sequenced 109 pools of T. japonicus (1090 individuals) from the 19 laboratory lines, resulting in 109,473,577 raw Illumina reads. In a preliminary PCA, two pools had unexpected behaviour. They may have been contaminated during sample preparation or DNA extraction: neither clustered with any other pool, and both were placed intermediate to other pools extracted from the same line and pools from other lines. We removed both before all subsequent analyses. After our final filtering of raw reads, the full dataset contained 1,889 SNPs across 95 pools with a mean SNP read depth of 64.9×, 11,360,339 filtered reads, and 3.82% total missing data. We used VCFtools to re-filter this dataset before running structure, using the same filtering parameters but sub-selecting near-equal sample sizes for pools from 18 lines in accordance with Puechmaille (2016). We did not include the 2021 line from Ontario, Canada in the full structure dataset because it only had two pools after filtering, leaving a total of 18 lines in the full structure dataset. In 16 of these 18 lines, we subsampled four pools; the other 2 lines (Walla Walla and Beijing-United States Department of Agriculture [USDA]) each had only three pools that passed our filters. Our final full structure dataset had 1,860 SNPs across 70 pools with a mean SNP read depth of 66.0×, 8,469,428 filtered reads, and 1.64% total missing data. In PCA, we identified 6–7 genetic clusters in the full dataset, and we found greatest statistical support for K=9 (LnP(K) method) and K=4 (ΔK method) in the full structure dataset (Fig. 2; Suppl. material 1: figs S1, S2).

Figure 2.

Principal component and structure analyses of SNP data from the full datasets comprising 18–19 Trissolcus japonicus laboratory lines. We present structure results with greatest LnP(K) and ΔK statistical support. Symbol shape depicts geographic area of collection: diamond = Switzerland, triangle = Pacific Northwest of North America, circle = China, square = Ontario, Canada (blue square = London, Ontario, Canada 2021 line). Colours behind laboratory line names correspond with geographic genetic clusters (Fig. 3). We present both modes of K=9 across its 10 replicate runs.

In the geographic dataset, we retained 36 pools across the 7 laboratory lines after filtering raw reads. The dataset contained 2,896 SNPs with a mean SNP read depth of 72.9×, 7,498,995 filtered reads, and 2.03% total missing data. We used VCFtools to sub-select four pools per line before running structure, ensuring equal sample size. Our final geographic structure dataset had a mean SNP read depth of 67.4×, 5,386,771 filtered reads, and 2.45% total missing data across 28 pools. We found greatest statistical support for K=11 (LnP(K) method) and K=10 (ΔK method) in this dataset (Fig. 3, Suppl. material 1: figs S3, S4, K=11 structure plot not shown due to similarity with K=10).

Figure 3.

Structure analysis of SNP data from the geographic structure dataset comprising seven Trissolcus japonicus laboratory lines. We present structure results with greatest ΔK statistical support. Symbol shape depicts geographic area of collection: diamond = Switzerland, triangle = Pacific Northwest of North America, circle = China, square = Ontario, Canada.

In the Pacific Northwest dataset, we retained 33 pools across the 6 laboratory lines after filtering raw reads. This dataset contained 1,976 SNPs with a mean SNP read depth of 78.8×, 5,069,119 filtered reads, and 2.15% total missing data. We used VCFtools to sub-select four pools per line before running structure, ensuring equal sample size; however we could only select three high quality pools from the Walla Walla laboratory line. Our final Pacific Northwest structure dataset had a mean SNP read depth of 71.3×, 3,186,750 filtered reads, and 2.83% total missing data across 23 pools. We found greatest statistical support for K=12 (LnP(K) method) and K=3 (ΔK method) in this dataset (Fig. 4; Suppl. material 1: figs S5, S6).

Figure 4.

Principal component and structure analyses of SNP data from six Trissolcus japonicus laboratory lines collected across the Pacific Northwest of North America. We present the structure results with greatest ΔK and LnP(K) statistical support. We present both modes of K=12 across its 10 replicate runs.

Overall, we found T. japonicus lines had the greatest genetic similarity when collected in close geographic proximity (Figs 1, 2), but some lines were exceptional (Figs 2, 3). In both PCA and structure analyses, lines collected across Switzerland were more closely related to each other than lines from other countries, as were those from Beijing and Hebei provinces of China (Fig. 2), but lines from Switzerland and Beijing/Hebei were genetically different from one another (Fig. 2). At K=4, the Heilongjiang line was included in the same genetic cluster as Beijing and Hebei, however, at K=9 it belonged to its own cluster. This could be evidence of isolation by distance across the native range of T. japonicus and evidence of one introduction event in Switzerland that did not originate from the areas of China we sampled. In contrast, Stahl et al. (2019) found that the Ticino, Switzerland and Beijing lines had only one base substitution difference in their partial COI sequences (i.e., high apparent similarity) and inferred that the Swiss populations could have originated from the Beijing area. However, our results suggest that neither population-level differentiation nor patterns of invasion history of T. japonicus can be evaluated using the mitochondrial genome and thus require a nuclear genome-wide investigation.

Among adventive T. japonicus populations in Canada, the lines from London, Ontario and Kelowna, British Columbia were the most similar to populations from China, the parasitoid’s native range, suggesting that these two populations may have originated from an area in proximity to our sampled Chinese populations. The two lines from London, Ontario were more closely related to the Beijing and Hebei cluster than was the Heilongjiang line in the full dataset (PCA and structure plots, Fig. 2), but the Ontario, Beijing/Hebei, and Heilongjiang lines each formed their own genetic cluster when proportionally represented in the geographic structure analysis (Fig. 3). It is likely that the seven Beijing/Hebei lines influenced the full dataset by over-representing the central China genotype, as each line is a genetic pseudoreplicate from the same population (see Puechmaille 2016). Nonetheless, this analysis suggests that the Ontario population could have originated from central China. Similarly, the Kelowna line was closely related to the Beijing and Hebei cluster at K=4 but formed its own cluster at K=9, suggesting it too may have been introduced from central China. Further genotyping of T. japonicus in its native range will be required to confirm the provenance of adventive populations.

We had expected that westernmost lines from Canada and the USA (Langley, Chilliwack, Agassiz, and Vancouver, WA) would be most closely related, as would those from the interior of BC and WA (Kelowna, BC and Walla Walla, WA). Instead, both lines from Washington State are members of the same cluster despite being separated by more than 350 kilometres, providing good evidence that T. japonicus in Walla Walla and Vancouver are either 1) descendants of a single introduction event in Washington or 2) two separate introduction events from the same region (Fig. 4). Likewise, the three lines from western British Columbia (Langley, Chilliwack, and Agassiz) form their own cluster that is genetically dissimilar from those collected in Washington and does not cluster near any other population we sampled. In addition, the newest detection in the Pacific Northwest, Kelowna, forms its own genetic cluster independent of both western BC and Washington State populations, indicating the population in interior British Columbia is likely an independent adventive introduction of unknown provenance and not a dispersion of the other adventive populations of T. japonicus in the Pacific Northwest. Strikingly, the amount of genetic variation within the Pacific Northwest lines alone is greater than populations separated by more than 1000 km in the native Chinese range. This suggests that there have been at least three distinct introductions of T. japonicus into the Pacific Northwest of North America, and that unintentional introductions of T. japonicus from different source areas may be happening with relative frequency. This may be consistent with the introduction history of its host, H. halys, in the Pacific Northwest, which appears to have resulted from multiple introduction events as suggested by the occurrence of at least three mitochondrial haplotypes (Abram et al. 2017). However, more extensive analyses using the same genotyping methods would be needed to adequately compare the invasion histories of H. halys and T. japonicus. In any case, future surveys for genetic and phenotypic variation in T. japonicus should not assume that populations in geographically proximal regions are necessarily the result of spread from adjacent regions.