Pectin Methylesterase genes influence solid wood properties of Eucalyptus pilularis.

This association study of Eucalyptus pilularis populations provides empirical evidence for the role of Pectin Methylesterase (PME) in influencing solid wood characteristics of Eucalyptus. PME6 was primarily associated with the shrinkage and collapse of drying timber, which are phenotypic traits consistent with the role of pectin as a hydrophilic polysaccharide. PME7 was primarily associated with cellulose and pulp yield traits and had an inverse correlation with lignin content. Selection of specific alleles in these genes may be important for improving trees as sources of high-quality wood products. A heterozygote advantage was postulated for the PME7 loci and, in combination with haplotype blocks, may explain the absence of a homozygous class at all single-nucleotide polymorphisms investigated in this gene.

Wood is primarily composed of lignin, cellulose, and noncellulosic polysaccharides (including neutral and pectic polysaccharides), all of which interact to form a complex structure that provides the foundation for the massive stature unique to the woody life form known as "trees." Wood properties are highly variable within and between tree species, a necessary attribute for adaptation to different environments (Carlquist, 1977). Quantitative genetic studies have revealed that wood properties are under strong additive genetic control in representative species across the Eucalyptus genus (Raymond et al., 2008;Stackpole et al., 2010aStackpole et al., , 2010b. Because wood phenotypes are complex traits, it is expected that they will be controlled by many genes, each with a small effect on the phenotype (Neale and Savolainen, 2004). Despite a basic knowledge of the biochemical pathways involved in lignin and carbohydrate (cellulose and noncellulosic polysaccharide) biosynthesis (Hertzberg et al., 2001), the relationship between physical wood properties and the activity of genes remains poorly understood.
Association studies are a means to bridge the understanding between phenotypic traits and underlying variation in the DNA sequence at specific candidate genes or dispersed genome wide. Differing from quantitative trait locus mapping, association studies use diverse individuals from natural populations, where generations of recombination have resulted in short stretches of linkage disequilibrium (LD; Butcher and Southerton, 2007). The rapid decay of LD in natural populations of forest trees means that genotypephenotype associations are identified with the causal quantitative trait nucleotide (QTN) or sites in close proximity to the QTN. Therefore, markers developed by association studies should be informative for selection across the full diversity represented in tree breeding populations (Butcher and Southerton, 2007;Grattapaglia, 2007). The rapid decay of LD in forest trees also means that a high marker density would be needed for genome-wide association testing; for this reason, a candidate gene approach has been preferred (Neale and Kremer, 2011). Candidate gene-based association studies restrict analysis to regions of the genome most likely to influence target phenotypes. Good candidates for wood properties include genes involved in lignin and carbohydrate biosynthesis (Hertzberg et al., 2001). The preferred DNA markers for association studies are single-nucleotide polymorphisms (SNPs), which are the most abundant source of variation in the genome and underlie much of the morphological variation observed between individuals (Rafalski, 2002).
This experiment targets two genes of the Pectin Methylesterase (PME) gene family (EC 3.1.1.11) that catalyze the demethylesterification of the most abun- Table I. Principal components analysis illustrating the high degree of correlation among the 40 wood properties used in the association testing The top three principal components (represented) explained 46% of the total variance in all wood properties analyzed. Loading values greater than 0.6 indicate strong correlations within components and are presented in boldface underlined numbers to assist with interpretation; strong negative correlations within component 3 are presented in italic underlined numbers.  Wolf et al., 2009). The degree and pattern of methylesterification of HG have strong influences on the water-holding capacity of pectin gels (Willats et al., 2001). In plant genomes, PMEs are a large gene family where paralogs have various biological roles, including wood formation, vegetative and reproductive processes, as well as plant-pathogen interactions (Pelloux et al., 2007). In wood formation, PME genes play roles in cell wall formation, stiffening, and remodeling (Moustacas et al., 1991;Al-Qsous et al., 2004;Wolf et al., 2009). Up-and down-regulation of PME in transgenic hybrid Populus tremula 3 Populus tremuloides (poplar) suggest an important role in wood formation, where it is localized in expanding wood cells and is involved in mechanisms determining cell/ fiber width and length (Siedlecka et al., 2008).
The candidate genes investigated (PME6 and PME7) were originally identified in southern blue gum (Eucalyptus globulus) where both genes are expressed in the intermediate and basal stem tissues, as well as the phloem and xylem (Goulao et al., 2011). PME6 and PME7 were also noted as good candidates for association studies because they are influenced by balancing selection with a high number of transsubgeneric SNPs shared between the Monocalyptus and Symphyomyrtus subgenera as well as excessively high SNP densities. It was proposed that balancing selection in these genes would be reflected by observed differences in the quantitative wood properties between Eucalyptus trees (T.R. Sexton, unpublished data).
The model species for this association study is Eucalyptus pilularis (blackbutt), which is an important hardwood species endemic to Australia. Timber from natural forests and plantations of E. pilularis is used for building construction, flooring, plywood, poles, railway sleepers, veneer, and pulp (Boland et al., 2006). Association testing demonstrated that SNPs in PME6 and PME7 were associated with variation in economically important wood properties of E. pilularis.

RESULTS
Relationships among wood properties were examined using principal components analysis (Table I). As expected, all traits were highly correlated within three components; for example, cellulose, Glc, and density had high loading values (greater than 0.6) in component 1. In this analysis of 40 wood phenotypes, three components accounted for 46% of the cumulative variation and were termed (1) cellulose and pulp yield traits, (2) shrinkage and collapse traits, and (3) lignin traits. The correlations within components 1 to 3 were consistent with published accounts for other Eucalyptus species; for example, acoustic wave velocity and pulp yield are positively correlated in both Eucalyptus dunnii and E. pilularis (Raymond et al., 2008(Raymond et al., , 2010. Considering the full population data set (n = 561) and correction for kinship, a total of 29 associations were identified to be significant at the individual test level (Supplemental Table S1). Of these, 10 were significant when P values were adjusted for multiple comparisons (Table II). Of the eight SNPs examined, five were significantly associated with at least one trait, and two of these, PME7 position 391 and PME6 position 204, had significant associations with multiple wood properties correlated in principal component 1, 2, or 3.

PME7 Is Primarily Associated with Lignin and Cellulose Traits
The SNP at position 391 of PME7 was significantly associated with lignin, cellulose, and pulp yield (Table  II; Fig. 1; Supplemental Table S1), traits correlated in principal components 1 and 2. Oddly, only two genotypes were observed (G/G and C/G). It seems unusual that the homozygous C/C genotype was not observed, despite large numbers of individuals with the alternate genotypes; consequently, there was a significant departure from Hardy-Weinberg equilibrium (HWE) at this locus (Table III). The C/G heterozygotes had a significantly higher percentage of Table II. Significant associations reported in the full population of 561 E. pilularis trees between wood traits and SNPs in PME6 or PME7 Test values reported were derived from a general linear model using 1,000 permutations and were significant following correction for multiple comparisons (P-adj). The P-perm value indicates significance at the permutation normalized SNP-wise level, while the P-adj value indicates significance at the experiment-wise level (corrected for multiple comparisons). NA, Not applicable. cellulose and pulp yield ( Fig. 1; Table II). In contrast, the G/G homozygotes had significantly more lignin. Another significant association was identified between a SNP in PME7 50 bp downstream at position 441 and the modulus of rupture (MOR; Table II; Fig. 1; Supplemental Table S1). MOR is the point at which wood fibers fail/break due to bending and is a trait positively correlated with cellulose, density, and acoustic velocity in principal component 1 ( Table I).
The alternate homozygote (G/G) was also missing, but, as indicated by the lower minor allele frequency (MAF), the G/G homozygote at position 441 was rare and not sampled. This was supported by no departure from HWE (Table III). Relative to the A/A homozygote, the heterozygous A/G genotype had a significantly higher MOR ( Fig. 1; Table II).
PME6 Is Primarily Associated with Cellulose, Pulp Yield, and Shrinkage Properties The SNP at position 204 of PME6 was significantly associated with wood shrinkage during drying, par-ticularly with tangential shrinkage in blocks and cores at 17% moisture content (Table II; Fig. 2; Supplemental  Table S1). Shrinkage is a trait influencing wood defect and is correlated with collapse in principal component 2. Only two individuals were observed with a homozygous T/T genotype. These two individuals were excluded from the association testing because of the small sample size. However, their phenotypes were plotted with the other genotypes for visual comparison (Fig. 2). All genotypes were observed at ratios expected under HWE (Table III). Relative to the trees with the homozygous C/C genotype, trees with the heterozygous C/T genotype had significantly higher levels of shrinkage (Fig. 2).
Another significant association was identified between the SNP in PME6 at position 268 and the ratio of tangential to radial shrinkage at 12% moisture content in cores (Table II; Fig. 2; Supplemental Table S1). Compared with trees with the homozygous T/T genotype, the heterozygous A/T trees had significantly higher levels of shrinkage (Fig. 2). The A/A homozygote was not sampled because the A allele was rare, with a MAF of 1.3%, and the genotypes showed no departure from HWE (Table III).
The last significant association was between the SNP at position 456 of PME6 and average maximum collapse (Table II; Table S1). Collapse is also a trait causing wood defect and is also correlated with shrinkage in principal component 2. Relative to the homozygous C/C genotype, the rarer heterozygous C/T genotype was associated with significantly higher levels of collapse, a trend consistent in the full population and all three subsets (Fig. 2). The T/T homozygous genotype was not sampled because the T allele was rare, with a MAF of 2.3%, and this locus showed no departure from HWE (Table III).

Associations Demonstrate Expected Genotype-Phenotype Links
In this study, associations were identified between two PME genes and a range of solid wood properties in E. pilularis. Associations in PME7 were primarily with cellulose, pulp yield, and lignin, traits correlated in principal components 1 and 3. In contrast, the associations with PME6 were primarily with shrinkage and collapse, traits correlated in principal component 2. The associations between PME6 and PME7 with multiple wood traits are consistent with the current understanding of wood formation and the biochemical activity of PME (Hertzberg et al., 2001;Willats et al., 2001;Kirst et al., 2004;Goulao et al., 2011).
In transgenic poplar, up-and down-regulation of PME localized in expanding wood cells has an effect on cell/fiber width and length (Siedlecka et al., 2008). Longer and thicker fibers increase the proportion of cellulose. Consequently, pulp yield is also increased, along with acoustic wave velocity (Raymond et al., 2010). The association between PME7 and cellulose might be explained due to competition for common precursor molecules, as the biosynthesis of pectin, cellulose, and hemicellulose all occur in the carbohydrate pathway (Hertzberg et al., 2001). The inverse association between PME7 and lignin is not as straightforward, given that lignin and carbohydrate biosynthesis occur in independent biochemical pathways; however, this inverse correlation may be explained by competition for carbon between the lignin and carbohydrate pathways (Kirst et al., 2004).
The associations between PME6 and shrinkage may be attributed to the water-holding (binding) capacity of pectic polysaccharides. Pectins are highly hydrophilic polysaccharides and, when demethylesterified, form gels that bind large volumes of water. The most abundant class of pectinaceous polysaccharides is HG, which is demethylesterified by PME (Wolf et al., 2009). Differing degrees and patterns of methylesterification of HG are known to influence the water-holding capacity of pectic gels (Willats et al., 2001). Trees with more water bound in their cell walls would be expected to show higher levels of shrinkage, as their timber is dried.
In addition, pit membranes facilitate the movement of free water between cells and contain high concentrations of pectin. The relative abundance of acidic versus methylesterified pectins of pit membranes is known to differ between species with low versus high ionic effect (Gortan et al., 2011). It has been proposed that weakening, distortion, and pit aspiration as timber is dried can explain differences in the drying behavior (shrinkage and collapse) of sapwood in Norway spruce (Picea abies; Rosner et al., 2010).
Although six SNPs seemed to interact with wood properties in some way, the two most intriguing associations were with PME7 position 391 and PME6 position 204. Neither of these SNPs resulted in clear alterations to the theoretically translated protein sequence of these enzymes. Because of the rapid decay of LD in Eucalyptus (Thumma et al., 2005;Thavamanikumar et al., 2011), the QTN-causative polymorphisms are likely to be located in close proximity to the SNP positions where associations have been demonstrated. Previously, these genes were sequenced, and multiple nonsynonymous (amino Table III. Genotypes, counts, and population parameters for 11 SNPs from PME6 or PME7 in E. pilularis (n = 561) Significant departures from HWE (P # 0.001) are designated by asterisks. Columns are as follows: gene and SNP position; AA, amino acid change; genotypes and respective counts; ObsHET and PredHET, observed and expected heterozygosity; HWpval, P value for HWE test; %Geno, percentage data set completeness; MAF; and whether the E. pilularis (P) SNP is shared with E. globulus (G) or E. pyrocarpa (Y). SNPs from PME7 at positions 30, 149, and 153 were excluded from the association testing.  Figure 2. A to D, Associations between PME6 position 204 and core radial scans at 12% moisture content (A), core tangential shrinkage at 17% moisture content (B), block tangential shrinkage at 17% moisture content (C), and average maximum collapse (D). The T/T genotypes were excluded from the association testing but included here for visual comparison. E, Association with PME6 position 268 and core tangential-radial shrinkage ratio. F, Association with SNP PME6 456 and average maximum collapse. Error bars represent SE, using a = 0.05; three asterisks indicate significance following correction for multiple tests, and one asterisk indicates significance at the single-marker test level. Trees were segregated into genotype categories, and corresponding phenotypes were plotted. [See online article for color version of this figure. ] ing haplotypes is an alternative approach capable of genotyping in highly heterozygous and polymorphic genotypes. Sequencing is an attractive alternative to traditional genotyping by primer extension, because SNP discovery and genotyping can be achieved in a single step, and this approach is now cost effective in light of rapid developments in second and third generation sequencing.
Strong Association between the SNP at Position 391 of PME7 and Wood Quality The SNP in PME7 at position 391 was significantly associated with a number of traits that were highly correlated in the principal components: (1) cellulose and pulp yield traits, and (3) lignin traits (Tables I and  II; Fig. 1). The homozygous G/G genotype was significantly associated with high lignin but also inversely correlated with lower cellulose and pulp yield (Fig. 1). The effect size for each of these associations was low (r 2 , 0.019), with phenotypic gains from the population mean to the G/G genotype mean being greater than 2 SE (Table II). The effect size was larger in the smaller unrelated subsets, for example, associations with Glc and director (acoustic wave) velocity 2 in subset B (n = 101), and both had high r 2 values of 0.06 (Supplemental Table S1). This might be attributed to the Beavis effect, where phenotypic variances of small populations/samples tend to be overestimated (Beavis, 1998). Given that E. pilularis is a species primarily grown for solid wood products, increasing the frequency or fixing the G allele in the breeding population is desirable to increase lignin and density.
It is intriguing that the alternate homozygous genotype of PME7 position 391 was not observed, despite the large sample size and high numbers of the alternate genotypes. PME orthologs do have pleiotropic functions, such as breaking seed dormancy (Ren and Kermode, 2000). Therefore, one explanation for the absence of the C/C genotype in our sample may be delayed seed germination, where slower germinating seeds would be inadvertently removed from the population in the nursery. Alternatively, the absence of the homozygous C/C genotype could be explained by a lethal recessive allele.
LD and the resulting haplotype blocks might explain why the SNPs sampled in PME7 only have two genotype classes. Significant departures from HWE meant that haplotyping and evaluation of LD in PME7 would be inappropriate (Table III). Forcing the analysis with all five SNPs in this gene did indicate strong LD (data not shown).
In PME6, detection of rare alleles better explained the lack of a homozygous genotype class, where each SNP in this gene had a low MAF and showed no significant departures from HWE (Table III). This sampling of three rare alleles in PME6 was not unusual; many populations of forest trees have been found to have a high abundance of rare SNPs (with low MAF), and it is postulated that this is attributed to large population sizes and an out-crossing mating system (Eldridge et al., 1993;Novaes et al., 2008;Grattapaglia et al., 2011;T.R. Sexton, unpublished data).

Strong Association between the SNP at Position 204 of PME6 and Shrinkage
The most highly represented associations in PME6 were with a synonymous SNP at position 204 and a number of traits that were highly correlated into principal components: (1) cellulose and pulp yield traits as well as (2) shrinkage and collapse traits (Tables I and II; Fig. 2). The homozygous T/T genotype was only sampled in two trees and therefore excluded from the association tests. When the phenotypes were segregated into genotypic categories for visual analysis, the alternate C/C homozygote class showed the lowest levels of shrinkage (Fig. 2). The heterozygous genotype displayed an intermediate phenotype, suggesting that there was additive gene action at this locus. The effect size of each association at PME6 position 204 was low (r 2 , 0.023), with predicted gains less than 1 SE above the population mean for the C/C genotype for most traits (Table II). Like associations with other SNPs, the effect size was larger in the smaller unrelated subsets; for example, associations with core radial shrinkage at 12% moisture content in subset C (n = 84) had a high r 2 value of 0.11 (Supplemental Table S1), although again, this may indicate overestimation due to the Beavis effect (Beavis, 1998). Given that E. pilularis is a species primarily grown for sawn wood products, fixing the C allele in the breeding population would be desirable to reduce the shrinkage of drying timber.

Balancing Selection in PME and Tree Breeding
Previous studies have suggested that both PME6 and PME7 are under balancing selection; therefore, these genes are likely to have adaptive functions that help the tree to cope with biotic and abiotic stresses (T.R. Sexton, unpublished data). The significant departures from HWE and the corresponding high abundance of heterozygous genotypes in PME7 suggested that a heterozygous advantage may be the mechanism balancing selection at this locus. The mechanism of balancing selection in PME6 is not as clear because the SNPs genotyped are rare and showed no departures from HWE.
This study suggests that the high nucleotide diversity maintained by balancing selection in these PME genes (T.R. Sexton, unpublished data) is also associated with the structural composition of wood. Tree breeders should be mindful that these PME genes may be under balancing selection; therefore, reductions of diversity at these loci may leave a population more vulnerable to disease or abiotic stresses. This concept needs further investigation, but this paper has identified genes and specific alleles that could be used as models to initiate further studies.

Implications for Breeding and Deployment Systems
We have shown how SNPs in two PME genes influence many wood properties of E. pilularis. The variance (r 2 ) explained by these associations is low but typical for complex traits and all associations reported in other forest trees (Thumma et al., 2005(Thumma et al., , 2009González-Martínez et al., 2007Neale, 2007;Ingvarsson et al., 2008;Dillon et al., 2010;Sexton et al., 2010b). Selections based on individual SNPs are expected to provide only small gains (Grattapaglia and Resende, 2010). Nevertheless, small gains may be economically advantageous even if only representing a very minor shift from the population mean. As demonstrated in the livestock industry (Goddard and Hayes, 2009), investments in marker development can be rapidly recouped by phenotypic gains among the massive numbers of individuals descended from genotyped parents.

SNP Discovery
SNPs were identified in a pooled sample of 30 trees representing the natural geographical range of Eucalyptus pilularis. Both PME6 and PME7 were amplified from this pooled DNA and sequenced on an Illumina GAIIx (T.R. Sexton, unpublished data).

Population and Phenotyping
We investigated 561 destructively sampled trees, which were descended from 284 open pollinated families and grown in a field trial at Hannam Vale, New South Wales, Australia (31°40#S, 152°33#E). These families were collected from well-separated mother trees from 37 provenances of E. pilularis (Fig. 3) and represented a large part of the natural geographical distribution of this species (Johnson and Stanton, 1993;Johnson and Nikles, 1997;Boland et al., 2006). Each of the 561 trees was destructively sampled at 112 months and had undergone extensive phenotypic assessment of quantitative traits, including basic density, wood stiffness, wood chemistry (cellulose, lignin, and hemicellulose), wood shrinkage, and collapse during drying Raymond et al., 2008Raymond et al., , 2009, as well as hardness, modulus of elasticity, and MOR (Henson and Smith, 2007).

Principal Components Analysis
Many of the wood traits were highly correlated, so principal components analysis was performed (using PASW statistics 18) in order to establish relationships between traits (Table I).

Sequenom Genotyping
SNPs previously characterized in PME6 and PME7 (T.R. Sexton, unpublished data) were genotyped across these same 561 individuals for which phenotypic data had been collected. The availability of high-resolution SNP discovery data (T.R. Sexton, unpublished data) permitted use of the standard Sequenom iPlex Gold Assay for genotyping. Most SNPs identified in these genes by sequencing could not be genotyped by primer extension because of the high SNP density; ignoring this information would lead to mispriming and subsequent assay failure (Sexton et al., 2010a;Grattapaglia et al., 2011). Assays were designed based on previously ascertained E. pilularis sequences for PME6 and PME7 (T.R. Sexton, unpublished data). All SNP positions in these genes were coded in the International Union of Pure and Applied Chemistry ambiguity code, and candidate SNPs desired for genotyping were converted to the format required for the Sequenom assay design (3.1) software.
A total of 11 SNPs were successfully genotyped with more than 75% dataset completeness. Only eight SNPs were used for association testing, as three SNPs from PME7 at positions 30, 149, and 153 were excluded because limited numbers of individuals were represented in the alternate genotype class or classes (Table III).

Association Testing
Five SNPs from PME7 and three from PME6 (for a total of eight SNPs) were used for association testing with 40 wood-quality traits. The sample size was nominally 561 related trees, but numbers in each test may be lower because of missing data for either genotypic or phenotypic data sets.
Population structure is known to lead to false positives in association studies (Pritchard and Rosenberg, 1999). However, in this study of E. pilularis, correction for population structure was unnecessary and not applied because the association population represents provenances previously classified as one structure group, north of Sydney (Shepherd et al., 2010).
Because Eucalyptus trees are diploid, each SNP will have only two nucleotide alleles in an individual tree (the genotype), and within the population three possible genotypes may exist: the heterozygote and two alternate homozygote genotypes. These molecular data are often unbalanced where a populations contains unequal numbers of individuals corresponding to each genotype; therefore, a general linear model of TASSEL 2.1 was used for association testing. This model uses a S-restricted model that is a good test of marker effects, even when the data are unbalanced (Bradbury et al., 2007;Buckler et al., 2009).
As the full sample (n = 561) included trees that were known to be half-sibs, a kinship coefficient was used as covariate data in the model, partitioning out covariance due to relationship. The degree of genetic covariance caused by polygenic effects was defined as positive for a pair of related individuals and zero for a pair of unrelated individuals (Yu et al., 2006). The kinship matrix was producing with SPAGeDi (Hardy and Vekemans, 2002), using genotypes from 132 SNPs, following an approach described previously (Loiselle et al., 1995;Yu et al., 2006;Bradbury et al., 2007).
Association testing using eight SNPs and 40 traits meant that 320 comparisons were made, so permutation testing was used to correct for multiple comparisons. Following the guidelines of Churchill and Doerge (1994), 1,000 permutations were applied for a 0.05 significance level. Two P values are reported. The first is a SNP-wise value (P-perm), which is a normalized value following the permutations (Supplemental Table S1). The second P value reported is an experiment-wise permutation test value (P-adj), which is corrected for multiple comparisons and is derived from a step-down MinP procedure (Ge et al., 2003;Buckler et al., 2009).

Validation of Association Testing
To ensure that the results of this study were not biased by relationships among trees, associations identified in a large population of 561 related trees were validated in at least one of three smaller subsets composed of unrelated family members: subset A (n = 284), subset B (n = 101), and subset C (n = 83). All subsets are derived from the larger population; subsets B and C are open-pollinated siblings of subset A. Tree phenotypes were segregated into genotypic categories based on the allele combination, and graphs were plotted for each subset, including the mean and SE. All test values from the full population and the subsets are presented in Supplemental Table S1.

Calculation of HWE
All 561 trees were included for tests of HWE using the software package GEVALT (Davidovich et al., 2007). Significant departures from HWE in PME7 suggested that these SNPs were inappropriate for the calculation of LD (Table  III). Only three SNPs were genotyped in PME6, but they showed little evidence of LD between them (data not shown).

Avoiding Paralogs
Based on studies in other plant species, PMEs are coded by a large gene family (Pelloux et al., 2007;Wolf et al., 2009). Therefore, precautions were necessary to minimize the potential to inadvertently confound paralogous gene sequences. During SNP discovery, primer specificity was carefully crossvalidated against the Eucalyptus genome version 1.0, 83 assembly (http:// eucalyptusdb.bi.up.ac.za/; T.R. Sexton, unpublished data). However, this problem could also exist within the genotyping data, where SNP assays utilized different primer pairs that amplify much shorter products. The unusual observation of 100% heterozygosity in PME7 at position 149 could be attributed to the sampling of gene paralogs, so this locus was removed from the analysis (Table III). Despite this potential problem, SNPs used in the association testing were not attributed to paralogs, as demonstrated by the observation of homozygous individuals in the population as well as by evidence from a SNP in PME7 at position 391, where a controlled back-cross between parents (G/G 3 G/C) resulted in a 1:1 Mendelian inheritance among six progeny.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Table S1. Results from association testing, including the population, traits, corresponding SNP, P value, corrected P value, and R 2 value for each test.