Computational protocol: Tracing the origin of the early wet‐season Anopheles coluzzii in the Sahel

Similar protocols

Protocol publication

[…] Quality of SNP scores was carefully evaluated and 619 SNPs were excluded because of poor scoring or a low genotype conversion rate across mosquitoes (<80%). To simplify concurrent analysis of males and females, the relatively few X‐linked loci (N = 65) were removed from the analysis. Following evaluation of A. coluzzii samples only, a further 40 SNPs were excluded because of extreme deviation from Hardy–Weinberg proportions (using a strict Bonferroni‐corrected threshold of α = 0.05/N SNPs), and a further 74 because of extremely low polymorphism (major allele frequency, MAF > 99%). This yielded a final set of 738 SNPs which were included in the present analysis. Similarly 14 mosquitoes (spread across samples: Kolimana [4], and the Thierola samples from August 2008, May 2009, October 2009 and January 2010 [3, 1, 4 and 2, respectively]), with over 10% genotyping failure rate were also removed, leaving a total sample size of 510 mosquitoes.Within‐sample genetic diversity was measured as the number of polymorphic SNPs (major allele frequency <95%); expected heterozygosity (H e) per locus; and mean number of alleles per locus. To accommodate variation in sample size, these statistics were also computed based on the smallest sample size (n = 40 chromosomes), selecting the first 20 mosquitoes in larger samples. The differences between samples (n = 40 A. coluzzii) in the mean heterozygosity and the mean number of alleles were evaluated using ANOVA, treating locus as a random variable and sample as a fixed parameter and using contrasts to estimate differences between samples according to a priori hypotheses (see predictions in Table ). The analysis used Proc Mixed (SAS 9.4). Exact goodness‐of‐fit tests were used to assess deviations of genotype composition from Hardy–Weinberg expectations (HWE) in each locus and sample, accommodating multiple tests as discussed below. The inbreeding coefficient, Fis, measuring the difference of observed heterozygosity from H e in units of H e (Nei, ; Wright, ) were used to describe such deviations.The Euclidian distance (SQRT(sum[MAFi − MAFj]2), where MAF denotes major allele frequency in the pooled A. coluzzii sample, while the subscripts i and j denote different samples) across loci was estimated as distance between samples. Unlike Fst, this genetic distance remains meaningful not only in comparisons between populations but also between samples representing different time points of the same population. The distance matrix among all sample pairs computed by Proc Distance (SAS) was used to generate a dendrogram using an UPGMA tree building algorithms as well as 2D plots produced by subjecting the distance matrix to multidimensional scaling analysis (MDS). Neighbour joining population trees were plotted using MEGA 4 (Tamura, Dudley, Nei, & Kumar, ) based on the distance matrices described above.To identify subgroups without prior information on species, geographical location or time point, we used STRUCTURE 2.3.3 (Pritchard et al. ). An admixture model assuming correlated allele frequencies was used throughout. Preliminary runs with the possible number of subpopulations (clusters, K) set between 1 and 10 employed 20,000 iterations, the first 10,000 as burn‐in. Probability plots suggested a maximum of K = 6 or 7, and also a rapid plateau in parameter values within the burn‐in period. Consequently, these were followed by test runs of 50,000 iterations, again with a burn‐in of 10,000, with K set between 1 and 7. Each value of K was replicated five times, very low variability in log probability values detected among runs following processing suggesting that additional replicates would not alter the optimum solution. Results were processed with Structure Harvester v0.6.94 (Earl & Vonholdt, ). The optimal value was determined using the deltaK method (Evanno, Regnaut, & Goudet, ). Subpopulation membership of individuals was determined by the highest ancestry proportion (Q) of each individual.To test whether genetic drift between samples collected during the DS (and immediately after the first rains) has occurred, we partitioned the temporal variation in allele frequency (F T, Krimbas and Tsakas ; Nei & Tajima, ; Pollak ) to its components: the experimental sampling variance and the genetic variance (Ne). Thus, if F T was larger than the experimental sampling variance, genetic drift was measurable, consistent with reproduction over one or more generations. On the other hand, if F T was smaller than (or equal to) the experimenter's sampling variance, genetic drift was not evident, consistent with no reproduction. The genetic component is more difficult to detect the larger Ne is (Nei & Tajima, ; Waples, ). Importantly, during the DS and immediately after the rains, density and thus Ne are thought to be at their lowest; hence, genetic drift should be evident if reproduction takes place. Under reproductive arrest (aestivation), no drift can be measured because no sampling of gametes took place, according to our prediction 1.4 (Table ). Instead, drift should be detected during the RS, when density and thus Ne are large yet reproduction is ongoing, facilitating the biological sampling of gametes between generations. While it might appear theoretically easier to detect drift during the DS, when density is the lowest, it should be impossible if there is no reproduction. On the other hand, under migration, temporal samples reflect reproductive cycles either at the origin or in situ, as well as founder effects, further increasing drift. For each pair of samples, F T and the experimental sampling variance (V T) were computed according to Nei and Tajima (): F T = (x i − y i)2/[(x i  + y i)/2 − (x i*y i)] and V T = (1/S xi) + (1/S yi), where x i and y i are the common allele in locus i (all di‐allelic) in the first (x) and second (y) samples, and S xi and S yi are the sample sizes (number of chromosomes) of the first (x) and second (y) samples from that locus (i). The Wilcoxon rank test was used to test whether the median of the distribution of the difference (F T i − V T i) over all loci was equal or larger than zero to determine whether it F T > V T for each pair of samples.Global tests were employed to evaluate significance of multiple tests. The sequential Bonferroni procedure (Holm, ) was applied to detect a single test‐specific departure (e.g., for a locus), whereas the binomial test (which estimates the probability of obtaining the observed number of significant tests at a given significance level given the total number of tests) can detect weaker departures across multiple tests, such as genomewide departures. Unless otherwise specified, calculations were carried out using programmes written by TL in SAS (). […]

Pipeline specifications

Software tools MEGA, Structure Harvester
Application Population genetic analysis
Organisms Anopheles coluzzii
Diseases Malaria
Chemicals Pyrethrins