*4.5. Phenotypic Variation Explained by All SNPs*

The phenotypic variation explained by all SNPs, denoted as *h* 2 *SNP*, was estimated for all traits based on the following mixed linear model [60] implemented in the GCTA software [61]:

$$\mathbf{y} = \mathbf{X}\boldsymbol{\mathfrak{f}} + \mathbf{g} + \boldsymbol{\varepsilon} \text{ with its variance } \mathbf{V} = \mathbf{A}\sigma\_{\mathcal{g}}^2 + \mathbf{I}\sigma\_{\varepsilon}^2 \tag{2}$$

where *y* is an *n* × 1 vector of phenotypes with *n* individuals in a population, *X* is the *n* × *n*.

structure matrix, *β* is a vector of fixed effects of population structure, including posterior probabilities of an individual assigning to a cluster in DAPC, *g* is an *n* × 1 vector of the total genetic effects of the individuals with *g* ~*N* (**0**, *Aσ* 2 *g* ), and *ε* is a vector of residual effects with *ε* ~*N* (**0**, *Iσ* 2 *ε* ). *A* is interpreted as the genetic relationship matrix (GRM) between individuals and estimated from SNPs. *σ* 2 *g* is estimated using the restricted maximum likelihood (REML) method based on the GRM estimated from all SNPs. Thus, SNP heritability *h* 2 *SNP* was estimated as

$$h\_{\rm SNP}^2 = \frac{\sigma\_{\rm \mathcal{S}}^2}{\sigma\_{\mathcal{S}}^2 + \sigma\_{\varepsilon}^2} \tag{3}$$
