Population structure, adaptation and divergence of the meadow spittlebug, Philaenus spumarius (Hemiptera, Aphrophoridae), revealed by genomic and morphological data

Sampling

Adults and nymphs of Philaenus spumarius were collected in 2010 and 2011 from eight populations (Fig. 1) across the distribution range of the species: Cerkes, Anatolia, Turkey (TUR); Mount Parnassus, Greece (GRE); Haapamäki-Keuruu, Finland (FIN); Fitou, South of France (FRAN); Gouveia/Fontanelas, Sintra, Portugal (POR); Aberdare, South Wales, United Kingdom (UK); S. Miguel island, Azores (AZO); Wonalancet, New Hampshire, United States of America (USA). P. tesselatus was sampled in Morocco (MOR), from three main localities: near Azrou, near Rabat and near Ceuta. In total, 170 specimens were sequenced, including 20–22 individuals from each sampling site except Morocco, from where only 7 specimens were included (Table S1). Since sampling in this last location was initially intended for a phylogeographic characterization (Rodrigues et al., 2014), only a few individuals were collected from each site.

Nymphs were hand collected from the spittle masses they produce. Adults were collected by sweeping the vegetation with an entomological net. Both nymphs and adults occur and feed on large numbers of host plant species and no particular hosts were selected during sweeping and hand collection. Efforts by others to show associations between common hosts and the genetically determined color polymorphism gave negative results (Halkka et al., 1967) and we know of no evidence suggesting host-specific genetic differentiation.

Insects were preserved in absolute or 96% ethanol until they were subjected to DNA extraction after up to one year.

Morphological characters

Philaenus species distinction based on the morphology of male genitalia was accomplished for a subset of 38 males from across all populations, except the AZO and FRAN, for which only immature individuals or females were collected. Nine additional males from MOR, POR and TUR were included to increase morphological sample size, but were not used for genetic analyses (Table S1). Preparation and measurements of male genitalia were done as detailed in Supplementary Information 1.

Five variables calculated from the nine measurements (Fig. 2) were used in the morphometric analysis: total length of aedeagus (TotLen), mean length of lower appendages (LowLen), mean length of middle appendages (MidLen), mean length of upper appendages (UpLen), mean curvature of upper appendages (UpCur). The mean value of measurements of paired structures was considered instead of both left and right measurements individually to reduce some of the variability and the number of specimens to be dropped out of the analysis due to missing values related to appendages that were occasionally broken or tilted during aedeagus removal.

A Principal Component Analysis (PCA) was used to evaluate if morphological characters of the aedeagus could separate Philaenus species and/or populations. PCA was applied to standardized variables (centred by the mean and scaled by the variance), since they were measured in different units. Three specimens were left out of the analysis due to missing values. The analysis was performed in R version 3.4.1 (R Core Team, 2017) using the prcomp function and figures were produced using the package ggplot2 version 3.2.1 (Wickham, 2016).

DNA extraction and mitochondrial DNA analyses

DNA was extracted from the head and thorax of each specimen using the DNeasy Blood & Tissue kit (Qiagen) following the manufacturer’s instructions and including an RNase A treatment step. Wings and abdomen were not used for DNA extraction to avoid extracting DNA of endosymbionts and parasites. The obtained DNA was assessed for the presence of a high molecular weight band on the agarose gel after electrophoresis, and it was quantified in Qubit 2.0 (Invitrogen), using Qubit dsDNA HS Assay kit.

A subset of 48 specimens from the nine areas were sequenced for mitochondrial DNA (Table S1) from which we amplified an 800 bp fragment of the 3′-end of the mitochondrial gene cytochrome c oxidase subunit I (COI) by polymerase chain reaction (PCR). Primers used were: C1-J-2195 (5′–TTGATTTTTTGGTCATCCAGAAGT–3′) and TL2-N-3014 (5′–TCCAATGCACTAATCTGCCATATTA–3′) (Simon et al., 1994). PCR was performed in a 12.5 μL reaction volume containing: 1×Colorless GoTaq Flexi Buffer, 2 mM MgCl₂, 0.2 mM dNTPs, 0.6 mg/mL of BSA, 0.5 μM of each primer, 0.0375 U GoTaq DNA Polimerase (Promega) and approximately 30 ng of DNA. PCR conditions were: an initial denaturation step at 95 °C for 5 min, followed by 35 cycles of denaturation at 95 °C for 45 s, annealing at 50 °C for 35 s and extension at 72 °C for 2 min, with a final extension period at 72 °C for 10 min.

Chromatograms were verified and edited using SEQUENCHER v. 4.0.5 (Gene Codes Corporation), they were aligned using CLUSTAL W on BIOEDIT v. 7.0.9 (Thompson, Higgins & Gibson, 1994; Hall, 1999) and subsequently trimmed to the same length. We followed the designation of haplotypes of Rodrigues et al. (2014). A median-joining haplotype network was constructed using POPART version 1.7 (Bandelt, Forster & Rohl, 1999; Leigh & Bryant, 2015).

RAD libraries preparation and sequencing

RAD libraries were prepared using a protocol by Etter et al. (2011), with modifications as in Rodrigues et al. (2016). The restriction enzyme used was SbfI (New England Biolabs). Six libraries were prepared, with 28 to 31 individually barcoded samples multiplexed. The libraries were sequenced on three lanes of an Illumina HiSeq 2000 in paired-end mode (2 × 100 bp) at Genepool (Ashworth Laboratories) (http://genepool.bio.ed.ac.uk/). The individuals from each population were distributed over the different libraries and lanes to avoid library or lane-specific biases.

Assembly and SNP calling

The sequence reads from each run were examined by process_radtags from STACKS version 1.45 (Catchen et al., 2013), to remove those with uncalled bases and low-quality scores (phred score lower than 10), to check that the barcode and restriction site were intact in each read and to demultiplex the samples based on the barcode identification. Reads were trimmed at the 3′end, using TRIMMOMATIC v. 0.38, to keep only 87 bases, since preliminary analyses using the entire read revealed a high number of (possibly false) SNPs at the 3′end after this number of bases (data not shown). This may be due to higher sequencing errors towards the end of the reads characteristic of Illumina sequencing (Dohm et al., 2008). The trimmed reads were de novo assembled into “stacks” (identical sets of reads, called loci) for each individual using the STACKS module ustacks. The minimum depth of coverage to build a stack (-m) was set to 10, the maximum number of nucleotide differences allowed between stacks to form a locus (-M) was set to 2. Then the STACKS module cstacks was used to build a catalog by merging stacks (loci) from multiple individuals, using the default options. The module sstacks was used to match loci from an individual against the catalog. Stacks with very high coverage were removed since they may represent highly repetitive regions and that may include non-orthologous sequences.

Finally, the populations module was used to create a Variant Call Format (VCF) file with the bi-allelic genotypes of each individual for each variable nucleotide position. The minimum number of populations a locus must be present in to process a locus (-p) was set to the number of populations analyzed (eight or nine, excluding or including Morocco, respectively–see below), the minimum percentage of individuals in a population required to process a locus for that population (-r) was set to 0.5 (50%). Only one SNP per locus was kept, using the option –write_random_snp. Other parameters (-M 2, 3 and 6, -m 5 and 10, -n 1 and 4, - -min_maf 0.05 or with no min maf) were tested and the differences were assessed by general patterns in the Principal Component Analysis.

The VCF file was then filtered using VCFTOOLS (v 0.1.14) (Danecek et al., 2011), excluding sites with less than 75% of individuals with genotype (–max-missing 0.75) and/or with minor allele count of 2 (–mac 2), to exclude singletons (Linck & Battey, 2019). In order to exclude overclustered loci, we filtered out those with a mean depth (across individuals) higher than 200×(–max-meanDP 200).

The filtered VCF file with the SNP genotypes was converted into the file formats needed for the different analysis programs using PGDSPIDER 2.0.4.0 (Lischer & Excoffier, 2012). For conversion of GESTE format to BAYPASS format, we used the script https://github.com/CoBiG2/RAD_Tools/blob/master/geste2baypass.py as of commit b99636e.

VCFTOOLS was used to calculate summary statistics of coverage and percentage of missing data. GENETIX v. 4.05.2 was used to obtain expected and observed heterozygosity, as well as F_IS in each population. Pairwise differentiation between populations (F_ST) were calculated in ARLEQUIN 3.5.1.3, and the significance of F_ST was obtained from permutation tests with 10,000 repetitions. Mantel tests between F_ST/(1-F_ST) matrices and the natural logarithm of the geographical distance (Rousset, 1997) were performed with ape package version 5.0 (http://ape-package.ird.fr/) in R version 3.4.0, using 9,999 permutations.

Population structure

Principal Components Analysis (PCA) was used as an exploratory tool of the population structure (Novembre & Stephens, 2008). Computations of PCA were performed in R using package SNPRelate version 1.12.0 (Zheng et al., 2012). Population structure was further examined using the model-based clustering algorithm implemented in STRUCTURE v. 2.3.4 (Falush, Stephens & Pritchard, 2003; Pritchard, Stephens & Donnelly, 2000). We obtained the coefficients of ancestry using the admixture model and assuming correlated allele frequencies among populations, and K from 1 to 9, with 10 replicate runs of each, applying 50,000 steps of burnin and 1,000,000 MCMC steps after burnin. STRUCUTRE_THREADER version 1.2.2 (Pina-Martins et al., 2017) was used to parallelize the runs and to find the K best explaining the data by calculating Delta K on STRUCTURE HARVESTER (Evanno, Regnaut & Goudet, 2005; Earl & vonHoldt, 2012). CLUMPP version 1.1.2 (Jakobsson & Rosenberg, 2007) was then used to obtain the optimal alignment of ancestry proportions, by permuting the 10 replicate runs of STRUCTURE for each value of K.

The complete dataset consisted of 9 populations and 133 individuals. We also analyzed a dataset excluding P. tesselatus individuals from Morocco, which consisted of 8 populations and 127 individuals of P. spumarius. In order to compare morphological variation and genetic variation in P. spumarius, we used a dataset of the 32 individuals for which we had data for both morphometry and RAD-seq. We performed PCA for both types of data, as above, and we calculated Spearman correlations between the Principal Component scores obtained from both PCAs using R.

Detection of selection–outlier analyses and environmental associations

In order to detect loci with signs of selection for the P. spumarius RAD-seq dataset (without Morocco), two approaches were taken: one that detects outlier loci departing from expectation under neutral demographic models (Foll & Gaggiotti, 2008; Vitalis et al., 2014), and another that detects loci associated with environmental variation between populations (Coop et al., 2010; Gautier, 2015).

Outlier analyses were carried out using BAYESCAN v. 2.1 (Foll & Gaggiotti, 2008) and SELESTIM v1.1.4 (Vitalis et al., 2014). BAYESCAN uses a Bayesian approach to estimate the posterior probability of two alternative models for each locus, with or without selection. Posterior odds are then obtained and False Discovery Rate calculated to control for multiple testing. The parameters of the chain and of the model were set to the default values. Outlier SNPs were defined to be those with q-values lower than 5%. SELESTIM v1.1.4 (Vitalis et al., 2014) estimates the intensity of selection at each locus and the posterior distributions of the locus-specific coefficients of selection are compared with a distribution derived from the genome-wide effect of selection using Kullback-Leibler divergence (KLD). KLD is calibrated with simulations from posterior predictive distribution based on observed data (Vitalis et al., 2014). A total of 50 pilot runs of length 1,000 were followed by a run of 1,000,000 with burnin of 10,000. The criterion for a candidate SNP for selection was defined to be the 99% quantile of the KLD distribution.

Environmental and geospatial variables used in the association analysis included 19 bioclimatic variables, as well as longitude and latitude. Bioclimatic variables were mined from WorldClim version 1.4 (release 3) (http://www.worldclim.org/) and the data was extracted for each location using DIVA-GIS 7.5.0 (http://www.diva-gis.org). Associations between SNP allele frequency differences and the environmental variables were assessed with BAYPASS v. 2.1 (Gautier, 2015), using the script Baypass_workflow.R (https://gitlab.com/StuntsPT/pyRona/blob/master/pyRona/R/Baypass_workflow.R) as of pyRONA v0.3.7 (Pina-Martins et al., 2017). Significant associations were assessed with Bayes Factor (BF) obtained with the auxiliary covariate model, considering a threshold for BF of 15. We did not exclude any variable at the start of the study based on their correlations, but we did reassess correlation between significantly associated variables. Spearman correlations between variables were calculated using R.

After finding the candidate SNPs for selection and environmental association, two datasets were created: a “neutral” dataset, for which we excluded the candidate SNPs, and a “candidate” dataset, that contained only the candidate SNPs. STRUCTURE analyses were also performed on these two datasets.

RAD tags with candidate SNPs were queried against the available P. spumarius partial draft genome and transcriptome (Rodrigues et al., 2016), using blastn with an e-value threshold of 1E–30. We obtained a longer sequence (100 bp extended from each end of the RAD tag) from the genome alignment, which was then queried against the NCBI nucleotide database (nr/nt) using BLASTN version 2.9.0 (Altschul et al., 1997), setting a threshold e-value of 1E–5.

Scripts used in these analyses are available at https://github.com/seabrasg/popgenom_Philaenus.git.

Morphology of male aedeagus

The analysis of male genitalia revealed strong differentiation of the three Morocco samples, which showed a characteristic P. tesselatus aedeagus, as originally described in Drosopoulos & Quartau (2002): with the upper appendages longer and weakly curved, extending beyond the lateral appendages and the lower appendages longer and more regularly curved than P. spumarius (Fig. S1). All the remaining samples showed P. spumarius-like aedeagi (Fig. S1). Morphometric analysis confirmed this distinction, with the segregation of Moroccan samples along the first component in the PCA (Figs. 3A and 3B). The variables most associated with this distinction are the mean lengths of upper and of lower appendages of the aedeagus (PCA loadings in Table S2A; boxplots in Fig. S2). The longer appendages in P. tesselatus are expected to be related to the longer body size in general in this species (Drosopoulos & Quartau, 2002) but, when controlling for the total length of the aedeagus, the relative size of the lower appendages remained larger in Moroccan samples (Fig. S2).

Within P. spumarius, we also noted variation across samples, mainly due to the morphometric variables of length and curvature of the upper appendages (Fig. S2; PCA loadings in Table S2B). In particular, there was some geographical structure, for example UK and Finland lying on one extreme of PC2 and Greece lying on the other, corresponding also to the extremes of latitude in this study. A less accentuated curvature and also smaller length of the upper appendages in the Finnish and British than in the Greek samples may be behind this differentiation (boxplots in Fig. S2).

Mitochondrial DNA

The fragment of COI spanned 540 bp and was analyzed for 48 specimens, revealing 25 haplotypes, 21 of which had already been described in Rodrigues et al. (2014). The remaining four (haplotypes UK15, UK18, GR18_13 and FIN9) differed from previously known haplotypes by 1 or 2 substitutions (Fig. S3). Two of these new haplotypes (GR18_13 and FIN9) lay in an intermediate position in the haplotype network between the previously defined “North-Eastern” and “Western” haplogroups. In fact, the three haplogroups are not completely distinct but we maintain their designation in order to more easily describe and visualize the mitochondrial variation in relation to the RAD-seq variation: “Eastern-Mediterranean” (EM) in red, “Western” (W) in green and “North-Eastern” (NE) in blue (Fig. S3). We also attributed similar colors to the groups resulting from RAD-seq for ease of visualization.

All seven specimens from Morocco (MOR) sequenced for mtDNA either showed the most common haplotype of the “W” haplogroup (H29) or a haplotype differing by only one substitution (H28 and H37) (Fig. S3). All haplotypes from the Azores (AZO) and continental Portugal (POR) belonged to “W”. France (FRAN) haplotypes belonged in “W” or in between “W” and “EM” (haplotype H49). Haplotypes from Greece (GRE) belonged to “EM” or in between “W” and “NE” (haplotype GR18_13). Haplotypes from Turkey (TUR) belonged to “NE”. In Finland (FIN), there were haplotypes from “EM”, “NE” and also one between “W” and “NE”. The USA population comprised haplotypes from “NE” and the UK population from “NE” and “W” (Fig. S3). The four new haplotypes were submitted to GenBank under accession numbers MT025773–MT025776.

RAD sequencing

A total of 838,730,936 reads was obtained from the Illumina sequencing. The process_radtags step in STACKS retained 647,870,180 reads. This corresponds to an average of 3,811,001 ± 3,524,799 (standard deviation) reads per individual. Thirty-seven individuals with lower numbers of reads (< 500,000 reads) or large amounts of missing data (> 60%) were excluded from the analysis (Table S1), leaving a total of 133 individuals, for which the number of reads ranged from 736,248 to 23,798,148 (average 4,507,853 ± 3,668,879 sd). Raw reads after demultiplexing were deposited in SRA database with accession PRJNA606428. The population STACKS module, followed by filtering, produced 1,691 SNPs, with a mean coverage of 105.5 reads per locus per individual (Fig. S4) and mean percentage of missing data per individual of 12.3 % (Table S1). We applied relatively stringent filtering criteria to avoid having large amounts of missing data per individual resulting from the large genome size in this species (2.58 Gbases; Rodrigues et al., 2016). This has produced a relatively small number of SNPs but that have a good representation across individuals and that are expected to be scattered across the genome. Since the draft genome is still incomplete and very scattered we were not able to assess this distribution.

Population structure

Principal Component 1 in the PCA clusters Morocco individuals away from the others (Fig. 3C and 3D). When testing other assembly and filter parameters we obtained similar patterns in the groupings of samples (Fig. S5). Also, STRUCTURE analysis gave support to a genetic group solely comprising Moroccan samples (the best K according to Evanno, Regnaut & Goudet (2005) was 4; Fig. 4). The average F_ST of Morocco versus other populations was 0.4, much higher than average F_ST of other populations’ comparisons (0.13) (Table 1). In all population-pairwise F_ST calculations involving Morocco, there were a considerable number of SNPs that were fixed or nearly fixed for one allele in Morocco and for the other allele in all the other populations, as seen in the relatively high frequency of high F_ST values on the histograms in all comparisons and in the high correlations between F_ST values among population pairs (Fig. S6A). There was neither such a high number of fixed SNPs nor such high correlations between F_ST values when considering the other pairs of populations (Fig. S6B). Moroccan samples are thus clearly differentiated, at the genome-wide markers, from the remaining eight populations here analyzed in contrast with mtDNA results that showed no differentiation between Morocco and the Iberian Peninsula (Fig. S3).

The relationship between geographical and genetic distances was significant when considering European populations (excluding from the dataset USA, Azores and Morocco) (Mantel test, z =19.49793, p =0.0177; Fig. 5). When considering the comparisons involving Morocco, a positive correlation is seen between genetic and geographical distances, but mainly because of the lower F_ST value obtained between Morocco and Portugal (F_ST =0.34), than between Morocco and the other populations (F_ST > 0.4; Fig. 5). This lower differentiation may be the result of some level of admixture, which was detected in the STRUCTURE analysis, where all individuals from the Portuguese population show a small contribution from the genetic group present in Morocco (Fig. 4).

Mean diversity (expected heterozygosity, H_E) ranged from 0.0373 (in Morocco) to 0.0808 (in Greece) and mean observed heterozygosity (H_O) from 0.0258 (in Morocco) to 0.0537 (UK). H_O values were generally lower than expected under Hardy-Weinberg equilibrium (HWE) in all populations (average F_IS of 0.346) (Table 1). For Morocco, interpretation of H_E should be carried out carefully, since individuals come from three different locations and thus are not necessarily expected to be in HWE. Additionally, we found a positive and significant correlation between observed heterozygosity and sequence read depth (r_S =0.686, p =0.0412; Fig. S7). This suggests that lower read depths may have led in some cases to allele dropout, contributing towards false homozygotes. However, in the case of Morocco, the mean read depths were not the lowest in this dataset, being similar to others (Fig. S4) and thus this should not be the main factor contributing to the low observed heterozygosity.

The dataset without Moroccan samples consisted of 127 individuals and 2,083 SNPs. For this dataset, PCA revealed two distinct clusters along PC1, generally separating Greece and Turkey from the remaining populations (Fig. 6A). The latter were separated along PC2 in three groups, one including mainly Portugal and France, another including mainly Finland, and a third one, in between these two, including USA, UK and the Azores (Fig. 6A). This structure had already been detected in PC3 of the analysis of the dataset that included Morocco (Fig. 3D). For this dataset, when excluding the USA and the Azores populations, there was again significant isolation-by-distance for the European populations (Mantel test, z =19.49793, p =0.0229). However, there were a few individuals that were genetically more similar to geographically more distant individuals, which is also seen in the STRUCTURE analysis (Fig. 6). In Greece and Finland, there were no admixed individuals between the two main clusters (“Eastern-Mediterranean” in red, and “North-Eastern” in blue)—they were either from one or the other group, with a few exceptions (Fig. 6B). An analysis of the Turkish population revealed the presence of possibly admixed individuals from these two groups (with a smaller contribution from the “North-Eastern” cluster). In the USA and UK populations, all individuals showed some level of admixture between the “North-Eastern” (blue) and the “Western” (green) clusters. The Azores allies a small contribution from the “Western” group to a major one belonging to the “North-Eastern” group and one individual was admixed between “North-Eastern” and “Eastern-Mediterranean” (also seen in PCA). The admixture in USA, UK and the Azores is also apparent from their intermediate position between the Portugal+France group and the Finland+Greece+Turkey group in the PCA (Fig. 6A). The best K in the STRUCTURE analysis, according to the method by Evanno, Regnaut & Goudet (2005), was 3.

The majority of specimens for which COI sequence was available, had a correspondence between the mtDNA and the genomic cluster. However, there were some specimens showing a mismatch consisting of a mtDNA haplotype belonging in a different genomic cluster (Fig. 6B and 6C; Fig. S3). For example, one individual from UK (UK6) bearing a mtDNA haplotype (H24) belonging to “Western” haplogroup (green), turned up “North-Eastern” (blue) in the genome analysis. Two individuals from France, bearing mtDNA H49 haplotype (intermediate between “Western” and “Eastern-Mediterranean”), came up as differentiated at genomic markers, one “Western” (green) and the other intermediate “Western”/“North-Eastern”. In Greece and Finland, COI sequenced individuals show both a mtDNA and genomic makeup belonging to either “Eastern-Mediterranean” or “North-Eastern”, except two individuals assessed as “North-Eastern” in genomic markers but with a mtDNA haplotype in intermediate position in the network, between “Western” and “North-Eastern” haplogroups. The four individuals from Turkey sequenced for COI belonged in “North-Eastern” haplotypes. While one of them had full ancestry from “North-Eastern” group, the other three had their largest proportion of ancestry from “Eastern-Mediterranean”, based on the genome-wide markers.

For samples for which both morphometric and RAD-seq data was available (N = 32), we computed Principal Components Analysis (Fig. S8) and calculated the correlation between PC1 and PC2 scores for both analyses. There was a significant correlation between PC1 from morphometry and PC2 from RAD-seq (r_S =0.63, p = 1E-04), while all the remaining were low and non-significant (r_S = −0.22, p =0.2194 between PC1 from each; r_S =0.1, p =0.5927 between PC2 from each; r_S =0.29, p =0.102 between PC2 from morphometry and PC1 from RAD-seq).

Detection of selection–outlier analysis and environmental associations

Candidate SNPs for positive selection were identified by detection of highly differentiated outliers: eight were detected by BAYESCAN (Fig. S9); and 25 by KLD (quantile 99% KLD 2.037087) in SELESTIM (Fig. S10 ; Table S3). Six outlier SNPs were common to both analyses. No outlier SNPs for balancing selection were detected in the BAYESCAN analysis (Fig. S9).

The BAYPASS analysis detected 163 SNPs associated with environmental variables (BF >15) (Tables S3, S4 and Fig. S11), one of them common to the candidate SNPs detected with Bayescan. Variables showing association were: Longitude, Temperature Annual Range, Precipitation of Driest Quarter, Precipitation of Wettest Quarter, Mean Temperature of Warmest Quarter, Mean Diurnal Range (Mean of monthly (max temp–min temp)). Spearman correlations between these 6 variables were generally low (absolute values below 0.6), with only two values above 0.6 (Table S4C).

When excluding these candidate SNPs (188 in total) from the “full” dataset, creating a “neutral” dataset, the main pattern of structuring was maintained, differing only in admixture proportions at higher values of K (4 and 5) (PCA in Fig. S12 and STRUCTURE in Fig. S13). When analysing the “candidates” dataset, the PCA showed a separation that corresponded generally to longitude variation along PC1 and to latitude along PC2 (Fig. S12). The STRUCTURE analysis, although artificial for the dataset of loci under selection, revealed similar structuring when compared to the other datasets, but with less admixture (Fig. S13). This is an expected outcome considering that these candidate loci have similar allelic variation within each population and different allelic variation between populations. The fact that there are still differentiated individuals within populations in this dataset, consistently assigned to the same groups as in the other datasets, is a reflection of the methods for detecting selection, based on population allelic variation.

Seventy-three candidate SNPs had hits (threshold e-value of 1E–30) with the draft genome of P. spumarius and seven with the transcriptome (Rodrigues et al., 2016). From these, nine had hits (threshold evalue of 1E–5) with predicted genes in the NCBI nucleotide database (Table S3).