Novel genomic approaches unravel genetic architecture of complex traits in apple

Background Understanding the genetic architecture of quantitative traits is important for developing genome-based crop improvement methods. Genome-wide association study (GWAS) is a powerful technique for mining novel functional variants. Using a family-based design involving 1,200 apple (Malus × domestica Borkh.) seedlings genotyped for an 8K SNP array, we report the first systematic evaluation of the relative contributions of different genomic regions to various traits related to eating quality and susceptibility to some physiological disorders. Single-SNP analyses models that accounted for population structure, or not, were compared with models fitting all markers simultaneously. The patterns of linkage disequilibrium (LD) were also investigated. Results A high degree of LD even at longer distances between markers was observed, and the patterns of LD decay were similar across successive generations. Genomic regions were identified, some of which coincided with known candidate genes, with significant effects on various traits. Phenotypic variation explained by the loci identified through a whole-genome scan ranged from 3% to 25% across different traits, while fitting all markers simultaneously generally provided heritability estimates close to those from pedigree-based analysis. Results from ‘Q+K’ and ‘K’ models were very similar, suggesting that the SNP-based kinship matrix captures most of the underlying population structure. Correlations between allele substitution effects obtained from single-marker and all-marker analyses were about 0.90 for all traits. Use of SNP-derived realized relationships in linear mixed models provided a better goodness-of-fit than pedigree-based expected relationships. Genomic regions with probable pleiotropic effects were supported by the corresponding higher linkage group (LG) level estimated genetic correlations. Conclusions The accuracy of artificial selection in plants species can be increased by using more precise marker-derived estimates of realized coefficients of relationships. All-marker analyses that indirectly account for population- and pedigree structure will be a credible alternative to single-SNP analyses in GWAS. This study revealed large differences in the genetic architecture of apple fruit traits, and the marker-trait associations identified here will help develop genome-based breeding methods for apple cultivar development.


Background
Until the end of the 20 th century, the lack of high throughput genotyping techniques and the limited development of high-density SNP arrays have hindered the advancement of genome-based breeding strategies for crop improvement. During the last 10 years, the genome sequences of about 20 plant species, including some from the Rosaceae family, were made publicly available [1], which offers opportunities for transforming breeding strategies to improve the yield and quality of major crops. Genome-wide association studies (GWAS) and genomic selection (GS) are among some new breeding tools proposed for crop improvement [2,3]. The underlying philosophy of both these strategies is to genotype enough markers across the genome so that at least one of the genotyped markers is likely in LD with the quantitative trait locus (QTL) alleles [4]. Both GS and GWAS can be conducted using the same genotypic and phenotypic data, but their objectives are different [3]. GS is used to predict phenotype from marker profiles alone, to reduce dramatically the length of the breeding cycle and the costs involved in phenotyping [5,6]. The objective of GWAS is to identify novel functional variation that can be deployed in cultivar development through markerassisted selection [2].
GWAS studies in humans have used two fundamentally different designs [7]: family-based and populationbased (that use unrelated individuals). The power of a GWAS of a quantitative phenotype using related individuals was shown to be slightly lower than that for a sample of unrelated individuals in a human study [8], but in crops and livestock controlled mating could make family designs more powerful than a population sample [4]. Both population-based [9] and family-based [10][11][12] designs have been used for GWAS in crops. The advantages of using relatives are manifold, including greater quality control, the ability to perform within-family tests of association that are robust to population stratification, and joint linkage and association analysis. A nested association mapping (NAM) population [10,11], which consists of multiple families derived from multiple inbred lines crossed to one or more reference inbred line, has been used for GWAS of different traits in maize. Multi-parent advanced generation inter-cross (MAGIC) population was first used in Arabidopsis [12]. Other family-based designs, such as parent-offspring, multigenerational pedigrees and multi-parent crosses, have historically been used in quantitative genetic studies. Thus, for plant populations, it is reasonable to consider large number of progenies derived from controlled crosses in various mating schemes for GWAS [4,11,13].
Population stratification and cryptic relatedness among studied individuals, when not taken into account, could lead to spurious genotype-phenotype associations in GWAS. Population stratification refers to the inclusion of individuals from different populations, while cryptic (or familial) relatedness refers to the presence of varying degree of genetic relationships between individuals within the study sample. GWAS methods based on the unified mixed linear model (MLM) were developed by [14] to account for population-level membership (to correct for structure) and individual-level relationships (to correct for cryptic relatedness). A realized individuallevel kinship matrix (G) calculated using molecular markers is more efficient than the pedigree-based kinship matrix (A) as it can account for Mendelian sampling and segregation distortion [4,15]. As family sizes in plant populations are much larger than those in other species, implementation of MLM was computationally very intensive. Therefore, the efficient mixed model association expedited (EMMAX) algorithm was developed to reduce this computational burden by reparameterizing the MLM likelihood function [16]. Further, a computationally more efficient and powerful compressed MLM (CMLM) that uses a group kinship matrix calculated from clustered individuals was developed recently [17]. Development of these methods has now made it much easier to analyze large amounts of data for GWAS. Unlike fitting each SNP in turn, which is a typical feature of GWAS, simultaneously fitting high density genome-wide SNPs could avoid the need to fit population and pedigree effects in MLM specifically [18].
In 2010, an international consortium published the first draft of the apple (Malus × domestica Borkh.) genome sequence using DNA from a popular apple variety 'Golden Delicious' [19], which led to re-sequencing of 27 apple cultivars that are the founders in global apple breeding programs. These efforts produced a huge reservoir of DNA markers, which helped the development of the first-generation apple Infinium SNP chip, comprising nearly 8,000 markers [20]. In the present study, we used this 8K SNP chip for GWA analysis of various fruit quality traits in a family-based design. Traits considered in this study relate to eating quality: fruit firmness (FF) and titratable acidity (TA); visual quality: red-flesh coverage (defined as weighted cortical intensity (WCI); see Methods); and susceptibility to physiological disorders: internal flesh browning (IB), bitter pit (BP) and fruit splitting (also termed cracking) (CR). To elucidate the relative contributions of different genomic regions, we implemented single-SNP analysis models, with and without accounting for population structure, and compared these with a model fitting all markers simultaneously. The statistical power of detecting SNP-trait associations was calculated using an expression derived in this study. The relative advantage of using realized relationships compared with pedigree-based expected relationships was also investigated. To our knowledge, this is the first large SNP array-based GWAS study to unravel the genetic architecture of quantitative traits for any major fruit crop.

Realized relationships and population structure
A plot of the first two principal components of the SNP genotypes data matrix grouped seedlings largely according to their familial relationships ( Figure 1). Some individuals did not cluster within their pedigree-assigned full-sib family groupings. For example, individuals in two families, namely A402 and A406, which have the same maternal parent, were clustered less tightly than the other five families. A break-away group of individuals from families A401 and A405, having the same maternal parent, apparently formed a separate cluster away from their respective full-sibs ( Figure 1). These patterns of clustering suggested some pollen contamination, so the actual number of pollen parents should be higher than that suggested by the mating design. Overall, a product-moment correlation of 0.65 was observed between pedigree-based (A matrix) and SNP-based estimates of pair-wise coefficient of relationships. The average pair-wise coefficient of relationships among all study individuals, obtained from the A and G matrices, were 0.36 and 0.50 respectively, reflecting that there are many more relationships not captured by the known pedigree records. The proportion of phenotypic variation explained R 2 LR À Á using the G matrix (in Equation 1) was higher than that using the A matrix for all traits (Figure 2). Results obtained after removing apparent contaminant seedlings, identified from PCA analysis ( Figure 1) and also by using PLINK software (http://pngu. mgh.harvard.edu/~purcell/plink), suggested that the magnitude of differences in R 2 LR values were almost identical (not shown) to those in Figure 2. Information presented in Figures 1 and 2 would suggest that using G would better account for population stratification than A, so only GWAS results (Equation 3) using G are presented here.

Linkage disequilibrium
The pattern of LD (r 2 ) decay in our GWAS population of 1120 individuals (arbitrarily called the first generation) was compared with that observed in the successive generation (i.e. second generation; see Methods). Results showed a high degree of LD even at longer distances between markers; for example, in the second generation the average r 2 for SNPs separated by 0.1 Mbp, 0.5 Mbp (approximately 1 cM in apple), and 1.0 Mbp was 0.28, 0.21, and 0.16, respectively ( Figure 3). This is somewhat lower compared to LD in the first generation (also reported earlier by [6]), but the pattern of LD decay was quite similar in both generations ( Figure 3).

Genome-wide associations
Goodness-of-fit of 'Q+K' (includes population structure and familial relationships) and 'K' (only familial relationships) models were compared to understand whether population structure could bias results. Different numbers of PCs of the SNP genotypes matrix constituted the Q matrix. The R 2 LR values of the 'Q+K' and 'K' models were identical for WCI, TA and BP, but were higher for 'Q+K' for the other three traits. Thus, the optimum number, as determined using the Bayesian information criterion (BIC), of PCs varied for different traits: 0 for WCI, TA and BP; 1 for IB; and 2 for FF and CR. However, results with or without incorporating Q in Equation 3 were not materially different, suggesting that accounting only for cryptic relatedness was sufficient to account for population stratification.
The profiles of p-values (in terms of -log 10 (p)) for all tested SNPs for each trait are illustrated in Figure 4. Uncorrected p-values of p < 5 × 10 -7 , which roughly equates to a genome-wide p-value of 0.00125 (= 2500 × 5 × 10 -7 ) as we tested a total of 2,500 SNPs, was used as    a significance threshold for individual SNP testing. The numbers of genome-wide significant SNPs detected for fruit firmness (FF), weighted cortical intensity (WCI), internal browning (IB), titratable acidity (TA), fruit cracking (CR) and bitter pit (BP) were 3, 36, 31, 18, 9 and 13 respectively. Most of the significant SNPs for any trait were clustered within a small genomic region, suggesting the presence of large-effect QTL at those positions. SNP-trait association signals for FF were identified on linkage groups (LG) LG3 and LG10; for WCI and IB on LG9 and LG16; for TA on LG8; for CR and BP on LG16 ( Figure 4). The SNP with the largest effect on FF was located on LG10, and this SNP is a T/G variant located within the first exon of the polygalacturonase (PG) gene (MDP0000232611), 20.833 kb from the top of LG10 (Table 1). A SNP with a massive effect R 2 LR À = 0.17) on WCI was located on LG9 ( Figure 4). This SNP on LG9 is a T/C variant and is located within the second exon of the MdMYB10 gene (MDP0000259616), 32.840 kb from the bottom of LG9. A cluster of SNPs with large effects on CR and BP, and moderate effects on WCI and IB, resides within the Leucoanthocyanidin Reductase (LAR1) gene (MDP0000376284) that is located between 1.496 kb and 1.669 kb on the top of LG16. The most significant SNPs described here are probably not the causative ones for our study traits due to extensive LD. For all traits except FF, the Q-Q plots ( Figure 5) showed a close adherence of the observed and expected -log 10 (p) values over most of the range, implying that the significant SNPs (highlighted in green colour) identified by the unified MMA are unlikely to be biased by population stratification.
The majority of SNPs individually explained only a small proportion of phenotypic variation (≈ 0.5%), while the largest-effect SNP explained 2, 17, 7, 6, 9 and 11% of the phenotypic variation for FF, WCI, IB, TA, CR and BP respectively (data not shown). The joint contribution of genome-wide significantly associated SNPs was also investigated. Because of LD, there were many significant SNPs within a small genomic region, so only the SNPs with the largest test statistics within 5 Mb regions were chosen. The joint R 2 LR calculated by fitting the chosen SNPs together in Equation 3, for FF, WCI, IB, TA, CR and BP were 0.03, 0.25, 0.11, 0.07, 0.11 and 0.12 respectively (Table 2), suggesting some improvement over single-SNP analysis. Fitting all 2,500 markers simultaneously (via SNP-derived G) captured nearly all genetic variance (i.e. heritability) for most of the traits ( Table 2). Correlation coefficients between SNP allele substitution effects (ASEs) obtained from single-SNP analysis and all-SNP analysis were about 0.90 for all traits, and largesteffect SNPs were generally common to both methods ( Figure 6).
Genomic regions with significant effects on the two pairs of traits (WCI and IB; and BP and CR) were further investigated by comparing LG-level estimated genetic correlations (r g ) for these two pairs of traits. Results suggested that r g values for the LGs harboring common significant regions were relatively higher than those for other LGs. Some of these LG-level correlations were quite different in magnitude as well as direction from the whole-genome correlation (Figure 7).

Power of the GWAS
The power of detecting marker-trait association for various QTL allele frequencies and trait heritabilities is shown in Figure 8. A LD value of 0.25 between a marker and QTL allele was assumed. For an unrelated sample size of 1120 individuals (the same size as in our study), the power of detection of an association with a locus explaining 2% of the phenotypic variation was 0.78, when marker and QTL allele frequency were 0.50. The power increased with reductions in marker and QTL allele frequencies. For a fixed sample size, the power of detecting SNP-trait associations declined with increasing relatedness among study individuals, but loss of power was minimal. The effect of trait heritability became more evident as the degree of relatedness increased (Figure 8).

Realized coefficient of relationships
The availability of genome sequence, the abundance of DNA markers, and high throughput genotyping platforms are providing a range of applications of molecular markers, including pedigree reconstruction, estimation of genetic parameters, and understanding relationships between genotype and phenotype in various species [5,[21][22][23]. Using a likelihood ratio based parameter R 2 LR À Á , our study showed that using SNP-based realized relationships in MLM could provide a better goodnessof-fit than using pedigree-based expected relationships. Our results also showed that for all traits, except IB, fitting all markers simultaneously could explain most or all the trait heritability. Similar results have been reported in studies on humans [23] and animals [24,25]. Also, approximate standard errors of G-based estimates of h 2 were considerably less than those for A-based because the former captured genetic relationships that are not accessible from pedigree records. Similar to studies on animals [22], our results suggest that the accuracy of artificial selection in plants species can also be increased by using more precise marker-derived estimates of genetic parameters.

SNP-trait genome-wide associations
With only a little loss of power, family-based designs in association studies offer various advantages compared with population-based designs [7,8]. The average SNPderived pair-wise coefficient of relationship in our study was about 0.50, but the loss of power to detect SNP-trait association was small (about 0.05) compared with that for an unrelated sample of the same size. Various factors such as sample size, high LD, minimal effect of population structure, and the presence of some largeeffect QTL provide high confidence in the SNP-trait associations identified in this study. The SNP array used in this study was designed to encompass SNPs in the coding region of predicted gene models and some candidate genes such as MdMYB10, MdPG, and MdLAR [20]. Peak association signals for FF, WCI and TA were located close to genomic regions that have been previously identified in bi-parental QTL mapping studies. The SNPs showing the largest effect on FF on LG10 (Figure 4) reside in the polygalacturonase (PG) gene, which depolymerizes cell wall pectin, and the involvement of this gene in the fruit softening process has been previously demonstrated in apple [26].
Red color in apple flesh results from a high concentration of anthocyanins. The role of the MdMYB10 gene on anthocyanin biosynthesis in red-fleshed apple was demonstrated using various approaches [27], and this gene has been mapped to LG9 [28]. The SNP marker Figure 6 Relationship between single nucleotide polymorphism (SNP) allele substitution effects obtained from single-SNP (y-axis) and all-SNP (x-axis) analysis. Correlation coefficient (r) is also shown for each trait. Estimates of narrow-sense heritability (h 2 ) obtained using pedigree-based relationships matrix (A). Approximate standard errors of h 2 A estimates varied from 0.12 (for TA) to 0.30 (for IB). 2 Proportion of the phenotypic variance explained jointly by genome-wide significant SNPs. 3 Simultaneously using all 2,500 markers in the form of the G matrix. Approximate standard errors of h 2 G estimates varied from 0.02 (for IB) to 0.04 (for WCI). associated with weighted cortical intensity (WCI) in our experiment is located in the second exon of MdMYB10, which is physically close to the motif in the MdMYB10 promoter that causally regulates transcription of MdMYB10 itself and thereby anthocyanin synthesis in apple flesh [27]. A cluster of SNPs at LG16 commonly associated with WCI and IB resides in the MdLAR1 candidate gene. LAR1 is a key enzyme in the flavonoid biosynthetic pathway, reducing leucoanthocyanidin into the flavanol compound catechin, a monomer of condensed tannins (also known as proanthocyanidins). Perhaps condensed tannins (CTs) act as co-pigments of cyanidin to create more intense red coloration in the fruit and hence the effect on WCI. Common genomic regions (especially MdMYB10 gene) found associated with WCI and IB are supported by results showing high genetic correlation (≈ 0.60) between these two traits [29]. Also, cold-stored fruit from all MdMYB10 transgenic lines of cultivar 'Royal Gala' showed varying degrees of symptoms of IB [29], suggesting a pleiotropic effect of MdMYB10 on WCI and IB. In our study, the estimated genetic correlation differed in degree as well as direction across different LGs, suggesting some possibility of breaking this undesirable correlation by means of carefully designed breeding and selection strategies, but further work on elucidating the genetic architecture of WCI and IB is required first.
The distribution of SNP effects for TA ( Figure 4) suggests one major QTL on LG8, supporting similar results from bi-parental QTL mapping studies [30,31]. However, there is no published report of candidate genes for TA on apple LG8. Segregation analysis approaches showed that inheritance of TA in a large apple population was better described by a mixed genetic architecture (a major gene and polygenes rather than polygenic or Mendelian inheritance [32]). Our results appeared to be in agreement with [32] in the sense that the largesteffect SNP accounted for about half the genetic variation, while the other half was accounted for by small-effect genes.
CR, which is a pre-harvest physiological disorder of apple, can be a serious economic problem for some cultivars [33]. However, little is known about the genetic architecture of this trait. Similarly to our results, estimate of h 2 from a previous study [34] indicated moderate genetic control of CR. While little is known on the physiological causes of genetic variability in CR susceptibility, it may be linked to the internal properties of the fruit during growth and differences in the elasticity of the peel when under particular stresses and strains caused by the developing parenchyma cells beneath it [33]. To our knowledge, there is no published report on marker-trait association for CR. A cluster of SNPs with a large effect on CR resides within the MdLAR1 candidate gene on LG16, which was previously reported to influence some other fruit quality traits such as astringency [6] and polyphenolic compounds [35].
BP is also a serious physiological disorder whose expression is generally observed in fruit after storage, but symptoms can also be observed on the fruit surface at harvest (sometimes called lenticel blotch [36]), and in our study they were classified as the same disorder. Genetic predisposition, calcium nutrition of the fruit, and environmental factors influence incidence of BP [36][37][38]. Based on the segregation ratio of resistant to susceptible seedlings, it was hypothesized that resistance to BP is controlled by two major dominant genes, named Bp-1 and Bp-2 [36]. Different segregation ratios (e.g. 1:1; 2:1 and 7:1) of resistance to susceptible seedlings were observed in various families in our study, suggesting a complex nature to this trait, which is further supported by our results showing that GWA-significant SNPs explained only about half the observed genetic variation (Table 2). Interestingly, the same cluster of SNPs on LG16 showed association with expression of BP and CR, and the direction and magnitude of LG-level genetic correlations were similar. Molecular, physiological and biochemical pathways that commonly contribute to the expression of these two traits are poorly understood, but our study provides a genomic hotspot for further investigations.
MLM that concurrently fits all available SNPs has been adapted in recent GWAS in animals [24,25] following an earlier study [18] that showed that provided highdensity SNPs are fitted simultaneously, admixed populations can be used to obtain reliable SNP effects even if pedigree structure and population structure have not been explicitly modeled. High correlations between SNP ASEs obtained from the 'Q+K' model (Equation 3) and all-marker analysis (Equation 4) for all six traits in this study (Figure 7) reinforces the findings of [18].

Application of SNP-trait associations
In our study population, the average LD between SNPs separated by 500 kb was high (r 2 = 0.25) largely because of relatedness (e.g. full sibs and half sibs) among seedlings and small effective population size. It is not uncommon to find different LD structures in different types of plant populations within a species. For example, in maize, LD decays within 1 kb in land races, within 2 kb in diverse inbred lines, and can extend up to 500 kb in commercial elite inbred lines (reviewed by [4]). Preliminary results (not shown) from an unrelated set of 125 individuals from a diverse apple germplasm collection showed that for a given distance (say, 500 kb) between markers, the extent of LD was almost one-third (r 2 = 0.08) of that in this study. One practical implication of these results is that marker-trait associations identified in advanced-generation crosses may not be repeated in relatively less improved breeding material (e.g. diverse germplasm). One of the key goals of GWAS is to identify large-effect marker-trait associations that can be deployed through marker-assisted selection (MAS) in subsequent generations of cultivar development populations. However, in order to conduct MAS using these SNPs in successive generations, strong LD between marker and QTL must persist across generations. Generally, only a selected set of individuals from the current generation are used as parents for the subsequent generation. Selection will lead to changes in allele frequencies at marker and trait loci, potentially reducing the LD between two loci, similar to that observed in the second generation of our study material (see Figure 3). As a result, the efficiency of GWAS-associated SNPs could be lower in the following generations. Nevertheless, except for WCI, the significant SNPs jointly explained less than 50% of the trait heritability in this study, which raises a question of how much variation in a quantitative trait needs to be accounted by a marker (or two) so that it would be worthwhile for MAS to be applied. Such a MAS scheme is generally viewed as cost effective compared with the genome selection, but this scheme does not bypass the phenotypic evaluation stage because there could be some quantitative traits for which no significant SNPs are identified. On the other hand however, a small SNP assay comprising GWA-significant SNPs for the key breeding traits could be used for pre-screening of seedlings before further field evaluation. Such an approach will not reduce the length of the breeding/selection cycle, but will shift the mean of the selection population. In order to keep the accuracy of MAS similar to that in the generation where SNPs were identified, periodical recalibration of SNP effects would be necessary [39].

Conclusions
The use of realized relationship matrix will provide higher accuracy of estimated genetic parameters, resulting in increased accuracy of artificial selection. There are apparently major differences in the genetic architecture of various traits in this study, i.e. for traits with similar heritability the distributions of SNPs effect were very different. The majority of SNPs individually explained only a small proportion of trait variation, but fitting all markers simultaneously captured most of the trait heritability for majority apple fruit traits. These findings suggest that genome-based methods could potentially replace the traditional apple cultivar breeding methods.

Plant material and phenotypes
A set of four female parents and two pollen parents were crossed in a factorial (4 × 2) mating design. One of the crosses was unsuccessful, leaving seven full-sib families. Seedling numbers varied between families, ranging from 40 to 350, with a total population size of 1,200. Seedlings were planted into the orchard (Havelock North, New Zealand) in July 2008 using a randomized block design. Further details of this experiment such as, parents involved, orchard management, harvesting and fruit storage protocols, were reported earlier [6]. Six traits were evaluated on the fruit samples using instrumental, sensory, or visual assessment methods. Phenotypic assessments for all traits were repeated for two years. Fruit splitting, observed as radial cracks in the stem end of the fruit (CR), and bitter pit (BP) were scored visually as presence or absence, with BP symptoms also noted if present within the fruit after cutting. Fruit from each seedling were cut in half across the equator and the proportion of the cortex area that was red (PRA) and the intensity of the red (RI) (=0 (none) to 9 (highest)) were scored. A weighted cortical intensity (WCI) was then calculated (PRA × RI) as an estimation of the amount of red pigment in the fruit. The proportion of the cortex area showing symptoms of cold-store-induced internal browning (IB) was recorded. Assessment protocols for fruit flesh firmness (FF) and titratable acidity (TA) were described in detail in an earlier study [6]. Individual fruit measurements (FF, WCI and IB) were first averaged for each seedling, and the average performance of each seedling over two years was used for testing genotypephenotype associations.

SNP Genotyping and LD estimation
Details of genotyping protocols for our study population were reported earlier [6]. Briefly, SNP genotypes were scored using the Genotyping Module (version 1.8.4) of the Illumina® GenomeStudio software (Illumina Inc.). The reliability of each genotype call was measured using the GenCall score, and SNPs were subsequently discarded using a sequence of criteria in the following order: GenCall score at the 50% rank (50% GC) < 0.40; cluster separation (ClusterSep) < 0.25; more than 5% missing calls; segregation discrepancy. Finally a high quality set of 2,500 SNPs was retained, and BEAGLE 3.1 software [40] was then used for imputing missing SNP genotypes.
Before looking at marker-trait associations, we calculated pairwise LD between SNPs, as a surrogate to LD between markers and QTLs, to evaluate the extent of LD in the study population (arbitrarily called the first generation) described above. These LD patterns were compared with those in an another population (second generation) comprised of 1600 seedlings derived from an incomplete factorial mating between six paternal parents (identified from the first generation) and four maternal parents (identified from previous selections). The degree of LD was quantified with the parameter r 2 [41] estimated using GOLD software [42].

Realized versus expected coefficients of relationship
The expected coefficients of relationship (i.e. the A matrix) based on pedigree records were compared with their realized counterparts (G) obtained using all available SNPs following [43]. A product-moment correlation was calculated between the elements (i.e. pair-wise coefficients of relationships) of the A and G matrices. We also compared the goodness-of-fit of a mixed model using realized or expected relationships: where y is the vector of observed phenotypic values of n seedlings; μ is an intercept, 1 n is a vector of 1s; Z is the known design matrices relating to a, the unknown vector of random additive genetic effects with a~N(0, A σ 2 a ) or a~N(0, G σ 2 a ). The scalar σ 2 a is the additive variance and ε is a vector of independent random deviates with variance σ 2 ε . Using A or G in Equation 1, we calculated and compared the R 2 LR values (which represent a likelihood-ratio based value of phenotypic variance explained) as follows [44]: where logL M is the maximum log-likelihood from fitting Equation 1; logL 0 is the maximum log-likelihood from fitting the intercept-only model. In addition to comparing R 2 LR values, we also compared estimates of heritability (h 2 ) of each trait obtained using the A or G matrices. Equation 1 was fitted using ASReml software [45].

Marker-trait association analysis
The unified mixed linear model (MLM) approach [14] that accounts for multiple degrees of relatedness (population structure and cryptic relationships) was used: where β is an unknown vector containing fixed effects, including a genetic marker, population structure (Q), and the intercept; X is the known design matrices relating to β. All other effects are same as in Equation 1. Equation 3 was fitted using GAPIT software [46], which uses computationally efficient and powerful methods, such as EMMAX [16] and CMLM [17]. To avoid spurious associations that could arise from population structure, we included principal components (PCs) derived from the genotypic data matrix (n × m) as covariates (i.e. Q matrix). The optimal number of PCs was determined by forward model selection using the Bayesian information criterion as implemented in GAPIT. In Equation 3, each SNP was tested in turn using a t-test (H 0 : No additive association between the SNP and trait), and p-values were obtained. Uncorrected comparison-wise p-value of p < 5 × 10 -7 , which is generally accepted to represent very strong proof of genome-wide association [47,48], was used to identify significant markertrait associations for all traits. A quantile-quantile (Q-Q) plot, which is commonly used for scrutinizing the population stratification in GWAS, was used to assess how well the model used for marker-trait association (Equation 3) accounted for population structure and familial relatedness. In this plot, the negative logarithms of the p-values from Equation 3 were plotted against their expected value under the null hypothesis of no association with the trait.
To compare the relative contributions of each SNP, we used R 2 LR values obtained by fitting Equation 3 with and without each SNP.
Allele substitution effects (ASE) at each SNP obtained from a single-SNP analysis (Equation 3) were compared with those obtained from a model (e.g. RR-BLUP) that concurrently fits all available SNPs. In our case, RR-BLUP is theoretically similar to using a SNP-derived relationship (G) matrix in Equation 1, i.e. a~N(0, G σ 2 a ) [43,49]. So the BLUP of a (i.e.α ) obtained using the G matrix in Equation 1 was used to estimate the vector of SNP ASEs (α) following [50]: where p i is the frequency of the A allele at the i th SNP (assuming three possible genotypes at each SNP were scored as AA, AB and BB); q i = 1 -p i ; elements of the i th column of M are 2q i , q i − p i and − 2p i for AA, AB and BB genotypes at a SNP locus. Product-moment correlations between ASE of SNPs obtained from Equations 3 and 4 are reported in this study. Common genomic regions showing a significant effect on a pair of traits were further investigated by comparing LG-level estimates of between-traits genetic correlations (r g ). For this purpose, correlation between the estimated BLUP-BV (from Equation 1) for two traits was termed as genome-level genetic correlation. BLUP-BVs of each seedling for each trait was then decomposed into LG-level BVs by using SNP ASEs (from Equation 4) and seedling's SNP genotypes at each LG. These LG-level BVs were then used to estimate LG-level between-trait genetic correlations.

Estimating the power of GWAS
When the significance of marker-trait association is tested by using the regression of phenotype on the number of copies of a SNP allele, the power of detecting association can be predicted from the probability [51]: where t v (δ t ) is a random variable with non-central Student's t-distribution with v (= n-2) degrees of freedom and non-centrality parameter δ t . The expression of δ t presented in [51] is that for unrelated samples, so given the genetic relationships among our study individuals, we derived a modified expression of δ t as: where R (≈ 1 − r 2 h 2 (1 − h 2 )) is the ratio of the approximate non-centrality parameter for genetically related individuals versus unrelated individuals, assuming that resemblance between individuals is due to additive genetic effects [8]; r is the coefficient of relationship and h 2 is the narrow-sense heritability; Γ(.) is a gamma function; b and σ b are the regression coefficient and its standard deviation respectively. For derivations of b and σ b , refer to [51].