We analyzed simulated data published by Lotterhos and Whitlock (2015). These data consist of 10,000 haploid, biallelic loci (representing SNPs: 9,990 neutral, 100 adaptive) on three replicate, quasi-continuous, square landscapes with 360 × 360 grid cells, each housing a deme. The simulated species has a large geographic range, high effective population size, and rapid linkage decay, with independently generated loci (Lotterhos and Whitlock, 2014). See the Supporting Information, Appendix S1 of Lotterhos and Whitlock (2014) for details of the simulator.

Variation in population demographic history was simulated under four models (Figure 1): island model (IM), isolation-by-distance (IBD), expansion from one refugium (1R) or two refugia (2R). The IM model served as a reference: its genetic structure is non-spatial, as migrants are drawn from a migrant pool independent of the distance between demes. The IBD model at equilibrium represents a stationary spatial process, whereas the 1R and 2R models are spatial and non-stationary, as allele frequencies are likely to show trend along clines. The carrying capacity per deme equaled 16, 71, and 124 for IBD, 1R and 2R, respectively. For all models except the IM, dispersal was modeled with a discretized Gaussian dispersal kernel with σ = 1.3 multiplied by cell width. The simulation controlled for global Fst by sampling at the time period when the global Fst was approximately equal to 0.05, which was reached after 1,000 generations for the 1R and 2R, 5,000 for the IM and 10,000 for the IBD model. This means that datasets are not directly comparable between demographic histories (Lotterhos and Whitlock, 2014, 2015).

Lotterhos and Whitlock (2015) sampled two sets of 6 and 20 individuals per deme, using three different sampling designs (random, transect, and pairs) with n = 90 sampling locations each. Here we used only the random samples of 90 sites. Note that the sampling locations were constant across demographic scenarios and replicate landscapes, and the spatial coordinates are available in Wagner et al. (2017).

For each combination of four demographic scenarios and three replicate landscapes (12 datasets), we determined deme-level allele frequencies among the 20 sampled individuals for each of the sampled 90 demes, using 9,900 neutral loci (full data). To assess the effects of genetic resolution, we repeated analyses with allele frequencies based on the sampling of 6 individuals per deme, and with a random sample of 3,300 or 500 loci. To assess the effects of the spatial sampling design, we repeated analyses for 10 random subsamples containing 60 or 30 of the original 90 sampling locations, using the same subsamples across datasets.

For each dataset, we used the R (R Core Team, 2020) package “hierfstat” (Goudet and Jombart, 2015) to estimate sample Fst with the function “basic.stats” and calculated two pairwise genetic distance matrices (Fst: pair-wise Fst, Dch: Cavalli-Sforza and Edwards Chord distance) with the function “genet.dist” of the R package “hierfstat.”

MEM spatial eigenvectors were derived from Euclidean distances between sampling locations with the function “mgMEM” of the R package “memgene” (Galpern et al., 2014). We used default settings, which implement the minimum spanning tree (MST) truncation, where distances exceeding the largest distance in the MST (dMST) are replaced by 4 × dMST.

We used the function “mgForward” of the package “memgene” to select those MEM spatial eigenvectors with positive eigenvalues (i.e., representing positive spatial autocorrelation) that were significantly associated with the genetic distance matrix. We applied default settings with 100 permutations and alpha = 0.05, which is sufficient for a one-sided test.

We considered two main response variables, the MEMgene adjusted R² and Moran’s I. While Fst quantifies the overall genetic structure, spatial or not, the MEMgene adjusted R² quantifies how much of the overall genetic variation is spatial (Peres-Neto and Legendre, 2010; Galpern et al., 2014), and Moran’s I quantifies the degree of positive spatial autocorrelation.

The MEMgene adjusted R², as returned by the function “mgForward,” is the correlation coefficient of a distance-based redundancy analysis (dbRDA) of the genetic distance matrix regressed against the set of significant MEM spatial eigenvectors with positive spatial autocorrelation (see above), adjusted for the number of predictors (Galpern et al., 2014).

A set of n = 90 spatial locations will result in 89 MEM orthogonal spatial eigenvectors (Dray et al., 2006). Any set of n− 1 variables will be able to explain 100% of the variation in any response variable observed at the n locations. Because the MEM spatial eigenvectors are orthogonal, the correlation of each one of them with the genetic data, and thus the R² explained by each MEM spatial eigenvector, does not depend on the other eigenvectors included in the model. The total genetic variation can thus be partitioned by the MEM spatial eigenvectors, resulting in a scalogram S, a vector with n− 1 eigenvector-specific R² values that sum to 1 (Dray et al., 2012). Each MEM spatial eigenvector k represents a synthetic spatial pattern, and the MEM eigenvalue associated with eigenvector k, rescaled through dividing by a constant, is equal to Moran’s I of this pattern (Dray, 2011). Note that MEM spatial eigenvectors are sorted by their eigenvalue, from largest to smallest, thus representing a gradient from large-scale patterns, with positive spatial autocorrelation, to finest-scale patterns, with negative spatial autocorrelation. We calculated Moran’s I of the genetic data as a weighted mean of all n-1 rescaled MEM eigenvalues, weighted by the scalogram S (Dray, 2011).

In the absence of spatial autocorrelation, Moran’s I has an expected value of E(I) = −1/(n – 1). Moran’s I values vary roughly between + 1 and −1, however, the exact bounds will depend on the sampling design. Specifically, the maximum value is given by the Moran’s I of the first MEM spatial eigenvector (k = 1), and the minimum by the Moran’s I of the last MEM spatial eigenvector (k = n – 1) (Dray, 2011). Thus, Moran’s I is largest if 100% of the genetic variation is explained by the first MEM spatial eigenvector. Because we used the same n = 90 sampling locations for all datasets, the maximum value of Moran’s I was constant across the full datasets. However, each replicate subsample with 60 or 30 locations will have its own maximum value. To account for this effect, we also calculated a rescaled Moran’s I (Ir) as Ir = I/max(I).

Standardized effect sizes (SES) are calculated by subtracting the mean effect size (ES_sim)of R simulated data sets from the observed effect size (ES_obs) and dividing by the standard deviation (σ_sim) of the simulated values (Gotelli and McCabe, 2002):

We simulated R = 200 datasets for each full dataset with 90 sampling locations, 20 individuals sampled per location, and 500 loci, and for each subset (with 30 or 60 sampling locations) thereof. A higher number of replicate simulations (e.g., 1,000) is generally recommended but was not feasible here due to computational constraints. The relatively low number of simulations will increase variability but is not expected to create bias in results. For each simulated dataset, genotypes of the sampled individuals were randomly permuted among the sampled locations. Then, we proceeded with the analysis as for the observed data (here, the data simulated by Lotterhos and Whitlock, 2015) to obtain an estimate of each response variable for each simulation run. The means and standard deviations among the 200 replicates were then used according to Eq. (1) to calculate SES.

We visually checked for deviations from normality of the distribution of ES_sim for the full data with 90 sampling locations under the IBD scenario. Normal probability plots for Moran’s I and rescaled Moran’s I showed no systematic deviations from normality, whereas the MEMgene adjusted R² showed a distinctly non-normal distribution with most values being zero and a few deviations in either direction. Note that the distribution did not improve when using unadjusted R² values. As deviation from normality may invalidate the use of SES (Botta-Dukát, 2018), we decided to report SES only for Moran’s I and rescaled Moran’s I.

To approximate a meta-analysis situation with independent analysis of multiple datasets, we performed MEMgene separately for each subsample with 60 or 30 sampling locations (meta-analysis mode). For n = 60 locations, any pair of subsamples will share at least 30 demes (50% of sampling locations), but their spatial configuration and thus MEM spatial eigenvectors will vary. While real meta-analyses will involve more independent datasets, this setting ensures that the subsamples are as comparable as possible and can be expected to exhibit similar spatial genetic structure as the full dataset with n = 90 sites, at least on average.

In a multi-species study where all species are sampled at the same locations, some species may be absent from some locations. As an alternative to performing MEMgene independently for each species (meta-analysis mode), we considered using the same, full set of 89 MEM spatial eigenvectors and dropping, for each species, the sites with missing values (comparative mode). This means that for each subsample, forward selection of MEM spatial filters with dbRDA was performed with the same predictor variables but including only the 60 observations with genetic data available. Note that because we used forward selection and included only those spatial eigenvectors with positive eigenvalues, this did not lead to overfitted models with more predictors than sites, and the adjusted R² could be calculated as above. For Moran’s I, we determined the scalogram S with a set of n – 1 simple regressions, one for each spatial eigenvector. Dropping sites is likely to render the spatial eigenvectors non-orthogonal, so that the sum of S may differ from 1. We chose not to normalize S (i.e., make it sum to 1) as preliminary results found that this led to decreased performance.

For each dataset and response variable, we compared the value obtained for the full sample of n = 90 sampling locations to the mean value of the 10 subsamples with n = 60, which we calculated separately for the meta-analysis mode and for the comparative mode. While we expected some variability among the 10 replicate subsamples, a significant deviation of their mean from the full sample will indicate systematic bias. We used paired t-tests to statically validate such bias across datasets and evaluated effect size using Cohen’s d (Cohen, 1988).

We observed considerable differences in the amount of spatial genetic structure identified for each demographic scenario by Fst, MEMgene adjusted R² (Figure 1), and Moran’s I (Figure 2). In all datasets, the demographic model where the greatest amount of genetic structure is explained by space was the 2R model, followed by 1R, and IBD. As expected, the IM showed Fst values similar to the other demographic scenarios but no spatial structure, and thus was dropped from further analysis. The difference in the strength of spatial genetic structure between the IBD, 1R, and 2R models was greater for Moran’s I than for adjusted R². This overall pattern was robust to the metric of genetic distance used, as we observed similar patterns to those based on pairwise Fst upon repeating the analysis using Dch distance (see Supplementary Figure 1).

Both adjusted R² and Moran’s I showed systematically lower values when calculated with six individuals per location, compared to 20 per location (Figure 3). In contrast, the number of loci had no discernible impact on the mean value of either response variable, though smaller numbers of loci increased variability somewhat (Figure 3). Fst was not sensitive to the number of loci or the number of individuals genotyped.

Across all metrics, datasets with all 90 sites had higher values than those with 60, with the lowest values for those with 30 sites (Figure 4). For the IBD demographic model, meta-analysis mode (where the MEM spatial eigenvectors are derived independently for each species based on those sites where the species occurred) increased this bias compared to comparative mode (where the same MEM spatial eigenvectors are used for all species) for adjusted R² and for Moran’s I, but the opposite trend was observed for scenarios with expansion from one or two refugia. Rescaling Moran’s I by its maximum value did not eliminate bias due to the number of sampling locations (Figure 4, bottom). When comparing the mean value of the then subsamples with 30 sites to the corresponding value for the full sample of 90 sites, paired t-tests showed large, statistically significant differences, with the largest effect size (Cohen’s d) for Moran’s I, the smallest for rescaled Moran’s I, and an intermediate effect size for the MEM adjusted R² (Table 1). In meta-analysis mode, standardized effect sizes (SES), reported for Moran’s I and rescaled Moran’s I (Figure 4), were equally biased as unstandardized values when comparing effect sizes between datasets with 90, 60, or 30 sampling locations. In comparative mode, SES performed considerably worse.

We found that expansion from one or two refugia resulted in much larger effect sizes for the strength of spatial genetic structure (i.e., MEM adjusted R² and Moran’s I values), although the datasets had been simulated with comparable overall genetic differentiation (Fst). These results suggest that we cannot accurately compare the effect of landscape features on gene flow between independent studies or within a multispecies study where there is variation in underlaying population demographic history. The underlaying population demography is typically unknown to the empirical researcher and may vary between landscapes for the same species and between species for the same landscape, which further increases the uncertainty of the potential comparison between studies or species.

The strong increase of adjusted R² and Moran’s I for scenarios with expansion from one or two refugia suggests that effect size in MEMgene is sensitive to non-stationarity. It is well known that several MEM spatial filters are required to model linear trend, hence one strategy is to detrend data (Borcard et al., 2004, 2018; Griffith and Peres-Neto, 2006; Beale et al., 2010). As the response is a matrix of genetic distances, detrending is not trivial. One approach would be to detrend the allele frequencies. While relative frequencies vary between 0 and 1, removing the mean (modeled as a function of spatial x- and y-coordinates) will result in some negative values, which makes it problematic to calculate genetic distances. On the other hand, one could partial out the spatial x- and y-coordinates during spatial filtering. However, this would likely render the spatial filters non-orthogonal. When we partial out the spatial coordinates from the correlation between the genetic distances and MEM spatial eigenvectors, we regress both on the coordinates and correlate their residuals (Legendre and Legendre, 2012). The residuals of the spatial eigenvectors are no longer uncorrelated. For the sampling design in Figure 1, for instance, the pairwise correlations among the residual eigenvectors range between −0.4 < r < 0.4.

Furthermore, trend in allele frequencies, such as along a cline, may represent biologically meaningful population genetic structure, so that its elimination could also distort comparisons between studies. Further research is needed to assess whether expansion from refugia creates other forms of non-stationarity beyond trend in allele frequencies. It is also important to recognize that this study only examined two metrics based on MEM and that other methods with different assumptions may be more robust between comparisons.

We found that both adjusted R² and Moran’s I were surprisingly robust to the number of loci, but sensitive to the number of individuals sampled per site. This is in contrast to Fst, which was robust to both changes in the number of samples per location and the number of loci. Our findings contrast with those of Landguth et al. (2012), who found that the partial Mantel r correlations between genetic and ecological distances (accounting for geographic distances) was sensitive to the number of loci and the level of polymorphism but robust to the number of individuals (with one individual sampled per location). This suggests that existing recommendations on resource allocation (i.e., when to invest into sampling more sites, genotyping more individuals per site, or increasing the number or polymorphism of genetic markers) for landscape genetic studies (Landguth et al., 2012; Oyler-McCance et al., 2013) may need to be revisited to clarify which factors are most important depending on the research question and the chosen analysis approach (Balkenhol and Fortin, 2015).

This study only examined neutral SNPs, and expanding into other common genetic markers (Hall and Beissinger, 2014), such as microsatellites, or markers under selection is beyond the scope of this study. However, in empirical datasets, other factors likely to occur (e.g., site polymorphism, strength of selection, linkage to sites under selection) may have an effect on genetic resolution. Our results suggest that the statistical power of MEMgene analysis may be more dependent on the number of individuals sampled per deme than on the number of loci. Further research should address whether sampling more individuals per location, or pooling nearby samples, increases the ability of MEMgene to detect cryptic spatial structure.

We found that smaller samples (i.e., subsamples with 60 or 30 of the original 90 sampling locations) resulted in systematically lower estimates of spatial genetic structure, both for the MEMgene adjusted R² and for Moran’s I. This sensitivity to the number of sampling locations suggests difficulty in comparing between landscape genetic studies in meta-analyses. Moreover, such bias in effect size may not be limited to the analysis of genetic variation with MEMgene but could affect other applications of MEM to univariate or multivariate ecological data (Dray et al., 2012). More generally, Moran’s I has been shown to underestimate spatial autocorrelation for small numbers of spatial locations (Carrijo and da Silva, 2017).

Analysis with comparative mode, where the same basis set of MEM spatial eigenvectors (based on all sites) are used for all species, somewhat reduced bias under the stationary scenario of isolation-by-distance but increased it under the non-stationary scenarios with expansion from one or two refugia. Where there is a small number of missing sites, it may be safer to restrict analysis to those sites where all species occurred. We do not recommend using comparative mode analysis, as the bias reduction was insufficient and dependent on stationarity. In addition, there are methodological concerns where the missing observations change the correlation structure of the MEM spatial eigenvectors resulting in scalograms that do not sum to 1.

Standardized effect sizes (SES) did not provide an improvement. First, non-normal distributions among permuted datasets precluded the calculation of SES for the MEMgene adjusted R². Second, SES for subsamples in meta-analysis mode, on average, showed the same bias compared to SES for the full dataset with 90 sampling locations as observed Moran’s I (or rescaled Moran’s I). Finally, comparative mode analysis further decreased the performance of SES.

This study only examined a random sampling design with randomized missing values. While we suspect similar problems will occur with other common sampling designs, such as transects or grids (Oyler-McCance et al., 2013), those analyses are beyond the scope of this paper.