Ver 3

To validate the phylogenetic trees, we assessed the degree of genetic differentiation among geographical isolates for the targeted genes by calculating an Analysis of Molecular Variance model (AMOVA) and pairwise fixation indices (F_ST statistic) between samples using Arlequin ver 3.5.2.2 (Excoffier and Lischer, 2010 (link)). A global locus by locus AMOVA, with 1000 permutations, was constructed using the gene sequence data from all 131 samples and variations within and between countries were estimated. In addition, a matrix of F_ST values, using 100 permutations and a significance level of 0.05, was measured to indicate the level of population divergence, with values ranging from 0 (no evidence of population divergence) to 1 (completely isolated). P-values in both the AMOVA and pairwise F_ST matrix were calculated to assess the level of significance. For each country, we also tested for positive selection among the gametocyte genes by the codon-based Z-test implemented in MEGA X. The Nei-Gojobori method and a bootstrap procedure of 100 replicates were conducted. To validate the results of Z-test, we further estimated the Tajima’s D statistic using Tajima’s D test and D, D^*, F, and F^* values using Fu and Li’s tests implemented in DnaSP (Rozas et al., 2017 (link)).

Allele frequencies, observed heterozygosity (H_obs), and expected heterozygosity (H_exp) were calculated using Microsoft Excel. Likelihood ratio test of linkage disequilibrium (LD) for the syntenic STR loci was performed with Arlequin ver 3.5^{22 (link)} (No. of permutations = 10,000, Significance level = 0.05).

A coalescence-based simulation was used to detect deviating SNP sites. With these methods, we expected to detect low-level differentiation sites under balanced selection (neutral loci) and high-level differentiation sites under directed selection (differentiation loci). The layering method we used was a modification of the FDIST method that was performed with Arlequin ver 3.5.2.2 [78 (link), 79 (link)].. We used a hierarchical island model with 50,000 simulations to calculate the relationship between the Fst and heterozygosity. A locus with an Fst value higher than the 0.99 limit of the neutral distribution was considered to be a putative outlier under divergent selection [59 (link)]. The remaining loci with nonsignificant Fst values were considered to be neutral SNPs. All procedures reduced bias and maintained highly differentiated loci between ecotypes individuals. We chose Fst values that were higher than the expected neutral distribution as the directional selection sites and Fst values that were lower than the expected neutral distribution as the balanced selection sites [60 (link)].

The Brownsea Island, Furzey Island, and Isle of Wight sequences were aligned to previously published data and used the British populations as defined by Hale et al. (2004). The final alignment has a length of 238 bp with 72 informative variants from 1,016 samples from across Europe (see references in Table 1). The numbers of haplotypes, haplotype diversity, nucleotide diversity, and neutrality tests were calculated using DNAsp (Librado & Rozas, 2009). F_ST and AMOVA calculations were performed using Arlequin ver. 3.5.2.2. (Excoffier & Lischer, 2010). A median‐joining haplotype network was constructed in PopART (Leigh & Bryant, 2015), and the Mantel test was calculated using R software and Ade4 package (Dray & Dufour, 2007).

All sequences used in analyses were aligned using MEGA 10.0.4 software [53 (link)]. The standard parameters of genetic diversity (the number of haplotypes, the number of polymorphic sites, haplotype diversity, nucleotide diversity, random match probability (RMP) and the mean number of pairwise differences (MPD)) were calculated using Arlequin ver. 3.5.2.2 software [54 (link)]. The RMP parameter is used for expressing the probability that two randomly sampled individuals from a population have a matching genotype and is calculated as the sum of square frequencies [55 (link)]. The MPD is a parameter that represents the measure of differences between all pairs of haplotypes in the sample. The same software was used for assessing genetic differentiation among populations by the analysis of molecular variances (AMOVA) and estimating pairwise population and overall F_ST values. The statistical significance of all performed tests was assessed with 10,000 permutations. Multi-dimensional scaling (nonmetric MDS) implemented in PAST 3.25 software was used for the visualisation of the matrix of pairwise population F_ST values [56 ]. Frequencies of different mtDNA haplogroups found in analyzed populations were calculated and used for the PCA to assess the contribution of each haplogroup to the genetic differentiation of the populations. PCA was performed using PAST 3.25 software.

Sequences of the 16S rRNA genes of the studied bacteria were analyzed with the Maximum Likelihood model using the MEGA 6.0 [32 (link)]. To compute the degree of the statistical support for branches in the phylogenetic tree, 1,000 bootstrap replicates were made. For an estimation of the relationships between individual strains of species, represented as nodes, minimum spanning trees (MSTs) were calculated; they consist of nodes connected by edges, the length of which corresponds to the distances between two individuals [33 (link)]. MSTs were calculated on the basis of nucleotide diversity according to the maximum parsimony method with the usage of the Arlequin ver. 3.5.2.2 software package [34 (link)].