Population genomic analysis

Using the SNPs identified from the 389 BMSB individuals, we conducted principal component analysis (PCA), neighbour-net tree and population structure analysis using fastSTRUCTURE to elucidate the population structure and genetic diversity among the BMSB populations. PCA was conducted based on 1775 high quality SNPs obtained from the 389 individuals using Plink 1.9 [57] (parameters used: --allow-extra-chr --pca) and plotted using R package ggplot [32]. The k-mean was 2 as elbow point. To further explore the genetic relatedness among the BMSB individuals, a Minimum Spanning Network (MSN) based on the SNP results of each individual was also reconstructed using R package poppr [33].

To test for the presence of population structure, pairwise F_ST values were generated from the ddRAD sequence data by combining the individual data from the same country as one geographic group. The analyses were implemented in Arlequin 3.5 [58] by converting the SNP profile (in Plink format) to Arlequin format using PGDSpider 2.1.1.5 [59]. The Fst and p values were obtained using 110 permutations in Arlequin 3.5. To visualize the genetic relationships, the obtained population average pairwise F_ST values were further used to construct a neighbour-net method tree using SplitsTree 4 [34].

The AMOVA calculated the variance components and their statistical significance levels for variations among and within the 12 geographic countries and the five genetic clusters. The AMOVA was conducted by using the distance matrix under 1000 premutation in Arlequin 3.5. To further test the genetic variation, the heterozygosity analysis was preformed using GenAlEx 6.5 [60, 61]. The genotype information of all sample was stored into a excel with GenAlEx 6.5 required format [60, 61]. All calculations were conducted based on the default setting.

To provide additional insight into the genetic variation and population differences, population genetic structure analysis was conducted. The population genetic structure was inferred using the SNPs profile of individual samples. In this process, fastSTRUCTURE 1.0 [62] was used to conduct model-based clustering of all individuals. The fastSTRUCTURE 1.0 was run with the default convergence criterion of 10^− 6, a simple prior, and ten replicate runs per dataset. The best K value (i.e. the number of populations or clusters that the samples are best divided into) was determined as 5 using a python script, chooseK.py within the fastSTRUCTURE software.

A roadmap of the BMSB invasive pathway was created using Tableau [63] by inputting a excel file with known BMSB background information. Then, the dotted lines and arrow were added in the map using Mac Preview 11.0.