Simulated Disperser Analysis: determining the number of loci required to genetically identify dispersers

Ethics statement

The collection of all samples followed the strict guidelines outlined in the ethics application (permit number: 05/011A) approved by The University of New South Wales ethics committee.

Sample collection and DNA methods

We investigated dispersal in Australian invasive starlings (Sturnus vulgaris), a highly vagile species that maintains an extensive distribution across southeast Australia. Samples were taken from three populations previously shown to have low levels of genetic differentiation (Rollins et al., 2009; F_ST = 0.02–0.07): Orange (New South Wales, 33°17′S, 149°06′E, N = 32); Mallala (South Australia, 34°27′S, 138°30′E, N = 32); and Munglinup (Western Australia, 33°42′29″S, 120°51′54″E, N = 30).

DNA was extracted using Gentra PureGene DNA extraction kit (Qiagen) following the manufacturer’s instructions. Microsatellites were developed using next generation sequencing on the GS-FLX 454 platform (Roche, Manheim, Germany) following methods described by Abdelkrim et al. (2009). QDD v 0.9.0.0 Beta (Meglécz et al., 2010) was used to identify microsatellites from the sequencing data and primers were designed using the program PRIMER 3 (Rozen & Skaletsky, 2000). A panel of 20 polymorphic markers was chosen (see Table S1). PCR reactions were multiplexed using universal fluorescently labelled primers (Neilan, Wilton & Jacobs, 1997). A step-down PCR protocol was used for each reaction consisting of ten cycles each at the following annealing temperatures: 70 °C, 64 °C, 58 °C, 54 °C, 50 °C. Reactions containing DNA from the same individual but with differently labelled universal markers (e.g., PET, NED, 6-FAM and VIC) were combined in equivalent amounts, so that all loci for each individual could be multi-loaded into three tubes prior to fragment analysis (see Table S1). Samples were genotyped using an ABI 3730 (Applied Biosystems, Foster City, CA, USA) using GS-500 (Liz) in each capillary as a size standard. Allele sizes were estimated on GENEMAPPER version 3.7 (Applied Biosystems). These data were combined with data from 11 microsatellite loci for the same individuals as detailed in Rollins et al. (2009) giving us a total of 31 loci.

Microsatellite markers were checked for departures from Hardy-Weinberg equilibrium (Arlequin version 3.5.1.2, Excoffier, Laval & Schneider, 2005) and linkage disequilibrium (GenePop version 4.0.10, Rousset, 2008). PIC was calculated for each locus using PICcalc (Nagy et al., 2012). Arlequin was used to calculate Pairwise F_ST values. Pairwise values for Shannon’s mutual information index (I, formerly called ^SH_UA, Sherwin et al., 2017) were calculated in GenAlEx (Peakall & Smouse, 2006) because mutual information is better than F_ST at handling a range of population sizes and dispersal rates (Sherwin et al., 2006; Dewar et al., 2011).

Ninety-four samples were successfully genotyped for 31 microsatellite loci. Of the 20 new species-specific loci developed here, two showed significant deviation from Hardy-Weinberg equilibrium (Svu002 and Svu010) and were removed from further analyses (see Table S3). None of the remaining 18 loci showed significant departures from linkage equilibrium after Bonferroni correction. Data from the other 11 loci previously have been shown to be in Hardy-Weinberg equilibrium in these populations by Rollins et al. (2009).

Simulations

Initially, we used GeneClass2 (Piry et al., 2004) to determine the probability of each individual originating in the population where it was sampled (see method below). One individual from each population was randomly assigned to be a ‘simulated disperser’; individuals with low probabilities of population membership were excluded from being simulated dispersers because they were potentially real dispersers. Simulated dispersers were moved in every possible combination relative to other simulated dispersers creating a total of 27 (3³) treatments (including a treatment of no movement, Table 1).

Loci were ranked by their Polymorphic Information Content (PIC) (see statistical analysis below). Each treatment was analysed using 57 genetic datasets. The first analysis started with the locus with the highest PIC and then consecutively adding the next highest PIC locus until all loci were added (N = 29, including the dataset having all markers; referred to as highest PIC hereafter). The second analysis started with the locus with the lowest PIC and then consecutively adding the next lowest PIC locus until all loci were added (N = 28, because the dataset with all markers was analysed in the previous analysis; referred to as lowest PIC hereafter).

The partial Bayesian method (Rannala & Mountain, 1997) implemented in GeneClass2.0 was used to detect simulated dispersers. We ran 10,000 MCMC simulations per population using the L_h∕L_max likelihood computation (Paetkau et al., 2004). Individuals were considered first-generation dispersers if they had an L_h∕L_max p-value that was below 0.01. An individual can incorrectly be identified as a disperser in two ways: if they are identified as a disperser when they are not one (False Positive), or, when they fail to be identified as a disperser when they are one (False Negative). Once identified, GeneClass2.0 assigns dispersers to the most likely source population.

In our assessment of the ability of GeneClass2.0 to detect simulated migrants, assignment tests were repeated for treatments and individuals, resulting in the non-independence of data. To account for this, we used Generalised Linear Mixed Models (GLMMs), specifying treatment and individual as random effects, and predictor variables (number of loci, highest/lowest PIC loci used, location the simulated disperser was moved to, and I) as fixed effects (Zuur et al., 2009). The response variable ‘dispersal status correctly identified’ was a binary indicator of the ability of GeneClass2.0 to correctly identify whether a simulated disperser was a disperser or not (0 = incorrectly identified, 1 = correctly identified) and was specified as having a binomial distribution of errors with a logit link function. For example, if an individual was identified as a resident in the population from which it was sampled, it was classed as ‘dispersal status correctly identified’. The response variable (dispersal status correctly identified) was modelled as a function of (a) the number of loci, (b) whether the highest/lowest PIC loci were used, (c) the population in which the simulated disperser was located (either its collection locality or one of the other two populations) and (d) the pairwise I value between populations when moved. The number of loci was a discrete numeric variable, while highest / lowest PIC loci used (PIC range = 0.251–0.862) was specified as a categorical variable with two levels, location was specified as a categorical variable where individuals were coded as either being from their collection locality or one of the other two populations (collection locality was used as the reference group with which other groups were compared), and pairwise I was specified as a continuous numeric variable. Model fit was estimated from marginal ($R_{\mathrm{GLMM∕(m)}}^2$) and conditional ($R_{\mathrm{GLMM∕(c)}}^2$) coefficients of determination, following Nakagawa & Schielzeth (2013). $R_{\mathrm{GLMM∕(m)}}^2$ estimates model fit using fixed effects only, while $R_{\mathrm{GLMM∕(c)}}^2$ estimates model fit including both fixed and random effects. By comparing these estimates it is possible to compare the contribution that random effects and fixed effects have on the response variable.

Literature search

To assess how studies have addressed parameters that affect power to detect dispersers (number of individuals, number of loci, levels of genetic differentiation) we used Google Scholar to conduct a literature search (on 19/12/2016) for articles that cited the publication that originally presented the GeneClass2 analysis (Piry et al., 2004). We used the key term ‘dispersal’ to narrow the search and of 1,280 articles identified, we randomly chose 160 articles to assess. If the article did not use GeneClass2.0’s migrant detection analysis we disregarded it; if the article did use this method we collected information on the number of loci used, global F_ST value, and general article meta-data such as author(s), year published, common name and species name. Of the 160 papers that were assessed, 132 matched the selection criteria, including a total of 136 datasets (Table S2).

Generalised Linear Models (GLMs) were used to determine the relationship between the global F_ST values of studies using genetic assignment tests and the number of loci that were used by these studies to identified dispersers. The number of loci used was specified as a continuous response variable, with a Poisson distribution, which was determined graphically and the global F_ST value was a fixed effect.

To demonstrate the continuing relevancy of this work (see De Barba et al., 2016), we searched within the 1,280 articles we identified above for those containing the search term ‘microsatellite’ by year since this technique was published (2005–2016) to determine whether the number of papers using this approach was decreasing over time.

Statistical analysis

Regression models and graphics were generated using the statistical programming environment R (version 2.15.2, R Development Core Team, 2012). GLMMs were generated using the function glmer within the package lme4 (version 0.999999-2, Bates, Maechler & Bolker, 2012). Graphics were generated using the package ggplot2 (version 0.9.3.1, Wickham, 2009).

Sample inclusion and genetic analysis

Populations from Orange (New South Wales) and Mallala (South Australia) showed the lowest levels of differentiation from each other (F_ST = 0.026, P value < 0.001; I = 0.081), while the population from Munglinup (Western Australia) displayed higher levels of genetic differentiation from Orange (F_ST = 0.082, P value ≤ 0.001; I = 0.152) and Mallala (F_ST = 0.081, P value ≤ 0.001; I = 0.139) indicating that Munglinup is the most genetically distinct area sampled.

One of 94 individuals was identified as a potential disperser within the original dataset, which included all loci and no simulated dispersal. We avoided using this individual as a simulated disperser.

Simulations

The simulation results showed that the ability to correctly identify a simulated disperser was strongly associated with an increase in the number of loci used and when the highest PIC loci were used rather than the lowest (Table 2, Fig. 1). For a subset of the simulation data (treatments where simulated dispersers were moved outside their collection locality), I values showed a strong positive relationship with correct identification of simulated dispersers. When a simulated disperser was moved to a location that had a relatively high pairwise I to its collection locality, it was more likely to be correctly identified as a disperser (Table 2). For two of our simulated dispersers, no False Positive errors were identified. However, False Negative errors were common.

Simulated disperser A (SD-A), artificially moved from its collection locality at Munglinup

SD-A was always identified as a resident when located in its collection locality, whether using the highest or lowest PIC loci. We were able to identify SD-A to its collection locality when using the three highest PIC loci available, whereas 11 loci were needed when using those having the lowest PIC values. The minimum number of loci required to identify SD-A to its collection locality 100% of the time did not vary between locations (Fig. 2).

Literature search

The studies we evaluated were published between 2005 and 2016. Studies that implemented first-generation migrant detection in GeneClass2.0 to identify dispersers used on average 11.2 ±5.1 loci. For a subset of these data (N = 72), where global F_ST was available, the average global F_ST was, 0.120 ±0.119. There was no relationship between the global F_ST of these studies and the number of loci used when conducting genetic assignment tests (GLM, t₇₁ =  − 1.312; Fig. 3).

Figure 4 demonstrates that, despite the recent increase in population genetic papers based on next-generation sequencing data, microsatellite data sets continue to be used for contemporary disperser identification at a similar rate over the past eight years (more than 100 papers per year during this period).