1,2-diarachidonoyl-glycero-3-phosphocholine

We used STRUCTURE [1 (link),2 (link)] as a benchmark for the performance of DAPC. We analysed all simulated datasets with STRUCTURE v2.1, using the admixture model with correlated allele frequencies to determine the optimal number of genetic clusters and to assign individuals to groups. Computations were performed on the computer resources of the Computational Biology Service Unit at Cornell University (http://cbsuapps.tc.cornell.edu/). For each run, results were based on a Markov Chain Monte Carlo (MCMC) of 100,000 steps, of which the first 20,000 were discarded as burn-in. Analyses were ran with numbers of clusters (k) ranging from 1 to 8 for the island and hierarchical island models (Figure 2a-b), from 1 to 15 for the hierarchical stepping stone (Figure 2c), and from 1 to 30 for the stepping stone (Figure 2d). Ten runs were performed for each k value. We employed the approach of Evanno et al. [57 (link)] to assess the optimal number of clusters. In order to assess assignment success, STRUCTURE was run by enforcing k to its true value. Individuals were assigned to clusters using CLUMPP 1.1.2 [58 (link)], which allows to account for the variability in individual membership probabilities across the different runs. To obtain results comparable to DAPC, individuals were assigned to the cluster to which they had the highest probability to belong.

The methodological approach presented in the paper is implemented in the adegenet package [6 (link)] for the R software [27 ]. The function find.clusters runs successive K-means for a range of k values, and computes the BIC of the corresponding models. The basic K-means procedure is implemented by the function kmeans in the stats package [27 ]. DAPC is implemented as the function dapc, and relies on procedures from ade4 [55 ,59 ,60 ] and MASS [61 ] to perform PCA (dudi.pca) and DA (lda). Both find.clusters and dapc can be used with any quantitative data, and have specific implementations for genetic data. The analysis of the four simulated datasets presented in Figures 4 and 5 can be reproduced by executing the example of the dataset dapcIllus. Similarly, analyses of the extended HGDP-CEPH and of the seasonal influenza (H3N2) data can be reproduced by executing the example of the datasets eHGDP and H3N2, respectively. Documentation and support can be found at the adegenet website (http://adegenet.r-forge.r-project.org/).

Let X be a n × p genetic data matrix with n individuals in rows and p relative frequencies of alleles in columns. For example, in the case of a locus with three alleles (A₁, A₂, A₃), a homozygote genotype A₁/A₁ is coded as [1, 0, 0], while a heterozygote A₂/A₃ is coded as [0, 0.5, 0.5]. We denote X^j the j^th allele-column of X. Missing data are replaced with the mean frequency of the corresponding allele, which avoids adding artefactual between-group differentiation. Without loss of generality, we assume that each column of X is centred to mean zero. Classical (linear) discriminant analysis seeks linear combinations of alleles with the form:
$f(\mathbf{v})={\displaystyle \sum _{j=1}^p{\mathbf{X}}^j{v}_j}=\mathbf{X}v$
(v = [v₁...v_p]^T being a vector of p alleles loadings, known as 'discriminant coefficients'), showing as well as possible the separation between groups as measured by the F statistic (Equation 3). That is, the aim of DA is to choose v so that F(Xv) is maximum.
Linear combinations of alleles (Equation 5) optimizing this criterion are called principal components, which in the case of the discriminant analysis are also called discriminant functions. Discriminant functions are found by the eigenanalysis of the D-symmetric matrix [51 ]:
$\mathbf{P}\mathbf{X}{(\mathbf{W})}^{-1}{\mathbf{X}}^T{\mathbf{P}}^T\mathbf{D}$
where P is the previously defined projector onto the dummy vectors of H, and W is the matrix of covariances within groups, computed as:
$\mathbf{W}={\mathbf{X}}^T{(\mathbf{I}-\mathbf{P})}^T\mathbf{D}(\mathbf{I}-\mathbf{P})\mathbf{X}$
This solution requires W to be invertible, which is not the case when the number of alleles p is greater than the number of individuals n. Moreover, this inverse is numerically unstable ('ill-conditioned') whenever variables are correlated, which is always the case in allele frequencies and can be worsened by the presence of linkage disequilibrium.
To circumvent this issue, DAPC uses a data transformation based on PCA prior to DA. Rather than analyzing directly X, we first compute the principal components of PCA, XU, verifying:
${\mathbf{X}}^T\mathbf{D}XU=\mathbf{U}\Lambda$
where U is a p × r matrix of eigenvectors (in columns) of X^TDX, and Λ the diagonal matrix of corresponding non-null eigenvalues. Note that when the number of alleles (p) is larger than the number of individuals (n), we can alternatively proceed to the eigenanalysis of XX^TD to obtain U and Λ [55 ], which can save considerable computational time. By definition, the number of principal components (r) cannot exceed the number of individuals or alleles (r ≤ min(n, p)), which solves the issue relating to the number of variables used in DA. Moreover, principal components are, by construction, uncorrelated, which solves the other issue pertaining to the presence of collinearity among allele frequencies.
DA is then performed on the matrix of principal components. At this step, less-informative principal components may be discarded, although this is not mandatory. Replacing X with XU into Equation 6, the solution of DAPC is given by the eigenanalysis of the D-symmetric matrix:
$\mathbf{P}\mathbf{X}U{({\mathbf{U}}^T\mathbf{W}U)}^{-1}{\mathbf{U}}^T{\mathbf{X}}^T{\mathbf{P}}^T\mathbf{D}$
The first obtained eigenvector v maximizes b(XUv) under the constraint that w(XUv) = 1, which amounts to maximizing the F-statistic of XUv. This maximum is attained for the eigenvalue γ associated to v (i.e., F(XUv) = γ). In other words, the loadings stored in the vector v can be used to compute the linear combinations of principal components of PCA (XU) which best discriminate the populations in the sense of the F-statistic.
However, it can be noticed that these linear combinations of principal components ((XU)v) can also be interpreted as linear combinations of alleles (X(Uv)), in which the allele loadings are the entries of the vector Uv. This has the advantage of allowing one to quantify the contribution of a given allele to a particular structure. Denoting z_j the loading of the j^th allele (j = 1,...,p) for the discriminant function XUv, the contribution of this allele can be computed as:
$\frac{z_j^2}{{\displaystyle \sum _{j=1}^p{z}_j^2}}$

We conducted a principal components analysis (PCA) and a discriminant analysis of principal components (DAPC) using the 301 K Set in ADEGENET [31 (link)] (functions: glPca, dapc) to investigate genetic differentiation between the two populations. We retained 30 principal components (PCs) following cross-validation procedures in ADEGENET (xvalDapc). Clustering procedures for the initial DAPC highlighted 12 individuals that grouped with the opposite capture location (function: find.clusters). We removed these individuals before further analysis and used SNP data for the remaining 104 dogs for calculating F_ST, measuring heterozygosity, and identifying outlier loci.

We used LOSITAN to identify outlier loci associated with directional selection [46 (link)]. Outlier analyses such as these frequently produce some degree of false positives [47 (link)]. To overcome this, we selected the 84 K set with a higher MAF threshold for this analysis to highlight the strongest signals of selection and biasing towards genomic regions with large effects. We included only the 104 dogs that grouped with their capture location in the DAPC for better identification of the differentiating loci in each population. We ran 1,000,000 simulations and applied a 95% confidence interval and a false discovery rate of 0.1. We considered loci significant when designated as ‘candidate positive selection’ and when the calculated F_ST was higher than the expected F_ST in all simulations (P(Simulated F_ST < sample F_ST) = 1).
The genome location of each significant SNP identified using LOSITAN was expanded to a 10 kb genomic interval of CanFam3 (5 kb either side of the SNP coordinates). We then surveyed each of the 10 kb regions for the presence of genes and identified corresponding gene ontology (GO) terms through the Mouse Genome Informatics Batch Query Search [48 (link)]. GO terms were evaluated for putative associations towards the exposures faced within the environment (e.g., GO:0,010,212 “response to ionizing radiation”).

The provenance assignment success of the DArTcap markers was tested with asssignPOP v. 1.2.2 (Chen, Marschall, et al., 2018 (link)) and rubias v.0.3.2 (Anderson et al., 2008 (link); Moran & Anderson, 2019 (link)). Assignment accuracy was tested with assignPOP, using both the Monte‐Carlo and K‐fold cross‐validation procedures to test the assignment of a hold‐out data set with 1000 iterations. We tested power of the markers by selecting a subset of loci with the highest F_ST values (5%, 10%, 50%, and 100% of all loci) to train the assignment model. Similarly, the assignment accuracy of simulated mixed groups, based on a reference leave‐one‐out dataset, was evaluated with rubias (Anderson et al., 2008 (link)). Known simulated proportions for each reporting unit were compared with the numbers estimated by rubias. Populations with a sample size of one (i.e., Sierra Leone, eastern Atlantic) were excluded from these analyses. We also examined the minimum number of informative markers needed to assign provenance by subsampling 5–500 markers based on loading contributions of each principal component from the DAPC analysis and testing the assignment accuracy with rubias.