Sylvain S Piry
email: piry@supagro.inra.fr

Marie-Pierre Pierre Chapuis

Bertrand Gauffre

Julien J Papaix

Astrid Cruaud

Karine Berthier

Julien Papaïx

Mapping Averaged Pairwise Information (MAPI): a new exploratory tool to uncover spatial structure

Keywords: exploring spatial structure pairwise information, geovisualisation, exploratory method, spatial genetics, landscape genetics, isolation-by-distance, barrier, gene flow, GIS, database

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Computation of metrics to estimate similarity, difference or flow between samples is a common approach to study a wide range of social and biological systems [START_REF] Borgatti | Network analysis in the social science[END_REF]Miele et al. 2015). Computing pairwise metrics allows the construction of networks with nodes (samples) being connected by edges that are associated with the value of the metric of interest. When geographic coordinates of samples are available, such networks are inherently informative about the spatial organization of the system under study [START_REF] Barthélemy | Spatial Networks[END_REF]. This organization can be visualised by superimposing the network on a geographic map, with the thickness or colour of the edges displayed according to the pairwise values. However, such graphical representations are difficult to interpret for large networks as the numerous nodes and edges strongly overlap and obscure both the network and the underlying geographic space (Hennemann 2013). To better visualise spatial information, variation surfaces can be produced using interpolation procedures (e.g. kernel density estimation, kriging). To do so, pairwise metric values have first to be attributed to unique geographical points (i.e. transformed into punctual values) as, for example, to the middle of the edges connecting the samples in the network (e.g. Miller 2005).

Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. ( 2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology and Evolution, 7 (12), 1463-1475. , DOI : 10.1111/2041-210X.12616

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

In this paper we present a novel method, MAPI (Mapping Averaged Pairwise Information), to produce variation surfaces from pairwise metrics computed between georeferenced samples. In essence, MAPI is a non-parametric smoothing procedure applied to pairwise values rather than punctual values. To smooth out pairwise values, the straight lines connecting the samples in the network are replaced by ellipsoidal polygons for which the foci are located on the two samples being connected. Then, similarly to kernel density estimators for punctual values, ellipses constitute a geometric shape used to average information between overlapping connections. This smoothing procedure produces a two-dimensional geographical layer that can be easily visualised and customised in a GIS. Although this method may be of interest in other research fields, we will hereafter focus on spatial and landscape genetics, where pairwise metrics are widely used but flexible visualisation tools are still lacking for these measures (see Miller 2005;Vandergast et al. 2010, Etherington 2011, Petkova et al. 2016).

Although new methods, such as spatially-explicit clustering, are increasingly used in spatial genetics to describe genetic variation and identify barriers to gene flow (Guillot et al. 2009), pairwise genetic measures, such as distance, remain appealing as: i) they can be computed between individuals or populations, ii) they facilitate handling massive datasets, which are increasing with the advance in high-throughput sequencing (Duforet-Frebourg and Blum 2014) and, iii) their regression against geographic distances may be informative with respect to isolation-by-distance (IBD), i.e. the increase of genetic differentiation with spatial separation due to restricted dispersal (Wright 1943).

The original IBD model has been further extended in landscape genetics to assess the impact of landscape heterogeneity on spatial genetic structure and gene flow. In this field, replacing straight line distances between samples by ecological distances has given rise to concepts of isolation-bybarrier, isolation-by-resistance and isolation-by-environment [START_REF] Balkenhol | Landscape Genetics: Concepts, Methods, Applications[END_REF]. Nowadays, various approaches using least cost paths, circuit theory or population graphs have been developed to assess the impact of landscape features on genetic structure and gene flow (e.g. Coulon et al. This article is protected by copyright. All rights reserved.

2004; [START_REF] Cushman | Gene Flow in Complex Landscapes: Testing Multiple Hypotheses with Causal Modeling[END_REF]McRae 2006;McRae & Beier 2007;Garroway et al. 2008;[START_REF] Dyer | Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure[END_REF][START_REF] Dyer | Population Graphs and Landscape Genetics[END_REF], Petkova et al. 2016).

Landscape genetics approaches often require a priori knowledge on species-landscape interactions to target relevant environmental variables, set cost/resistance values and account for demographic effects such as population size (e.g. [START_REF] Broquet | Genetic isolation by distance and landscape connectivity in the American marten (Martes americana)[END_REF]Weckworth et al. 2013). These methods can be difficult to apply when little is known about the species under study. In such a situation, we still rely on exploratory analyses to identify environmental variables of interest and draw up hypotheses that can be further tested using appropriate sampling schemes (Kelling et al. 2009;Richardson et al. 2016). Although the visualisation of pairwise genetic measures may be an obvious first step in such an approach (Etherington 2011), there are still very few tools allowing visualisation of such measures along with environmental layers without, first, attributing features such as cost or resistance values to habitat types. Among those tools, Allele In Space (Miller 2005) and the GIS toolbox of Vandergast et al. (2010) offer the rare possibility to produce variation surfaces of genetic distances that can be mapped over landscape layers (Wood et al. 2013;[START_REF] Adams | Influence of ecological and geological features on rangewide patterns of genetic structure in a widespread passerine[END_REF]. Both tools rely on an inverse-distance-weighted interpolation of pairwise measures that are first spatially attributed to the middle of the segments linking the samples in the network.

Confounding effects due to IBD can be somewhat limited by i) using residuals from the regression of genetic distances against geographic distances (Miller 2005) or, ii) limiting the network to the nearest neighbours for each sample (Vandergast et al. 2010).

Disentangling the relative contributions of IBD and environmental features in shaping spatial genetic patterns remains a central issue in landscape genetics [START_REF] Bradburd | Disentangling the effects of geographic and ecological isolation on genetic differentiation[END_REF]Wang et al. 2014). This question is often addressed using Mantel and partial Mantel tests between matrices of population-or individual-based genetic distances, Euclidean distances and ecological distances (e.g.

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved. [START_REF] Cushman | Gene Flow in Complex Landscapes: Testing Multiple Hypotheses with Causal Modeling[END_REF]Hagerty et al. 2011). There is an ongoing debate on the statistical performance of the Mantel test in spatial and landscape genetics (see Raufaste & Rousset 2001;Legendre & Fortin 2010;[START_REF] Cushman | Spurious correlations and inference in landscape genetics[END_REF]Guillot & Rousset 2013;[START_REF] Cushman | Re-evaluating causal modeling with Mantel tests in landscape genetics[END_REF]) and, recently, alternative methods using multivariate or geostatistical techniques have been proposed to disentangle the effects of geographic distance and environmental heterogeneity on spatial genetic structure (e.g. [START_REF] Duforet-Frebourg | Nonstationary patterns of isolation-by-distance: Inferring measures of local genetic differentiation with bayesian kriging[END_REF][START_REF] Bradburd | Disentangling the effects of geographic and ecological isolation on genetic differentiation[END_REF]Galpern et al. 2014;[START_REF] Botta | Enhanced computational methods for quantifying the effect of geographic and environmental isolation on genetic differentiation[END_REF]. Moreover, results from lot of methods do not include spatially-explicit variation surfaces that can be mapped over landscape layers to facilitate their interpretation.

In this general context, MAPI provides an approach to i) visualise pairwise genetic measures without confounding effects resulting from IBD, ii) test for spatial genetic discontinuity through a nonparametric randomisation procedure and, iii) explore relationships between observed genetic patterns and environmental heterogeneity, notably by using MAPI results for further statistical analyses. In this work, we assessed the efficiency of MAPI to accurately detect spatial discontinuity in pairwise genetic metrics by applying permutation tests on controlled simulated datasets of panmictic populations and populations under isolation-by-distance, including scenarios of separation by a linear barrier to gene flow. We also appraised the potential of MAPI in a landscape genetics framework by analysing genotypes simulated under landscape constraints with various spatial configurations of favourable and unfavourable habitats. We used a Bayesian conditional autoregressive model to analyse the relationships between MAPI results and landscape variables. To illustrate MAPI we also re-analysed published data from three biological models: 1) microsatellite genotypes from a rodent population under IBD, 2) DNA sequences from a plant virus exhibiting spatial genetic discontinuities and, 3) microsatellites and landscape data from both a forest specialist and a generalist species of ground beetles. Finally, we used controlled and observed datasets to explore MAPI sensitivity to parameter setting and sampling scheme.

Grid of cells.

A grid of hexagonal cells is superimposed on the area defined by the convex hull of the sampling points. Based on the Nyquist Frequency concept, optimized cell size should be at most half the average distance between closest samples. The half-width (h w ) of the cells is computed as:

ℎ = ² √ .
with N the number of sampling points, A the surface's area defined by their convex hull and a parameter depending on their spatial dispersion: 0.5 for regular sampling and 0.25 for irregular sampling (see Hengl 2006). The grid is automatically generated by informing the parameter. The final number of cells (n c ) in the grid is then computed as: ≈ . . . Alternatively, a user-defined grid of cells can be imported in PostgreSQL (see user manual for details on how to build or import a grid).

This article is protected by copyright. All rights reserved.

Network of ellipses. Network edges are formed by ellipsoidal polygons (hereafter referred to as ellipses) for which foci are the geographical locations of the two samples being connected. The shape of these polygons can be adjusted by two parameters: 1) the eccentricity of the ellipses which controls the smoothing intensity and must be greater than 0 (infinite circle) and smaller than 1 (straight line) and 2) the radius of the error circle that controls for uncertainty on sample coordinates (error_circle_radius). In addition, two optional parameters (min_distance and max_distance) limit the analysis to a given range of between-sample distances. Effects of these four parameters on the shape of the network are detailed in section 1 of the Appendix.

Surface. The ellipses receive the value of the metric computed between the samples they connect and the cells of the grid receive the weighted arithmetic mean (m w ) of the n e ellipses intercepting their geographical extent, computed as: = ∑ , with v i and a i the metric value and area of the ellipse i, respectively, and s w , the sum-of-weights of the ellipses defined as: = ∑ . This weighting procedure limits long distance effects as long and inflated ellipses participate far less than short ellipses to the computation of m w . Cells not intersected by ellipse have no m w value and are not included in the final result.

Test of significance for spatial structure. We implemented a nonparametric randomisation procedure, as schematised in the Figure 2, to test whether the pairwise metric values associated to the ellipses are independent of the sample locations (i.e. under the null hypothesis that pairwise metric values, and then m w cell values, are randomly distributed in space). Sample locations from the observed dataset (Figure 2A1) are permuted n p times (Figure 2B1) (parameter n_permutations). At each permutation, new m w values are computed and stored to build a cumulative null distribution for each of the n c cells of the grid (Figure 2B2). Each m w cell value from the observed dataset is then ranked against its null distribution (Figure 2C1). For each cell, the proportion of permuted values that are smaller than the observed value provides a lower-tailed test (LT) p-value. An upper-tailed

Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. ( 2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology andEvolution, 7 (12), 1463-1475

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

test (UT) p-value is also computed for each cell as: 1 -(LT p-value). As the probability to find significant cells only by chance (i.e. type I error) increases with the number of tests performed (i.e. n c tests), the False Discovery Rate (FDR) procedure proposed by [START_REF] Benjamini | The control of the False Discovery Rate in multiple testing under dependency[END_REF] is applied to account for multiple testing under positive dependency conditions (i.e. spatial autocorrelation between cells) (Figure 2C2). The significance level at which FDR is controlled can be set by users through the parameter my_alpha. For example, when my_alpha is set to 0.05, this means that 5% of the cells detected as significant can be false positives. Finally, for each test, significant cells that are spatially connected are aggregated together. When using metrics estimating differences, such as distances, significant highest m w values localise areas of higher dissimilarity than expected by chance (hereafter referred to as discontinuous areas) while significant lowest m w values localise areas of higher similarity than expected by chance (hereafter referred to as continuous areas).

Data output.

Cell-specific information, such as their geometry and associated m w and p-values, are stored as a spatial PostgreSQL table which can be exported as a shapefile or text file and imported into GIS software or the software R using the rgdal package. These data can be used for further statistical analyses of relationships between genetic measures and landscape variables.

Visualisation. On the final graphical output (Figure 2D), the observed cell values (m w ) are displayed using a colour scale. Cells with significant extreme values after FDR correction are visualised as black contours mapped as an additional layer. When significant cells are spatially connected, only the perimeter of the polygon they form is represented. Vector or raster maps, satellite imagery or teledetection products can be inserted below the MAPI layers for landscape interpretation.

This article is protected by copyright. All rights reserved. Analyses of simulated test datasets with MAPI were performed using a grid of hexagonal cells defined by setting the parameter to 0.5 and 0.25 for regular and irregular samplings, respectively (Table 1). The network was built using the genetic distance a r (Rousset 2000), an error radius for sample location of 0.01 and ellipses with an eccentricity value of 0.975. Cells with extreme high m w values were detected using 1,000 permutations of the sample locations and a significance level of 0.05 for computing significance for each cell. As the exact location of the barrier was known, we computed the proportion of the barrier covered by cells with significantly higher values

Simulated data in spatial genetics

This article is protected by copyright. All rights reserved.

(hereafter referred to as Coverage). We also computed the proportion of the study area covered by cells with significantly higher values having no contact with the simulated barrier, which may result from stochasticity, edge effects or IBD (hereafter referred to as Unexpected). We also investigated the sensitivity of MAPI to the eccentricity value (0.900, 0.975 and 0.999) and number of cells (setting the parameter to 0.15, 0.25, 0.5 and 0.75). To this aim, we used a subset of datasets simulated under both large regular sampling and small irregular sampling (see section 2.1 in the Appendix to identify test datasets).

Simulated data in landscape genetics

Using a forward-in-time (individual-based) simulation algorithm, we simulated datasets of 10 microsatellite genotypes for diploid individuals distributed in a landscape raster of about 50x50 cells.

Three landscape models, which consider two habitats with contrasted carrying capacity (20 and 2 for the favourable and unfavourable habitat, respectively) were analysed (20 simulated datasets each):

(1) a spatial transition from the favourable to unfavourable habitat with a high level of interpenetration (landscape test datasets 1; Gradient); (2) fragmentation of the favourable habitat in small areas isolated by a prevalent unfavourable habitat (landscape test datasets 2; Fragmentation);

(3) a random distribution of the favourable and unfavourable habitats over the study area (landscape test datasets 3; Random) (see Figure 4 to visualise the simulated landscapes). For the latter, the detection of the habitat effect is unexpected since the spatial scale of autocorrelation in habitat is very small. Details on simulation algorithms, parameter settings and basic measures of genetic variation can be found in sections 2. This article is protected by copyright. All rights reserved.

number of individuals was randomly sampled from each square of a 3x3 grid encompassing the landscape raster.

Analyses of simulated datasets with MAPI were performed using grids constituted of squared cells that matched the landscape raster pixels (side length=1), the genetic distance a r (Rousset 2000), an error radius for sample location of 0.5 to consider that individuals can be anywhere within a cell, ellipses with an eccentricity of 0.975 and a minimal distance between samples of 0.1 to exclude intra-cell connections as individuals simulated within a same cell have the same geographical coordinates. For each simulation and sampling scheme, we extracted from each cell both, the m w value and habitat type, which was expressed as a factor according to favourability (hab): 1 for the favourable habitat and 2 for the unfavourable habitat. We fitted a regression model

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

Biological datasets

To illustrate MAPI behaviour in the presence of IBD, we analysed microsatellite genotypes from a rodent population that was previously characterised as a single genetic unit, only structured by IBD (Gauffre et al. 2008), using the clustering method GENELAND (Guillot et al. 2005). MAPI efficiency in detecting genetic discontinuities is further illustrated using DNA sequences from different strains of plant virus. These data were previously analysed [START_REF] Desbiez | Emergence of new strains of Watermelon mosaic virus in South-eastern France: Evidence for limited spread but rapid local population shift[END_REF]Joannon et al. 2010) using phylogenetic analyses, the clustering method SAMOVA [START_REF] Dupanloup | A simulated annealing approach to define the genetic structure of populations[END_REF]) and the maximum difference Monmonier's algorithm (Monmonier 1973). Finally, MAPI was applied to genetic distances computed from microsatellite genotypes of two species of forest ground beetles with contrasted level of habitat specialisation (i.e. specialist and generalist). Previous work using Mantel tests between genetic, geographic and landscape distances showed that open field areas were a stronger barrier to gene flow for the forest specialist [START_REF] Brouat | Fine-scale genetic structure of two carabid species with contrasted levels of habitat specialization[END_REF]. See section 3 in the Appendix for further details on the datasets and related works.

RESULTS

All simulation setups and MAPI graphical outputs are available online at: https://www1.montpellier.inra.fr/CBGP/software/MAPI/.

Sensitivity to isolation-by-distance

Simulation test datasets. MAPI did not detect unexpected areas of genetic discontinuity under strict IBD, regardless of its strength (see test datasets 1 in Table 2; Figures 3E,3I). When models combine IBD and barrier effects, MAPI still did not detect unexpected areas of genetic discontinuity, i.e. which did not overlap with the simulated barrier (see test datasets 3 in Table 2; Figures 3F-H

, 3J-L). Under

This article is protected by copyright. All rights reserved. large regular sampling, the average percentage over replicates of unexpected discontinuous areas was less than 0.15% whatever the simulation setting considered (see 'Unexpected' in Table 2), with a maximal percentage across the 625 simulated datasets of 1.3%. Under small irregular sampling, the percentage of unexpected discontinuous areas was still minute, with an average lower than 0.53% whatever the simulation setting (Table 2) and a maximum across the 625 simulated datasets of 5.8%. Thus, gaps between samples did not drastically increase the proportion of unexpected significant discontinuous areas, even under strong IBD (g=0.25, slope = 0.055). The single exception was the combination of strong IBD and genotype datasets of 100 loci, for which the small and irregular sampling led to a more significant (but still relatively low) increase in the percentage of unexpected discontinuous areas, with an average over replicates of 2.75% and a maximum of 10.8% (Table 2).

However, the spatial distribution of the unexpected significant cells hardly suggested the presence of a barrier to gene flow (see online illustrations).

Biological dataset.

When applied on rodent microsatellite dataset, MAPI did not find any spurious significant area of genetic discontinuity that could be interpreted as a barrier to gene flow despite IBD, i.e. slope=0.005, p-value=0.005 (see section 3.1 in the Appendix).

Detection of a barrier to gene flow

Simulation test datasets. We found that MAPI was efficient in detecting strong to moderate barriers to gene flow, i.e. N e m≤1; F ST ≥0.1 with a barrier coverage ≥ 75% (see test datasets 2 in Table 2; Figures 3C-3D). When permeability increased (N e m>1; F ST <0.1), the method lost its accuracy to detect the barrier, i.e. the proportion of undetected barriers became large (see test datasets 2 in Table 2; Figure 3B). However, our simulations showed that IBD increased the power of the method to detect the linear barrier, especially for high levels of gene flow. For example, for N e m=2.5, even This article is protected by copyright. All rights reserved.

weak IBD (g=0.750; slope=0.005) increased the barrier coverage from 7% to 79%. As a result, MAPI performed well in identifying barriers to gene flow given that N e m<10 (see test datasets 3 in Table 2; . When individual samplings were small and irregular (N=75), the performance of MAPI to detect a weak barrier to gene flow (N e m>1) decreased by 30 to 70% (Table 2), though in these situations there was still often a graphical signal for a barrier (see Figure S3 in the Appendix). As expected, increasing genotyping effort improved further the detection of barriers for high levels of gene flow (Table 2). For example, for N e m=2.5, increasing the number of loci from 10 to 100 increased the barrier coverage from 7% to 98% in absence of IBD and for N e m=10, 30% of the barrier was still recovered.

Biological dataset.

In line with previous results published on the plant virus sequence dataset, MAPI identified a major area of genetic discontinuity bisecting the study area from North to South (see section 3.2 in the Appendix).

Sensitivity to parameter settings. The eccentricity and number of cells have little effect on the detection of unexpected significant areas (i.e. Unexpected in Figure 5). The worst situation (i.e. 3%) occurred when using very narrow ellipses (e=0.999; see Figure S1 of the Appendix) combined with a very high number of cells (i.e. β = 0.15). The combination of narrow ellipses (e=0.999), which result in a lower number of ellipses intercepting each cell, and a very low number of cells (i.e. β = 0.75) decreased MAPI efficiency to detect a barrier to gene flow (i.e. Coverage in Figure 5; Figure S3 in the Appendix). An eccentricity value of 0.975 associated with a β value of 0.5 for regular sampling or 0.25 for irregular sampling ensures high barrier coverage, low rate of false positives and reasonable computational time. Lower values of eccentricity (0.8-0.95; inflated ellipses) provide same quality results as a value of 0.975 but at the cost of expended computational time and stronger smoothing effect (see broader aggregates of significant cells for a value of 0.9 in the Figure S3 in the Appendix).

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

The effect of the eccentricity parameter is further illustrated on virus data (Figure S8 in the Appendix). In all cases, the major area of genetic discontinuity was still uncovered by MAPI.

Assessment of landscape effects

Simulation test datasets. As expected, when the two habitats were distributed randomly, there was no obvious correspondence between the spatial variation in m w and landscape cell values (Figure 4; see also Figure S4 in the Appendix to visualise examples of MAPI graphical outputs). This was confirmed by the results of the Bayesian conditional autoregressive model which showed that most of the posterior probabilities for a habitat effect were not significant (Table 3 and Figure S5 in the Appendix). A difference between habitats was, however, still detected from a few simulated test datasets, especially when using random sampling. For the gradient and fragmented landscapes, MAPI graphical outputs showed a relative convergence between variation in m w cell values and spatial habitat configuration (Figure 4 and Figure S4 in the Appendix). Accordingly, the intercept was significantly less for the favourable habitat (∝ ) than for the unfavourable habitat (∝ ) for all datasets regardless of the sampling scheme and size for the fragmented landscape and for most of the datasets when large (500 genotypes) and random or gridded samplings were used for the gradient landscape (Table 3). When the sampling size was smaller (200 genotypes) the proportion of significant posterior probabilities for ∝ < ∝ decreased of 25 and 50% for the random and gridded samplings, respectively. Interestingly, under the gradient configuration, the balanced sampling led to very poor results and even, in a few cases, to false positives (i.e. ∝ < ∝ ) (Table 3 andFigure

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

Biological dataset. On ground beetle microsatellite dataset, MAPI analyses supported previously published conclusions by identifying a significant area of genetic discontinuity corresponding to a large open field for the forest specialist only (Figure 6). When using m w cell values in the Bayesian conditional autoregressive model with the proportion of trees as an explanatory variable, we found a significant negative relationship only for the forest specialist (posterior probability=1). This means that, on average, between-individual genetic distances were lower within highly forested areas. No significant pattern was found for the generalist (see section 3.3 in the Appendix).

DISCUSSION

In this work, we presented a new method to translate networks of pairwise relationships into variation surfaces. MAPI is essentially a smoothing procedure using the overlap between ellipses as a way to share information between spatial connections. The surface produced is a grid of cells that provide information on the average intensity of the pairwise relationships crossing at the cell locations. The variation of the cell values over the surface allows to localise areas where the majority of the crossing connections correspond to very low or very high pairwise values. The significance of these areas can be assessed using a nonparametric randomisation procedure.

MAPI can be applied to genetic data to detect areas of high genetic continuity and discontinuity. When using neutral markers (e.g. microsatellites) and genetic differentiation measures, continuous and discontinuous areas reflect areas where gene flow is the highest and the lowest, respectively. Variation in gene flow intensity can result from different processes such as spatial heterogeneity in population density or migration success (Richardson et al. 2016). Here, using controlled simulations, we determined that MAPI performed well to detect genetic discontinuities resulting from a physical barrier as long as it is not too permeable to gene flow. When gene flow was This article is protected by copyright. All rights reserved.

high, MAPI's performance was substantially improved with sampling effort in the number of loci and individuals. Under such spatial genetics models, the sensitivity analysis to eccentricity setting showed that a value of 0.975 ensured high barrier coverage, low rate of unexpected significant cells, high spatial resolution and reasonable computational time. This result might, however, differ in more complex scenarios than a linear barrier to gene flow (e.g. highly fragmented landscape).

A central feature of MAPI is its relative insensitivity to IBD, which is a critical issue for spatial and landscape genetics analyses (Guillot et al. 2009;Thomassen et al. 2010;[START_REF] Bradburd | Disentangling the effects of geographic and ecological isolation on genetic differentiation[END_REF]. In an ideal situation, with perfect regular sampling and no edge effects, all cells should theoretically display the same m w value under strict IBD. When analysing highly irregular samplings, the method detected a few unexpected significant areas of genetic discontinuity as the cells located within the spatial gaps were only informed from long-distance connections. These highly discontinuous areas were, however, small and generally located on the border of the study area. As for many spatial exploratory analyses, regular individual-based samplings are more likely to provide reliable and interpretable results (Oyler-McCance et al. 2013;[START_REF] Balkenhol | Landscape Genetics: Concepts, Methods, Applications[END_REF].

In landscape genetics, assessing the effects of environmental features on spatial patterns of genetic variation often requires going beyond the detection of barriers to gene flow. In this context, our controlled simulations and ground beetle datasets [START_REF] Brouat | Fine-scale genetic structure of two carabid species with contrasted levels of habitat specialization[END_REF] showed that overlaying MAPI graphical outputs on landscape layers may provide information on which environmental variables are potential candidates to explain observed genetic patterns (see Figure 4 and Figure S9 in the Appendix). These candidate variables can be further explored by using m w cell values in post-MAPI analyses. To illustrate this possibility we used a regression model accounting for spatial autocorrelation in m w cell values and fitted in a Bayesian framework. From our simulated datasets, the regression model successfully retrieved the expected relationships (i.e. smaller m w values within the favourable habitat) for the fragmented and gradient landscapes but not when the two habitats This article is protected by copyright. All rights reserved.

were randomly distributed. These results can be explained by the difference in the scale at which habitats are spatially autocorrelated in the three simulated landscapes. As long distance effects are limited in the computation of m w by using an inverse-area-weighting procedure, the cells mainly reflect the average intensity of the shorter connections crossing at one location. Consequently, significant relationships between m w cell values and landscape variables are expected to be detectable when the landscape is spatially structured in such a way that a cell is mostly influenced by short-distance connections occurring between samples located within a relatively homogeneous landscape aggregate. Contrarily, when the spatial scale of autocorrelation in landscape features is smaller than the resolution of the sampling, a cell reflects between-as much as within-habitat connections and habitat effect are likely to become undetectable (e.g. random landscape). In the datasets simulated under a fragmented landscape model, sampling size and strategy did not affect the detection of habitat effect. Under the gradient configuration, the habitat effect was well detected especially when a large number (i.e. 500) of genotypes were sampled following a random, or even better, a gridded sampling strategy. Contrarily, the balanced sampling strategy was inefficient to detect an effect. This result can be due to the multi-scale patterns of autocorrelation of the habitats that have different consequences. First, as the two habitats strongly interpenetrate, the central part of the study area looks like the random landscape model (i.e. small habitat patches).

Within this zone, the balanced sampling produces numerous between-habitats connections at short distance that blur the relationship between the m w cell values and habitat type. Second, over the whole study area, the scale at which the habitats are autocorrelated is quite large and, subsequently, autocorrelation in m w cell values can be expected to be large as well. In such a situation, increasing the neighbouring structure from two to five cells in the conditional autoregressive model significantly improves the detection of the habitat effect from the balanced sampling (i.e. mean posterior probability = 0.909 and 0.927; percentage of significant probability = For each simulated dataset and sampling scheme, we extracted cell-specific m w values and habitat type (1 or 2) to fit the following regression model: =∝ + , where ∝ is an intercept depending on the habitat type in cell i and error term ε i includes a structured and an unstructured component. Here are presented the posterior probability that ∝ <∝ (P(α 1 < α 2 )) averaged over 20 replicates ( ) as well as the proportion of P(α 1 < α 2 ) ≥ 0.95 (P 95% ) and ≤ 0.05 (P 5% ). Note that for all simulated test datasets, the regression model was well adjusted with a mean predictive p-value very close or equal to 0.5 and residuals not overly structured in space (see landscape genetics illustrations on the website).

NA

Landscape model

Sampling scheme

P

  . , DOI : 10.1111/2041-210X.12616

  Figure 3 to visualise the simulated lattice). A null model of one panmictic population without a

  2 and 2.3 of the Appendix. For each of the simulated datasets, we also investigated MAPI sensitivity to sampling effort by sampling 200 and 500 individuals under three sampling schemes : (1) random sampling, in which individuals were sampled anywhere in the landscape regardless to the habitat; (2) balanced sampling, in which an equal number of individuals was sampled from each habitat; (3) gridded sampling, in which the same

  is an intercept depending on the habitat type in cell i. To account for spatial correlations it was assumed that the error term included two components: a structured component (i.e. spatially correlated) and an unstructured component (i.e. independently distributed). The unstructured heterogeneity term was assumed to be centred and normally distributed while the structured heterogeneity term was assumed to have a conditional intrinsic Gaussian autoregressive (CAR) distribution[START_REF] Besag | Bayesian image restoration, with two applications in spatial statistics[END_REF] with first-and second-order neighbours as neighbouring structure. The model was fitted in a Bayesian framework with an INLA approach[START_REF] Blangiardo | Spatial and spatio-temporal Bayesian models with R-INLA[END_REF] using the R-INLA package(Martins et al. 2013) of the R software(R Core Team, 2015). The fit of the model to the data was assessed using posterior predictive checking(Gelman et al., 2004), and we systematically checked that residuals were not overly structured in space. Bayesian inference resulted in posterior densities for the parameters ∝ (favourable habitat) and ∝ (unfavourable habitat). As we used a genetic distance, we expect a high posterior probability for the intercept of the favourable habitat (∝ ) to be less than the intercept of the unfavourable habitat (∝ ). The significance was assessed with a threshold of 0.05. Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. (2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology and Evolution, 7 (12), 1463-1475. , DOI : 10.1111/2041-210X.12616

  Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. (2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology and Evolution, 7 (12), 1463-1475. , DOI : 10.1111/2041-210X.12616

  S5 in the Appendix). Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. (2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology and Evolution, 7 (12), 1463-1475. , DOI : 10.1111/2041-210X.12616

Figure 1 .

 1 Figure 1. Schematic view of MAPI components.

Figure 2 .

 2 Figure 2. Steps in the test for spatial structure. Computation of m w values for the observed (A1, A2)

Table 2

 2

to identify test datasets).

Table 2 . Performance of MAPI under different spatial models.

 2 Parameters for each spatial model are indicated, including: the product of the effective population size and migration rate per generation (N e m); the shape parameter of the geometric distribution (g); the number of simulated populations (Ksim); and the number of loci (Nloc). Performance is measured as the average percentage, over 25 replicates, of the barrier covered by significant discontinuous areas (Coverage), and the average percentage, over 25 replicates, of significant discontinuous areas having no contact with the barrier (Unexpected). Large sampling refers to the sampling of 200 genotypes distributed regularly on the lattice, and small sampling to the sampling of 75 genotypes drawn randomly from the lattice. The datasets with 100 loci are indicated in italics. Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. (2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology and Evolution, 7 (12), 1463-1475. , DOI : 10.1111/2041-210X.12616This article is protected by copyright. All rights reserved.

	not used	minimum distance
		set to 1

Table 3 . Performance of MAPI under different landscape models.

 3

This article is protected by copyright. All rights reserved.

ACKNOWLEDGEMENTS

The authors thank T. Jombart, G. Guillot, N. Verzelen, R Leblois and A. Estoup for helpful comments on previous versions of the manuscript as well as the associate editor and four referees for their very constructive remarks. The authors also thank C. Desbiez and C. Brouat for providing the virus and forest ground beetle data, respectively. K. Berthier acknowledges the INRA Metaprogram "Sustainable Management of Crop Health" for funding the COPACABANA project. This work is dedicated to late Serge Meusnier.

MAPI methodology MAPI is implemented as an open-source SQL extension within open source database PostgreSQL 9.x and Postgis 1996-2013, The PostgreSQL Global Development Group: http://www.postgresql.org/ -http://postgis.net/). Source code, user manual, and training datasets can be downloaded from: https://

DATA ACCESSIBILITY -Rodent microsatellite data: DRYAD entry doi.org/10.5061/dryad.jf7sn, -Plant virus sequence data: Genebank accession numbers EU660581-EU660590, HM044202-HM044215, -Ground beetle microsatellite data: Figshare entry https://figshare.com/articles/GroundBeetles_zip/3482804 and French National Institute for Agricultural Research, https://www1.montpellier.inra.fr/CBGP/software/MAPI/ (file: MAPI_examples.zip). -Simulated test datasets and R scripts: French National Institute for Agricultural Research, https://www1.montpellier.inra.fr/CBGP/software/MAPI/ (file: MAPI_examples.zip).

Version postprint

Comment citer ce document : Piry, S. (Auteur de correspondance), Chapuis, M. P., Gauffre, B., Papaix, J., Cruaud, A., Berthier, K. ( 2016). Mapping averaged pairwise information (MAPI): a new exploratory tool to uncover spatial structure. Methods in Ecology and Evolution, 7 (12), 1463-1475. , DOI : 10.1111/2041-210X.12616

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

Analysis of MAPI cell values conjointly with landscape variables is a data-driven exploratory approach that can help identify candidate landscape features. Adequate sampling schemes and analyses aiming to test for the effects of those candidate variables should be conducted at an appropriate spatial scale before drawing conclusions (Richardson et al. 2016).

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

-Etherington, T. R. (2011). Python based GIS tools for landscape genetics: visualising genetic relatedness and measuring landscape connectivity. Methods in Ecology and Evolution, 2(1), 52-55.

-Galpern, P., Peres-Neto, P. R., Polfus, J., & Manseau, M. (2014). MEMGENE: Spatial pattern detection in genetic distance data. Methods in Ecology and Evolution, 5(10), 1116-1120.

-Garroway, C. J., Bowman, J., Carr, D., & Wilson, P. J. (2008). Applications of graph theory to landscape genetics. Evolutionary Applications, 1(4), 620-630.

-Gauffre, B., Estoup, A., Bretagnolle, V., and Cosson, J. F. (2008). Spatial genetic structure of a small rodent in a heterogeneous landscape. Molecular Ecology, 17:4619-4629. -Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2014). Bayesian data analysis (Vol. 2). Taylor & Francis.

-Guillot, G., Leblois, R., Coulon, A., & Frantz, A. C. (2009). Statistical methods in spatial genetics.

Molecular Ecology, 18(23), 4734-4756. -Guillot, G., Mortier, F., & Estoup, A. (2005). Geneland: a computer package for landscape genetics.

Molecular

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

-Kelling, S., Hochachka, W. M., Fink, D., Riedewald, M., Caruana, R., Ballard, G., & Hooker, G. (2009).

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved.

Accepted Article Accepted Article

This article is protected by copyright. All rights reserved. discontinuous areas having no contact with the barrier (Unexpected, right panel). Eccentricity is on the x-axis. We considered large regular sampling schemes (200 samples, top row) and small irregular sampling schemes (75 samples, bottom row). We used a subset of ten datasets simulated under spatial genetics models (see the section 2.1 in the Appendix to identify test datasets). Figure 6. MAPI graphical output for the forest specialist species of ground beetle. Cell-specific m w values appear as a colour scale. The hatched area corresponds to the location of the significant area of genetic discontinuity. Variation in tree coverage follows a greyscale from a minimum of 11% to a maximum of 75%. Sampling locations are illustrated by white circles proportional to sampling sizes.


Dirk-Jan De Koning

Xinyue Zhang
email: xinyue@uga.edu

Daniela Lourenco

Ignacio Aguilar

Andres Legarra

Ignacy Misztal

Weighting Strategies for Single-Step Genomic BLUP: An Iterative Approach for Accurate Calculation of GEBV and GWAS

Keywords: genome-wide association, SNP window, WssGBLUP, BayesB, BayesC

Genomic Best Linear Unbiased Predictor (GBLUP) assumes equal variance for all single nucleotide polymorphisms (SNP). When traits are influenced by major SNP, Bayesian methods have the advantage of SNP selection. To overcome the limitation of GBLUP, unequal variance or weights for all SNP are applied in a method called weighted GBLUP (WGBLUP). If only a fraction of animals is genotyped, single-step WGBLUP (WssGBLUP) can be used. Default weights in WGBLUP or WssGBLUP are obtained iteratively based on single SNP effect squared ( u 2 ) and/or heterozygosity. When the weights are optimal, prediction accuracy, and ability to detect major SNP are maximized. The objective was to develop optimal weights for WGBLUP-based methods. We evaluated 5 new procedures that accounted for locus-specific or windows-specific variance to maximize accuracy of predicting genomic estimated breeding value (GEBV) and SNP effect. Simulated datasets consisted of phenotypes for 13,000 animals, including 1540 animals genotyped for 45,000 SNP. Scenarios with 5, 100, and 500 simulated quantitative trait loci (QTL) were considered. The 5 new procedures for SNP weighting were: (1) u 2 plus a constant equal to the weight of the top SNP; (2) from a heavy-tailed distribution (similar to BayesA);

(3) for every 20 SNP in a window along the whole genome, the largest effect ( u 2 ) among them; (4) the mean effect of every 20 SNP; and (5) the summation of every 20 SNP. Those methods were compared to the default WssGBLUP, GBLUP, BayesB, and BayesC. WssGBLUP methods were evaluated over 10 iterations. The accuracy of predicting GEBV was the correlation between true and estimated genomic breeding values for 300 genotyped animals from the last generation. The ability to detect the simulated QTL was also investigated. For most of the QTL scenarios, the accuracies obtained with all WssGBLUP procedures were higher compared to those from BayesB and BayesC, partly due to automatic inclusion of parent average in the former. Manhattan plots had higher resolution with 5 and 100 QTL. Using a common weight for a window of 20 SNP that sums or averages the SNP variance enhances accuracy of predicting GEBV and provides accurate estimation of marker effects.

INTRODUCTION

Genomic best linear unbiased prediction (GBLUP) and Bayesian methods are routinely used for genomic selection (GS) and genome-wide association studies (GWAS) in animal and plant breeding, and in human disease studies. The main objective of GWAS is to identify the quantitative trait loci (QTL) with major effect, whereas GS focuses on using genomic information to better evaluate the genetic potential of individuals. Genomic estimated breeding value and single nucleotide polymorphism (SNP) effect or variance can be easily derived from both GBLUP and Bayesian methods. GBLUP assumes a normal distribution for SNP effects, and calculates their effect from phenotypes and a genomic relationship matrix. When not all individuals in a population are genotyped, single-step GBLUP (ssGBLUP) is the method of choice [START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF]. According to [START_REF] Legarra | Single step, a general approach for genomic selection[END_REF], the extra information on non-genotyped animals, the ability to account for pre-selection, and the independency of pseudo-phenotypes are partially responsible for gains in accuracy over other genomic methods. Bayesian methods such as BayesA, BayesB, and BayesC [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Kizilkaya | Genomic prediction of simulated multibreed and purebred performance using observed fifty thousand single nucleotide polymorphism genotypes[END_REF][START_REF] Habier | Extension of the Bayesian alphabet for genomic selection[END_REF] are non-linear approaches that usually assign to SNP effect a heavy tailed prior distribution, and then sample from the posterior distribution via Markov chain Monte Carlo (MCMC). GBLUP and ssGBLUP usually assumes equal weights for all SNP [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Goddard | Genomic selection based on dense genotypes inferred from sparse genotypes[END_REF]. This assumption is biologically incorrect but makes the statistics robust by limiting the number of unknown parameters [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF]. Nonlinear methods such as BayesA and BayesB assume heterogeneous variances of SNP effects, with emphasis on the SNP that links to major effects [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Meuwissen | Mapping multiple QTL using linkage disequilibrium and linkage analysis information and multitrait data[END_REF]. The performance of these methods has been proven to be better than GBLUP approaches in simulation studies assuming a few QTL with large effects and many QTL with small or null effects [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF][START_REF] Meuwissen | Mapping multiple QTL using linkage disequilibrium and linkage analysis information and multitrait data[END_REF][START_REF] Lund | Comparison of analyses of the QTLMAS XII common dataset. I: genomic selection[END_REF][START_REF] Guo | Comparison between genomic predictions using daughter yield deviation and conventional estimated breeding value as response variables[END_REF]. However, experiences with real dairy cattle data indicate that Bayesian regressions have resulted in reduced accuracy because of ignoring SNP with small effects [START_REF] Cole | Distribution and location of genetic effects for dairy traits[END_REF][START_REF] Su | Preliminary investigation on reliability of genomic estimated breeding values in the Danish Holstein population[END_REF] and that GBLUP approaches performed well for the majority of the traits of interest in livestock, mainly because of their polygenic nature [START_REF] Hayes | Invited review: genomic selection in dairy cattle: progress and challenges[END_REF][START_REF] Vanraden | Invited review: reliability of genomic predictions for North American Holstein bulls[END_REF][START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Chen | Genome-wide marker-assisted selection combining all pedigree phenotypic information with genotypic data in one step: an example using broiler chickens[END_REF][START_REF] Forni | Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information[END_REF][START_REF] Wang | Genome-wide association mapping including phenotypes from relatives without genotypes in a single-step (ssGWAS) for 6-week body weight in broiler chickens[END_REF].

One way to account for locus-specific variance in GBLUP-based methods is to include different weights for SNP. If those weights are known, weighted GBLUP (WGBLUP) provides genomic estimated breeding value (GEBV) similar to those of a Bayesian procedure using the same weights [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF]. WGBLUP and WssGBLUP were developed to allow the estimation of weights within GBLUP or single-step GBLUP (ssGBLUP), respectively. Many studies have found advantages of WGBLUP-based methods compared to unweighted GBLUP [START_REF] Snelling | Partial-genome evaluation of postweaning feed intake and efficiency of crossbred beef cattle[END_REF][START_REF] Gao | Comparison on genomic predictions using three GBLUP methods and two single-step blending methods in the Nordic Holstein population[END_REF][START_REF] Tiezzi | Genomic prediction using a weighted relationship matrix to account for trait architecture in US Holstein cattle[END_REF]. [START_REF] Sun | An iterative approach for efficient calculation of breeding values and genomewide association analysis using weighted genomic BLUP[END_REF] proposed two iterative procedures for calculating weights in WGBLUP, where the iteration was used to mimic the MCMC sampling in which both the prior and the posterior distributions kept being updated. In the first procedure, the weights were calculated as w (i+1) j = û(i)2 j , where w (i+1) j is the weight of SNP j at iteration i + 1 and û(i) j is the effect of SNP j at iteration i. This procedure is effective for identifying top QTL but excessively shrinks small SNP effects; thus, the accuracy of GEBV is reduced. In the second procedure the weights were calculated as w (i+1)

j = û(i)2 j + t, where t = σ 2 g 2 m j=1 p j q j
, σ 2 g is the genetic variance; p and q are the major and minor allele frequencies at locus j, respectively, and m is the number of SNP. This procedure introduced a constant to avoid SNP with no effect and brought the accuracy of GEBV close to that by BayesC but yielded Manhattan plots with lower resolution.

Recently, it was found that assigning a common weight to markers on a chromosomal region yielded more accurate estimates of GEBV under GBLUP-based methods [START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF][START_REF] Su | Comparison of genomic predictions using genomic relationship matrices built with different weighting factors to account for locus-specific variances[END_REF]. [START_REF] Su | Comparison of genomic predictions using genomic relationship matrices built with different weighting factors to account for locus-specific variances[END_REF] used group-marker variance from BayesR as a weighting factor for GBLUP in the evaluation of a dairy cattle population. Reliabilities for 4 production traits and mastitis were, on average, up to 1% higher and bias was reduced by 11% when using the mean variance of 30-SNP window compared with single SNP weighting. [START_REF] Xu | Genetic mapping and genomic selection using recombination breakpoint data[END_REF] demonstrated improved predictability in diploid plant QTL mapping using an artificial bin of linkage disequilibrium (LD)-linked neighboring markers. [START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF] evaluated WssGBLUP with simulation data using

d i(t+1) = u 2 i(t) 2p i (1 -p i ) , where d i(t+1) is the weight of SNP i at iteration t+1, u 2 i(t)
is the effect of SNP i at iteration t, and p i is the minor allele frequency (MAF). They iterated either on SNP alone or on GEBV and SNP. The first procedure gave a good identification of top QTL, and the second procedure provided a slightly higher accuracy of GEBV compared to BayesB, but only at the second iteration. The same study also found that the correlation of simulated and estimated SNP effects was lower than the correlation between simulated SNP and a sum of the effect of 8 SNP clustered around simulated QTL.

We hypothesized that efficient SNP weighting methods could considerably increase accuracy of predicting GEBV and help to better estimate SNP effects under WssGBLUP models. Therefore, the objective of this study was to present new procedures to calculate SNP weights individually or for genomic regions (windows) that would improve genomic predictions at animal and SNP level in WssGBLUP compared to other genomic models.

MATERIALS AND METHODS

Data Simulation

To test our hypothesis, one trait with a mean of 1.0, a phenotypic variance of 2.0, and an heritability of 0.5 was simulated using QMSim [START_REF] Sargolzaei | QMSim: a large-scale genome simulator for livestock[END_REF]. All the genetic variance was caused by QTL. A total of 20 chromosomes were created with an average length of 82 cM and containing 45,000 evenly distributed SNP with MAF ≥ 0.05. Three scenarios were considered involving different numbers of randomly placed QTL (5, 100, and 500) to mimic simple traits defined by major effects and complex traits affected by numerous minor effects. The scenarios were named 5-QTL, 100-QTL, and 500-QTL. All QTL were selected among SNP, and their effects were sampled from a gamma distribution with a shape parameter of 0.4. Both SNP and QTL were bi-allelic with no overlap between their positions, and a mutation rate of 2.5 × 10 -5 was assumed for markers and QTL per generation per locus. The simulated population originated from a historical population with 1000 generations of random mating to create linkage disequilibrium (LD) and establish mutation-drift balance. After that, a recent population was established from 200 males and 2600 females out of the historical population. The simulation was carried out for 205 generations. In every generation, the above number of males and females were selected to mate and produce 1 offspring, forming a population with an effective size (Ne) of 743. For all the analyses, generations 200-204 were treated as a training population and generation 205 as a validation population, with 1240 and 300 genotyped animals respectively. The complete dataset contained 18,400 individuals in the pedigree, of which 13,000 were phenotyped and 1540 were genotyped. Phenotypes were the sum of general mean, true breeding values (sum of QTL effects) and random residuals. Each simulated scenario was replicated 10 times.

WssGBLUP and New Weighting Methods

The ssGBLUP uses SNP to construct the genomic relationships. The genomic relationship matrix in our study was created as in [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF]:

G = ZDZ ′ 2 p i (1 -p i ) ,
where p i is the MAF of SNP i, Z is a matrix of centered genotypes, and D is a diagonal matrix of weights, where d ii is the weight for SNP i. In regular GBLUP-based methods, D = I, which gives a weight of 1 to all SNP. [START_REF] Strandén | Technical note: derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[END_REF] derived conversions of GEBV into SNP effects under GBLUP models and [START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF] extended the approach to ssGBLUP. In this approach, weight for SNP is based on SNP effect, and the latter is a best linear unbiased prediction derived from GEBV:

û = DZ ′ G -1 ĝ,
where û is a vector of estimated SNP effects and ĝ is a vector of GEBV. The weight for SNP i can be calculated as u 2 i or u 2 i 2p i (1p i ), and in the current study we used the first way.

Improvements in the SNP weights can be obtained iteratively either by recomputing only the SNP effects or by recomputing SNP effects and GEBV [START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF]. The latter also improves predictions of GEBV and was chosen for this study. The default WssGBLUP calculates and uses SNP-specific weights. Following [START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF], the iterative steps in WssGBLUP were:

a) Set t = 1, D (t) = I and G (t) = ZD (t) Z ′ 2 p i (1-p i ) b) Compute ĝ(t) using the ssGBLUP approach c) Compute SNP effects as û(t) i = D (t) Z ′ G (t) -1 ĝ(t) d) Calculate SNP weight as d ii (t+1) = û2 i , for all SNP i e) Normalize D (t+1) f) G (t+1) = ZD (t+1) Z ′ 2 p i (1-p i ) g) t = t + 1 h) Iterate from b) until t -1 = 10
Whereas the approach derived by [START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF] calculates SNP effects based on GEBV which combines direct genomic value (DGV; or sum of SNP effects weighted by SNP content), parent average (PA), yield deviation (YD), progeny contribution (PC), and pedigree prediction (PP). Lourenco et al. (2015a) showed that DGV is a more appropriate starting point for calculating SNP effects, because genotyped populations may comprise animals with different levels of accuracy. Therefore, in the step (c) GEBV was replaced by DGV. For all the 5 new procedures, the changes were done in step (d) where SNP weights are calculated based on SNP effects. Therefore, instead of using d ii = û2

i as SNP weight, the 5 new procedures we propose are: 1) Constant: still using the concept of SNP-specific weights, the weight was calculated as

d1 ii = û2 i + c where c = max( û2 i(0)
) is a constant chosen as the top SNP weight in the first iteration; 2) Nonlinear A: using similar approach to VanRaden (2008), SNP-specific weights were calculated as

d2 ii = û2 i /v |s i -2|
where ν is a scale parameter standing for the departure from normality and s i is the number of standard deviations from the mean for each 2 p i (1p i ); 3) Largest window: uses the concept of SNP-window weights, with weights for a group of 20 SNP as

d3 i,i = d3 i+1,i+1 = . . . = d3 i+19,i+19 = max(u 2 i , . . . , u 2 i+19 )
where d3 ii with ii from i to i+19 was used as weight for all SNP in the window comprised of 20 SNP; it uses non-overlapping windows; 4) Mean window: the weight for a group of 20 SNP was calculated as

d4 ii = d4 i+1,i+1 = . . . = d4 i+19,i+19 = n i=1 u 2 i n
where n is the size of the non-overlapping window; 5) Summed window: the weight for a group of 20 SNP was calculated as

d5 ii = d5 i+1,i+1 = . . . = d5 i+19,i+19 = n i=1 u 2 i

Models and Computation

Quality control of genomic data retained SNP with call rates >0.9, minor allele frequencies >0.05, and departures from Hardy-Weinberg equilibrium (difference between expected and observed frequency of heterozygous) <0.15. An average of 36,000 SNP remained after quality control, and the average LD (r2 = correlation between loci pair) in the last generation was about 0.29. For WssGBLUP, the model included population mean, animal effect, and random residual error. All phenotypes (except for validation animals), pedigree, and genotypes were accounted for. GEBV was obtained by using the software BLUPF90 [START_REF] Misztal | BLUPF90 and related programs (BGF90)[END_REF] with the simulated variance components; SNP effect was obtained by postGSf90 [START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Wang | Genomewide association mapping including phenotypes from relatives without genotypes[END_REF]. A total of 10 iterations of BLUPF90 combined with postGSf90 were used for the default WssGBLUP and WssGBLUP with new weighting methods. Each run of postGSf90 updated weights for SNP, whereas each run of BLUPF90 used the updated weights to constructed G matrices and, consequently, to improve GEBV estimates.

For GBLUP, the model was exactly as the one used in WssGBLUP, however, only genotypes and phenotypes for genotyped animals (except validation animals) were used together with the simulated variance components. The software BLUPF90 [START_REF] Misztal | BLUPF90 and related programs (BGF90)[END_REF] was used assuming the animals were unrelated through the pedigree.

The model for BayesB and BayesC had similar effects to the model used for WssGBLUP, but the animal effect was replaced by SNP effects. The phenotypic information was estimated breeding value (EBV) de-regressed (EBV DP ) following the approach described in [START_REF] Garrick | Deregressing estimated breeding values and weighting information for genomic regression analyses[END_REF], assuming that 0.05 of the genetic variation was not explained by the markers. The marker effects were assumed to follow a normal distribution with u i ∼ N(0, σ 2 u i ), and σ 2 u i the variance of the i th SNP. The proportion of SNP having no effect was set to 50%, 90% or 99%. Priors for σ 2 u i were calculated as

σ 2 a 2 p i (1-p i ) ,
where σ 2 a is the total simulated genetic variance. Degrees of freedom for the SNP and residual variances were set to 4 and 10, respectively. The software GenSel [START_REF] Fernando | GenSel-User Manual. Mapping Genes for Complex Traits in Domestic Animals and Their Use in Breeding Programmes[END_REF] was used to calculate SNP effects under the BayesB and BayesC approaches. A Monte Carlo Markov Chain (MCMC) was run for 41,000-50,000 and 1000-10,000 iterations were discarded as burn-in. For BayesB, 10 Metropolis-Hasting iterations were run per MCMC. Estimates of SNP effects were based on the posterior means of the MCMC.

For WssGBLUP, the output for genotyped animals in the validation population was GEBV that contained DGV, PA, PP; whereas for GBLUP, GEBV contained only DGV. The output for BayesB and BayesC is in the format of SNP effect, which can be translated into DGV by multiplying a row vector of SNP content by a column vector of SNP effects. In this way, GEBV from Bayesian methods for validation animals contained only DGV.

Model Comparison

The 5 new weighting procedures for WssGBLUP were compared with GBLUP, BayesB and BayesC with the proportion of markers with no effect (π) set to 0.5, 0.9, and 0.99. The simulated values were used as benchmark. Two comparisons were made:

(1) accuracy of predicting GEBV defined as the correlation between true simulated breeding value (TBV) and GEBV in the validation population ("accuracy" hereinafter); (2) ability to detect simulated QTL, done by visual inspection of Manhattan plots and by the amount of genetic variance explained. We did not derive significance thresholds or P-value for the latter comparison method because the objective was to identify the ability of the new weighting procedures in tracking the position simulated QTL, disregarding the significance of their effect on the trait. In addition, several studies reported difficulty in obtaining a measure of SNP significance when using shrinkage methods or Bayesian methods [START_REF] Servin | Imputation-based analysis of association studies: candidate regions and quantitative traits[END_REF][START_REF] Wakefield | A Bayesian measure of the probability of false discovery in genetic epidemiology studies[END_REF][START_REF] Wakefield | Bayes factors for genome-wide association studies: comparison with P-values[END_REF].

RESULTS AND DISCUSSION

Estimation of GEBV

Tables 1-3 show accuracies of GEBV and standard errors for the 5 new procedures to calculate SNP weights in WssGBLUP under three QTL scenarios, along with default WssGLBUP. Table 4 shows the accuracies of GEBV for GBLUP, BayesB, and BayesC. The average accuracies of the 5 new procedures were 0.87, 0.80, and 0.77 under 5-, 100-, and 500-QTL scenarios, respectively, and the standard deviation among 10 iterations ranged from 0.02 to 0.07. With default weighting, the accuracy increased initially but declined after some iterations depending on the number of simulated QTLs. As the number of QTLs increased, the inflection point came earlier (on iterations 4, 3, and 2 for 5-, 100-, and 500-QTL scenarios, respectively). The decline in accuracy with iteration was the result of continuously adding weight to SNP with large effects while shrinking SNP with small effects. Consequently, the accuracy of GEBV gradually decreased with iteration because the number of SNP with small effect increased.

For early iterations (≤5), largest, mean, and summed windows were the most accurate methods. Using the mean weight of 20 SNP in the windows as weight for the entire window improved accuracy by 0.01-0.09 compared to the other new procedures under the 5-QTL scenario. For 100-QTL, the procedure that added a constant outperformed the other new procedures, whereas for 500-QTL both constant and nonlinear A achieved the greatest accuracies that persisted for 10 generations. Overall, we observed that window procedures performed better than procedures with single SNP weights because the uncertainty was smaller [START_REF] Su | Comparison of genomic predictions using genomic relationship matrices built with different weighting factors to account for locus-specific variances[END_REF] and it also avoided extremely small values for SNP weights. A window size of 20 SNP was chosen over 5, 10, 50, and 100 based on accuracy (results not shown). Many factors, including the size of the reference population and population structure, influence the optimum window size [START_REF] Su | Comparison of genomic predictions using genomic relationship matrices built with different weighting factors to account for locus-specific variances[END_REF]. Window procedures maintained high accuracy with 5 QTL but lost the superior performance in late iterations with more QTL. The largest window scenario decreased in accuracy fastest among all window procedures, especially under the 500-QTL scenario, because it gave the greatest weight to the windows with large SNP effects and least weight to those with small SNP effects. This overand under-weighting introduced bias into the solutions. In regard to real genetic evaluations of massive data, the performance of iterations higher than 3 may not matter because one iteration usually takes from several hours (simple models) up to several days (complex model and population structure). The greatest accuracy at the earliest iteration is in bold face.

One of the purposes of investigating new weighting methods was to avoid sudden drops in accuracy after the peak was reached. We observed that most of the new methods were able to fulfill this requirement, especially when the number of simulated QTL was small. When the number of QTL increased, the method that adds a constant to û2 was the only one able to hold prediction accuracies over iterations. The main reason is that the constant value was chosen as the greatest û2 in the first iteration and remained the same for all iterations. This indicates that the shrinkage occurring along iterations is an important cause for the drop in reliability when SNP weights are calculated iteratively. The best constants for 5-, 100-, and 500-QTL scenarios were 8, 40, and 13, respectively, as they showed the highest accuracy averaged among 10 iterations compared to other tested constants. These values avoided SNP with no effects, which could reduce the accuracy of GEBV, while not deviating large effects significantly. Although adding a constant did not give as high an accuracy at early iterations as the window procedures, the accuracy remained stable after the peak was reached at iteration 5 (for 5-and 100-QTL). The accuracy for the default procedure was exceeded by the constant only in the 100-QTL scenario (+0.01). Adding a constant to avoid under-weighting in the 5-QTL scenario, where most SNP did not have effect on the trait, was redundant and counterproductive. Still, in the procedure that added a constant, the plateaued accuracies exceeded GBLUP by 0.11, 0.06, and 0.14 under the 5-, 100-, and 500-QTL scenarios (Tables 1234). On average, the increase was greater than the 0.06 reported in [START_REF] Sun | An iterative approach for efficient calculation of breeding values and genomewide association analysis using weighted genomic BLUP[END_REF] that used WGBLUP in a similar simulation of 10,000 SNP and 33 QTLs. The drawback of this approach is that the mechanism behind picking the right constant is still empirical and unclear; e.g., the average genetic variance t derived from GBLUP in [START_REF] Sun | An iterative approach for efficient calculation of breeding values and genomewide association analysis using weighted genomic BLUP[END_REF] was possibly too small for ssGBLUP. Theoretically, a threshold between zero and the peak SNP effect increases the bottom line of the absolute value for SNP with no effects. This threshold should both guarantee high accuracy of EBV and differentiate SNP effects. Number and distribution of QTL effects are related to this threshold, but in reality these are unknown.

The nonlinear A weighting method used in our study was adapted from VanRaden (2008) who defined the weight of SNP i as 1.25 |s i -2| , where s i is the number of standard deviations from the mean, and 1.25 represents the departure from normality. In our study, s ranged from 1.06 to 1.12. This procedure gave more weight to SNP with smaller effects, thus preventing the drastic decrease in accuracy. Its inferior performance compared with other procedures occurred for two reasons. First, oligogenic traits are correlated with few large QTL and the mean effect of all QTL is close to zero. Secondly, nonlinear A assigns more weight to SNP with effects but not to those with no effects; thus, it introduces bias into GEBV. This study showed that nonlinear A only performed as well as some of the other procedures (maximum accuracy of 0.81) under the 500-QTL scenario, whereas for 5-QTL scenario it had the lowest accuracy among WssGBLUP approaches.

The accuracies of GEBV from WssGBLUP were compared with those from Bayesian methods (Tables 1234). Except for the 5-QTL scenario, all WssGBLUP procedures under all scenarios surpassed BayesC and BayesB in accuracy before iteration 7. Accuracies for BayesB and BayesC π = 0.99 were the same under the 5-QTL scenario. For this scenario, BayesB with 3 different values of π and BayesC π = 0.99 outperformed WssGBLUP with nonlinear A weights. BayesC with π = 0.99 was 7 percentage points lower compared with the peak accuracy of the default procedure under the 100-QTL scenario (0.76 vs. 0.83) and 13 percentage points lower under the 500-QTL scenario (0.68 vs. 0.81). The biggest difference between default WssGBLUP and BayesB was under the 500-QTL scenario when π was 0.99 (0.81 vs. 0.41). This is consistent with previous studies [START_REF] Daetwyler | The impact of genetic architecture on genome-wide evaluation methods[END_REF][START_REF] Van Den Berg | Across population genomic prediction scenarios in which Bayesian variable selection outperforms GBLUP[END_REF], which indicated that Bayesian methods perform well when the number of QTL is small, whereas GBLUP-based methods perform better when the number of QTs is large, because it assumes an infinitesimal model. Single step GBLUP also includes pedigree relationships that contribute to the accuracy of GEBV [START_REF] Legarra | A relationship matrix including full pedigree and genomic information[END_REF][START_REF] Aguilar | Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF], whereas DGV from Bayesian methods exclude parent average [START_REF] Garrick | Deregressing estimated breeding values and weighting information for genomic regression analyses[END_REF]. The accuracy of GBLUP (Table 4) that had both pedigree and parent average removed was most of the times outperformed by BayesB or BayesC. [START_REF] Sun | An iterative approach for efficient calculation of breeding values and genomewide association analysis using weighted genomic BLUP[END_REF] also found that accuracies for BayesB and BayesC were 2.5 and 9.9% higher compared with GBLUP. In our study, posterior variances from BayesB and BayesC were used to weight G in WssGBLUP, but it did not improve the accuracy of GEBV (results not shown). When SNP variance from the literature was used as SNP weight for GBLUP, [START_REF] Zhang | Improving the accuracy of whole genome prediction for complex traits using the results of genome wide association studies[END_REF] showed considerable increase in accuracy for most of the analyzed traits in dairy cattle. [START_REF] Lourenco | Accuracy of estimated breeding values with genomic information on males, females, or both: an example on broiler chicken[END_REF] showed that GEBV from ssGBLUP is composed of parent average, phenotypic information, contribution from progeny, DGV, and pedigree prediction (pedigree relationships among genotyped animals). The accuracy of GEBV depends, therefore, on the weights each of these components receive. On the other hand, the accuracy from BayesB, BayesC, and GBLUP depends only on DGV, but if the same variance is used for all SNP in GBLUP, weights for DGV will not be optimized.

When methods that use SNP window weighting (largest window, mean window, and summed window) were compared to single SNP weighting (constant and nonlinear A), we observed an average increase in accuracy of 5 percentage points for 5-QTL (Tables 1234). For 100-and 500-QTL, there was a decrease in accuracy when window weighting was applied (0.80 vs. 0.81 and 0.74 vs. 0.81, respectively). When looking at the third iteration, which was the inflection point for most of the SNP window methods, the difference between windows and single SNP increased to 8 percentage points in 5-QTL. [START_REF] Hassani | Accuracy of prediction of simulated polygenic phenotypes and their underlying quantitative trait loci genotypes using real or imputed whole-genome markers in cattle[END_REF] observed that a SNP window of 1 SNP had the worst performance in terms of accuracy, and an increased SNP window helped to better capture the QTL signal when few QTL are simulated. In our 5-QTL scenario, largest and mean window of 20 SNP had better performance than the default weighting method, especially after the second iteration, although the difference was small (1% point). Mean window 20 was also able to maintain accuracy along the 10 iterations.

Differences between the first iteration (D = I) and the iteration where the peak accuracy was reached reduced as the number of simulated QTL increased, confirming that weights are more important when traits are influenced by fewer QTL. Thus, for the majority of the traits of interest in livestock breeding, ssGBLUP is able to attain the greatest accuracies without using weights for SNP. For specific traits influenced by few QTL or traits known for being influenced by important QTL (e.g., DGAT1 in milk fat and milk protein in dairy cattle) using weights that are common to a SNP window may help to better capture the signal from the QTL in that region. Therefore, using weights for SNP windows helps to improve the accuracy when a relatively small number of markers lose the ability to capture the possible QTL effects. Possibly, differential SNP weighting better depicts the covariance structure of G and reflects the actual contribution of SNP to additive genetic variance of a trait [START_REF] Su | Comparison of genomic predictions using genomic relationship matrices built with different weighting factors to account for locus-specific variances[END_REF].

QTL Identification

The average number of QTL with major effect that account for 50% of the total variance in 5-, 100-, and 500-QTL scenarios were 1.5 (0.5), 6.0 (2.0), and 26.1 (2.8). Figures 123show Manhattan plots for SNP effects in 5-, 100,-and 500-QTL. Figures 1A,2A, 3A show effects for default and proposed weights for WssGBLUP, whereas Figures 1B, 2B, 3B show effects for default WssGBLUP, BayesB, and BayesC. Under all QTL scenarios, using weights for SNP window produced clearer Manhattan plots than the single SNP weights. Although up to 20% of the simulated QTL were not correctly detected, most QTL with large effects were identified and few peaks were spurious associations. The good performance of window procedures resulted from assigning equal weights to neighboring SNP. This is especially instrumental to scenarios with many effects that are likely to be correlated with the majority of the QTL. The procedure with a constant reduced the difference between large and small SNP effects; hence, the plot was of lower resolution and had many small peaks compared to the default WssGBLUP. This was expected because SNP with small effect had a boost due to the extra weight added. A similar pattern was found in the nonlinear A procedure because the weighting factors had a limited range for all SNP. We observed that the shrinkage of SNP effects was bigger for the constant and nonlinear A methods under all scenarios. In the 500-QTL scenario, single SNP weighting caused an underestimation of SNP effects and SNP window weighting caused an overestimation, whereas default WssGBLUP returned SNP effects of the same scale as the simulated QTL. These results confirm that using weights for SNP when the trait is influenced by several QTL is not beneficial for better estimating SNP effects. If the trait is influenced by few QTL, using the maximum, average or sum of effects in a window of 20 SNP would help to capture clear signals from QTL. Bayesian methods, especially BayesB, best detected the simulated QTL under the 5-QTL scenario, however, with a different scaling of the SNP effect. The scaling of the SNP effect was very variable under the 5-QTL scenario. Under the 100-and 500-QTL scenarios, BayesB captured <1% of the SNP and assigned extremely high variance to them, especially in the 500-QTL scenario. This is because of the extreme shrinkage and avoidance of SNP with small effects, which caused bias in estimating SNP variances [START_REF] Calus | Accuracy of breeding values when using and ignoring the polygenic effect in genomic breeding value estimation with a marker density of one SNP per cM[END_REF]. Furthermore, BayesB and BayesC patterns depended on the choice of π.

Applying common weights for groups of SNP can help to increase the resolution of Manhattan plots; however, the resolution may be limited by the effective number of independent markers as proposed by [START_REF] Li | Evaluating the effective numbers of independent tests and significant p-value thresholds in commercial genotyping arrays and public imputation reference datasets[END_REF]. In their human genetics study, they showed that, when the ratio of effective number of markers to total number of markers is high (i.e., less redundant SNP), the threshold to declare a SNP as significant is more stringent, which indicates that the resolution worsens. [START_REF] Pocrnic | The dimensionality of genomic information and its effect on genomic prediction[END_REF] showed that the number of independent markers relates to the effective population size (Ne), where small populations share less independent markers. Perhaps, a Manhattan plot with high resolution (with noisy pattern) is a feature of smaller genotyped populations.

CONCLUSION

New procedures for calculating SNP weights in WssGBLUP can be effective in improving both the accuracy of GEBV and SNP effects. GEBV from WssGBLUP are more accurate than those from BayesB and BayesC, although different priors and π for the latter can change the ranking of the methods. Procedures that consider weights for a SNP window may be the best choice compared to single SNP given that the true number of QTL may not be known in real data. In addition, considering a group of SNP in the same genomic region may be the most appropriate way to capture the signal of an unknown QTL. The WssGBLUP method is especially useful for GWAS and GEBV estimation when only a fraction of the population is genotyped. WssGBLUP with 3 iterations may be enough to reach maximum predictivity of GEBV and SNP effects when a trait is influenced by few QTL.

FIGURE 1

 1 FIGURE 1 | Continued

FIGURE 2

 2 FIGURE 2 | Continued

FIGURE 3

 3 FIGURE 3 | Continued

FIGURE 3 |

 3 FIGURE 3 | Proportion of variance (%) explained by simulated QTL and SNP for different methods under 500-QTL simulation. (A) (a), true QTL; (b), default; (c), constant; (d), nonlinear A: weights as ν |s-2| , where ν is a scale standing for the departure from normality, and s is number of standard deviation from mean for each u 2 i ; (e), largest window; (f), mean window; (g), sum window. (B) (a), true QTL; (b), BayesC π = 0.5; (c), BayesC π = 0.9; (d), BayesC π = 0.99; (e), BayesB π = 0.5; (f), BayesB π = 0.9; (g), BayesB π = 0.99.

TABLE 1 |

 1 Average accuracy of genomic estimated breeding values (GEBV) and standard deviation across 10 replicates per iteration, for the 5-QTL scenario under the default weighting and 5 new weighting procedures for single-step Genomic Best Linear Unbiased Prediction (ssGBLUP).

	5-QTL

TABLE 2 |

 2 Average accuracy of genomic estimated breeding values (GEBV) and standard deviation across 10 replicates per iteration, for the 100-QTL scenario under the default weighting and 5 new weighting procedures for single-step Genomic Best Linear Unbiased Prediction (ssGBLUP).

	100-QTL

TABLE 3 |

 3 Average accuracy of genomic estimated breeding values (GEBV) and standard deviation across 10 replicates per iteration, for the 500-QTL scenario under the default weighting and 5 new weighting procedures for single-step Genomic Best Linear Unbiased Prediction (ssGBLUP).

	500-QTL

TABLE 4 |

 4 Average accuracy of genomic estimated breeding values (GEBV) for BayesB and BayesC under 3 simulation scenarios with varying π.

	Method	π a	5-QTL		100-QTL		500-QTL
			Accuracy	SE	Accuracy	SE	Accuracy	SE
	GBLUP	0.00	0.77	0.08	0.77	0.07	0.67	0.10
	BayesC	0.50	0.66	0.05	0.61	0.04	0.63	0.10
		0.90	0.75	0.06	0.67	0.03	0.65	0.09
		0.99	0.87	0.04	0.76	0.07	0.68	0.07
	BayesB	0.50	0.88	0.07	0.65	0.06	0.42	0.17
		0.90	0.89	0.07	0.68	0.07	0.41	0.16
		0.99	0.90	0.07	0.72	0.07	0.47	0.15

a π is the proportion of SNP that have no effect.

Frontiers in Genetics | www.frontiersin.org

August 2016 | Volume 7 | Article 151

ACKNOWLEDGMENTS

We thank Suzanne Hubbard and Sreten Andonov for revising the manuscript. Preliminary results of this study were presented at the10th World Congress on Genetics Applied to Livestock Production in Canada, 2014.

FUNDING USDA's National Institute of Food and Agriculture (Washington, DC; Agriculture and Food Research Initiative competitive grant 2015-67015-22936).

AUTHOR CONTRIBUTIONS

XZ implemented the research, contributed ideas to the algorithms, and wrote the article. DL contributed ideas to the algorithms and revised the article. IA programmed the original version of software, interpreted the algorithms, and revised the article. AL contributed ideas to the algorithms and revised the article. IM lead and designed the project, contributed ideas to the algorithms, and revised the article. All authors approved the final version to be published, and agreed to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.