Characterisation of microsatellite and SNP markers from Miseq and genotyping-by-sequencing data among parapatric Urophora cardui (Tephritidae) populations

GBS

DNA was digested with the restriction enzyme EcoR1. Based on in silico digestion of c. 1 Mb of the tephritid Ceratitis capitata (GenBank assembly accession: GCA_000347755.2) and Bactrocera tryoni (GCA_000695345.1) genomes, we estimated that EcoR1 would have c. 15,000 fragments between 200 and 500 bp. DNA for the GBS library was extracted with the Roche template preparation kit (Roche Diagnostics, Mannheim, Germany). DNA quality and its interaction with GBS library quality (sequencing success) was assessed for extractions of (1) whole individuals cut longitudinally and treated with RNAase (Qiagen, Venlo) after DNA extraction, and (2) for DNA extracted from ultrasonically lysed (Bandelin Sonoplus HD70) tissue pellets without RNAase treatment. GBS adapter preparation followed Elshire et al. (2011). We used three pools with six barcodes (individuals) each (Appendix S1). Individual DNA was digested separately in 20 µl end volume consisting of 2 µl NEB buffer (10×), 1 µl EcoR1 (NEB) and 100 ng/µl DNA (variable µl), and filled up to 20 µl with RNase free water (variable µl). The mixture was gently stirred by flicking, centrifuged for 30 s (8,000 rpm), incubated at 37 °C for 15 min and heat inactivated at 65 °C for 20 min. Ligation was done for each individual separately in 50 µl end volume consisting of 20 µl digest (above), 5 µl NEB ligase buffer (10×), 10 µl barcoded adapter mix (=concentrated stock), 13.5 µl RNase free water, 1.5 µl NEB T4 DNA ligase. Reactions were centrifuged for 30 s (8,000 rpm) and incubated at 22 °C for 60min. The ligase was heat inactivated at 65 °C for 30 min. The ligation products of six individuals (each 10 µl) with different barcodes were pooled (=60 µl), added 250 µl binding buffer (Roche) and eluted in 30 µl RNase free water. Each pool was amplified in 6 parallel PCRs (with one specific TruSeq primer per pool). PCR purification was done on three PCR reactions simultaneously using 3 × 250 µl = 750 µl binding buffer, and eluted in 30 µl RNase free water. Finally, the two 3× purifications á 30 µl were pooled, resulting in 60 µl end-elution product. DNA concentrations were measured with Qubit 3.0 (Thermo Fisher Scientific, Waltham, MA, USA), DNA quality was assessed with the 260/280 ratio using Nanodrop (PEQLAB, Erlangen, Germany). Library validation was done with the Agilent 2,100 bioanalyzer (Santa Clara, CA, USA). Libraries were sequenced as 100 nt paired ends on Illumina Hiseq 2000 (1/5 lane, three pools on six individuals) at the Institute of Molecular Genetics, Mainz University.

GBS pools were de-multiplexed with pyRAD 3.0.3 (Eaton, 2014). The paired end reads were merged with PEAR 0.9.5 (Zhang et al., 2014) before being processed further with pyRAD 3.0.3, using standard settings for paired end reads. Data processing was done on the High Performance Computing (HPC) Mogon-cluster (Mainz University). We assessed the number of recovered loci, unlinked SNPs and parsimony informative sites by varying combinations of the parameters “minimum coverage of individuals” (mincov = 10 and 14), and “maximum number of individuals with heterozygotic sites” (maxSH = 4, 6, 8 and 10).

Genetic analyses

Genetic diversity indices were estimated in three populations situated north (Vildbjerg, 56.187°N, 8.771°E), within (Frøslev, 54.827°N, 9.331°E) and south (Neumünster 54.136°N, 9.907°E) of the transition zone. Microsatellite loci were tested in 19–20 individuals per population. GBS SNP variation was quantified for 6 individuals per population. Cross-species amplification of microsatellite loci was tested in U. stylata (Fabricius 1775) that were sampled in Christiansfeld, Denmark (55.354°N, 9.476°E) and Waldgrehweiler, Germany (49.664°N, 7.737°E). For microsatellites, we tested for deviations from Hardy–Weinberg proportions (per locus and population), linkage disequilibria (LD) within populations, genetic diversity per population (allele number and expected heterozygosity) with Genepop on the web (Raymond & Rousset, 1995; Rousset, 2008). Genetic differentiation, F_ST, was estimated locus-wise and among populations with Arlequin 3.5.2.2 (Excoffier & Lischer, 2010). Microsatellite null alleles, stuttering and large allele dropout were evaluated with Micro-Checker version 2.2.3 (Van Oosterhout et al., 2004). For GBS data, we converted the Structure file of unlinked SNPs generated by pyRAD (Eaton, 2014) into an Arlequin file (PGDspider 2.1.0.1, Lischer & Excoffier, 2012). For both microsatellites and GBS, we assessed population affiliations with Structure 2.3.4 (Pritchard, Stephens & Donnelly, 2000) using the admixture model, both with and without prior population information (10 × 50,000 burn-in and 100,000 iterations, K = 1–5). The highest level of genetic clustering was calculated in Clumpak (Kopelman et al., 2015) according to delta K (Evanno, Regnaut & Goudet, 2005) and probability of K, ln Pr(X|K) (Falush, Stephens & Pritchard, 2003).

GBS

Neither DNA preparation/extraction methods nor treatments with or without RNase affected GBS library quality and sequencing success. Although treatment with RNase significantly improved the mean ratio of 260/280 by approaching ≈1.80 (no RNase = 2.066, with RNase = 1.864; t = 9.59, 16 df, P < 0.001), the mean number of reads per pool was higher, though not significantly, for DNA without RNase treatment (Appendix S2). PyRAD analysis quantified 645,164–2,640,743 reads and 7,615–10,002 loci per individual. The mean depth of clusters per individual with depths greater than 7 was 31.7–88.1. Assembled GBS reads have NCBI accession numbers SAMN06241878–SAMN06241895. The parameter settings “minimum coverage of individuals” (mincov) = 10 and 14 and “maximum number of heterozygote samples” (maxSH) = 4, 6, 8 and 10 identified between 1,177 (14 mincov–4 maxSH) and 2,347 (10 mincov–10 maxSH) unlinked SNPs and 1,750 and 4,469 parsimony informative sites, respectively (Fig. 2). Regression analysis showed that the frequency of polymorphic sites, 0.0034–0.0040, did not increase significantly with the number of reads (t =  − 0.61, df = 16, P = 0.55) whereas expected heterozygosity did (t = 2.46, df = 16, P = 0.03). The latter was caused by a significant positive relationship between the number of reads and the number of loci recovered for analysis, t = 3.59, df = 16, P > 0.01. Post hoc analysis showed that this result was influenced by the three individuals with the lowest number of recovered loci (ca. 200–300 loci less than the other individuals) and implied that a plateau was reached above c. 1,000,000 reads using in our protocol. The summary results are presented in Appendix S3.

Structure analyses were performed on loci recovered with the parameter settings 14–6 (i.e., mincov = 14, maxSH = 6), 14–10, 10–6 and 10–10. The parameter set 14–6 recovered 1,492 unlinked SNPs (the second lowest number of SNPs found). For 14–6, the methods of Evanno, Regnaut & Goudet (2005) and Falush, Stephens & Pritchard (2003) estimated K = 2 (Fig. 3A) and K = 4 (Fig. 3B) whereas membership proportions identified three clusters corresponding to the three populations (Fig. 3C), thus repeating the observed pattern found for microsatellites without inclusion of population prior. Increasing the number of loci (=unlinked SNPs) by reducing the number of recovered individuals/locus (e.g., 10–6) and/or allowing more heterozygotic sites across samples (e.g., 14–10) resulted in an identical number of optimal clusters, K = 3, identified by the three clustering methods (Figs. 3D, 3E, 3G, 3H, 3J and 3K). Genetic differentiation among populations (AMOVA) based on parameter set 14-6 (1,492 unlinked SNPs) included 508 loci with less than 0.05 missing data: F_ST = 0.197, pairwise estimates: F_{ST(north-transition)} = 0.165, F_{ST(north-south)} = 0.208, F_{ST(transition-south)} = 0.220, all P < 0.001. Locus-wise AMOVA based on 701 polymorphic loci found 99 loci (14.1%) with significant differentiation, 0.20 < F_ST < 0.91, P < 0.05). Expected heterozygosity/population across 1,492 loci was: North = 0.317, transition zone = 0.331, South = 0.327.