Association Genetics of Wood Physical Traits in the Conifer White Spruce and Relationships With Gene Expression
Jean Beaulieu, Trevor Doerksen, Brian Boyle, Sébastien Clément, Marie Deslauriers, Stéphanie Beauseigle, Sylvie Blais, Pier-Luc Poulin, Patrick Lenz, Sébastien Caron, Philippe Rigault, Paul Bicho, Jean Bousquet, John MacKay

Abstract

Marker-assisted selection holds promise for highly influencing tree breeding, especially for wood traits, by considerably reducing breeding cycles and increasing selection accuracy. In this study, we used a candidate gene approach to test for associations between 944 single-nucleotide polymorphism markers from 549 candidate genes and 25 wood quality traits in white spruce. A mixed-linear model approach, including a weak but nonsignificant population structure, was implemented for each marker–trait combination. Relatedness among individuals was controlled using a kinship matrix estimated either from the known half-sib structure or from the markers. Both additive and dominance effect models were tested. Between 8 and 21 single-nucleotide polymorphisms (SNPs) were found to be significantly associated (P ≤ 0.01) with each of earlywood, latewood, or total wood traits. After controlling for multiple testing (Q ≤ 0.10), 13 SNPs were still significant across as many genes belonging to different families, each accounting for between 3 and 5% of the phenotypic variance in 10 wood characters. Transcript accumulation was determined for genes containing SNPs associated with these traits. Significantly different transcript levels (P ≤ 0.05) were found among the SNP genotypes of a 1-aminocyclopropane-1-carboxylate oxidase, a β-tonoplast intrinsic protein, and a long-chain acyl-CoA synthetase 9. These results should contribute toward the development of efficient marker-assisted selection in an economically important tree species.

ACQUISITION of a perennial woody growth habit has had an enormous impact on the evolution of plants and of life on land (Groover 2005). Trees have shaped many terrestrial ecosystems and provide wood and fiber used by diverse industries worldwide. The adaptive importance of wood (secondary xylem) is indicated by its remarkable diversity among species, its significant variation within species, and its developmental plasticity. Although many of the genes involved in wood growth and development are known (e.g., Whetten et al. 2001; Aspeborg et al. 2005; Pavy et al. 2008a), the genetic basis underlying the variability of wood structure is only partly understood.

Revealing the genetics of wood and growth traits in forest trees initially proceeded through mapping quantitative trait loci (QTL) (Groover et al. 1994; Bradshaw and Stettler 1995; Grattapaglia et al. 1996; Plomion et al. 1996; Sewell et al. 2000; Wheeler et al. 2005; Ukrainetz et al. 2008). Major QTL were identified and reported to explain up to a total of ≥25% of the naturally occurring variation depending on the species and the trait (Kirst et al. 2004). One limitation with linkage mapping is its requirement for large families with known relatedness, i.e., full-sib families in plants or extended pedigrees in humans (Myles et al. 2009). Furthermore, QTL do not replicate well across populations of different genetic background and across environments (Mackay 2001; Pelgas et al. 2011; Ritland et al. 2011).

In contrast to QTL mapping, genetic association mapping [or linkage disequilibrium (LD) mapping] can readily be applied to natural or breeding populations to identify marker–trait associations that remain linked following a large number of recombinations over a population's history. LD mapping has led to research advances in human, animal, and plant genomics and is increasingly being adopted as the method of choice to dissect quantitative traits (Hirschhorn et al. 2002; Cheung et al. 2005; Zhu et al. 2008; Goddard and Hayes 2009; Grattapaglia et al. 2009; Rafalski 2010; Neale and Kremer 2011). By utilizing past recombination events, association mapping increases the resolution of marker–trait associations. It is advantageous for species with long generation times because it may capture a broader sampling of allelic variation in a single study. Its application in undomesticated species like forest trees also holds potential to enhance and accelerate genetic resource management activities, including gene conservation and marker-assisted selection (MAS) for quantitative traits (Haussmann et al. 2004).

Two approaches can be taken in association mapping to identify marker–trait associations. Randomly chosen markers may be distributed evenly across the genome to capture associations through LD with the genes controlling the traits. The alternative is to select a subset of candidate genes that are hypothesized to contribute to the quantitative variation on the basis of prior knowledge. Because conifers have immense genomes (Murray 1998) and low LD within gene limits (Neale and Savolainen 2004), genome-wide coverage of randomly selected markers is not practical. On the other hand, the rapid accumulation of cDNA sequences and expression data in conifers (e.g., Pavy et al. 2008a) enable large-scale candidate gene approaches.

Association studies of wood traits in a handful of tree species suggest the approach is promising but results have been species dependent. In Eucalyptus nitens (H. Deane & Maiden) Maiden, allelic variation in a cinnamoyl CoA reductase gene explained a small but significant proportion of the variation in the cellulose micofibril angle (Thumma et al. 2005), and a variant of a cobra-like gene was linked to cellulose content and pulp yield (Thumma et al. 2009). In Pinus taeda L., four cell wall genes, i.e., S-adenosyl methionine synthase 2, cinnamyl alcohol dehydrogenase, the water-stress inducible protein lp3-l, and α-tubulin, explained ∼3% of the variation in wood specific gravity, in the proportion of latewood to earlywood and in cellulose microfibril angle (González-Martínez et al. 2007). In P. radiata L., 10 single-nucleotide polymorphisms (SNPs) located on nine genes were associated with several traits (Dillon et al. 2010) but none of the genes overlapped with previous studies. This observation and the small number of genes tested to date suggest that these associations likely represent a small fraction of the genes that influence wood trait variation. Larger-scale analyses are needed to resolve this question and have become feasible with the availability of genome sequences for Populus and Eucalyptus (available at http://www.phytozome.net/), as well as large-scale cDNA resources in conifers such as Picea and Pinus (Mackay and Dean 2011).

The objective of the present study was to test for associations between SNP markers and wood quality and radial growth traits in white spruce [Picea glauca (Moench) Voss]. White spruce is widely distributed in North America (Nienstaedt and Teich 1972) and colonizes a wide range of sites covering most of Canada and a part of the United States (Nienstaedt and Zasada 1990). It exhibits high genetic variability for neutral markers (Jaramillo-Correa et al. 2001; Bousquet et al. 2007; Namroud et al. 2008), for growth and phenology (Li et al. 1993, 1997), and for wood traits (Corriveau et al. 1991; Beaulieu 2003; Lenz et al. 2010). As in many other tree species, spruce breeding aimed at wood quality has seen limited application because wood traits are time consuming and expensive to assess and change over the course of tree development. It is expected that DNA markers that explain a significant portion of the variation in wood traits would help to develop more rapid and effective breeding methods. Our specific objectives were (1) to identify single-marker associations with 25 wood traits using SNPs from several hundred candidate genes, (2) to determine the proportion of the variation in wood traits explained by these associations, (3) to test whether RNA transcript accumulation varies between the different genotypes for significant associations, and (4) to confirm the potential of a candidate gene approach for characterizing the genetic basis of wood traits.

MATERIALS AND METHODS

Association population:

The discovery population consisted of 492 30-year-old trees representing 165 open-pollinated families from 40 provenances sampled in Quebec, Canada. For each family, 3 trees were selected to cover as much as possible of the existing range of diameters at breast height (dbh) observed in the family. Needle tissue as well as 12-mm increment cores taken at breast height were collected, stored on ice, and transported to the Canadian Forest Service facilities where they were stored in a freezer at −10° until further treatment.

Phenotypic data:

Pith to bark profiles of wood physical attributes were obtained using the SilviScan technology (Evans and Ilic 2001) at FPInnovations (Vancouver, BC, Canada). SilviScan is designed to rapidly measure wood attributes using a combination of X-ray densitometry, X-ray diffractometry, and image analysis. Measurements were taken from the radial surface of 2 × 7-mm wood flitch samples prepared from the 12-mm cores. The thawed samples were acetone extracted and conditioned at 40% relative humidity. Air-dry (20°) density was measured in 25-μm steps using X-ray densitometry. Microfibril angle (MFA) was measured in 1-mm steps using X-ray diffractometry, and wood modulus of elasticity (MOE) was calculated from the densitometry and diffractometry data. Cell wall thickness as well as radial and tangential cell dimensions (∼1.5 μm) were obtained using optical microscopy and image analysis. From the raw data, early-, late-, and total wood averages were weighted by calculated annual ring area for each of the attributes (Table 1).

View this table:
TABLE 1

Eight SilviScan wood physical attributes measured for earlywood (EW), latewood (LW), and total wood (TW) as well as the percentage of earlywood in white spruce

Wood traits (supporting information, File S1) were submitted to principal component analysis (PCA) to create a small number of orthogonal composite wood characters that could explain a large proportion of the variation observed in the raw characters (Legendre and Legendre 1998). PCA was carried out using PROC PRINCOMP (Sas Institute 2008) with standardized data. Three components were found to be significant (Frontier 1976), and altogether they explained 83% of the total variation. The first principal component was mainly influenced by wood stiffness, density, and cell wall thickness. The second one was composed of coarseness and cell diameters in both the tangential and the radial directions. Microfibril angle was the trait having the largest influence on the third one.

Selection of candidate genes:

White spruce candidate gene sequences were identified from >16,500 white spruce consensus sequences that were assembled from 50,000 partially sequenced cDNA clones produced by the Arborea project (Pavy et al. 2005a). A total of 747 candidate genes were identified for SNP discovery and classified on the basis of expression data or information available for each gene. The expression data included an analysis of vascular tissues and needles in P. glauca with a custom cDNA array based on 9503 gene clusters (Pavy et al. 2008a) and three such studies of secondary xylem carried out in pine species (Plomion et al. 2000; Pavy et al. 2005b; Paiva 2006). The functional annotations considered white spruce transcription factors, as well as sequences known to be involved in cell wall and lignin synthesis, as well as sequences reported to be potentially involved in wood formation.

SNP discovery by resequencing genomic DNA:

PCR primers for amplification and resequencing were designed using Primer3 software (Rozen and Skaletsky 2000). Whenever possible, one of the primers was anchored outside of the coding regions to ensure specificity of amplification (Pavy et al. 2008b). Primer pairs were typically designed upon either one or two regions per gene.

Genomic DNA was extracted from 100 mg of fresh foliage of each of 25 individuals with the DNeasy Plant Mini Kit (QIAGEN, Mississauga, ON, Canada). DNA concentrations were assessed by GeneSpec spectrophotometer (MiraiBio, Alameda, CA). For the detection of polymorphism, DNA samples were pooled in equal amounts (4 ng of DNA per individual) before polymerase chain reaction (PCR) amplification following Pelgas et al. (2004). Genomic DNA was also extracted from ground haploid megagametophyte tissue to detect paralogs in SNP discovery (Pelgas et al. 2006; Pavy et al. 2008b). PCR amplifications were performed in volumes of 30 μl containing 10–20 ng template DNA, 20 mM Tris-HCL (pH 8.4), 1.5–2.5 mM MgCl2, 220 μM each dNTP, 0.25 μM each primer, and 0.04 unit Platinum Taq DNA polymerase (Invitrogen, Carlsbad, CA). PCR reactions were performed on a Peltier Thermal Cycler PTC200 (DNA Engine, DYAD; MJ Research, Waltham, MA) with an initial denaturation of 4 min at 94°, followed by 12 cycles of touchdown PCR (a 30-sec denaturation step at 94°, an annealing step with temperature decreasing from 66° to 60°, and a 3-min extension step at 72°) and by 35 cycles of a 30-sec denaturation step at 94°, then a 30-sec annealing step at 58°, and a 3-min extension at 72°, with a final extension of 10 min at 72°.

Single-locus products were sequenced for both DNA strands using BigDye Terminator v3.0 cycle sequencing reaction kits (Applied Biosystems, Foster City, CA) and an automated ABI Prism 3700 Genetic Sequencer (Applied Biosystems). DNA sequence strands were analyzed with SEQMAN software (Swindell and Plasterer 1997). SNPs were detected from the polymorphic positions indicated by double peaks or frameshifts in sequence chromatograms of the DNA pool (Pelgas et al. 2004). For each SNP detected, haploid DNA sequences from individual megagametophyte tissue were used to eliminate any heterozygous nucleotide position due to paralogy (Pelgas et al. 2006; Pavy et al. 2008b). A total of 809 of the SNPs discovered using genomic sequences were used to construct part of a 1536-SNP GoldenGate genotyping array (Illumina, San Diego).

In silico SNP discovery in expressed sequence tags:

In silico SNPs were detected using expressed sequence tags (ESTs) obtained from Sanger sequencing of 17 white spruce cDNA libraries (Pavy et al. 2005a, 2006). Essentially, the search for SNPs was made using ∼6500 contigs derived from at least two sequenced clones. SNP prediction was made using a proprietary method developed by Gydle (Quebec City, QC, Canada) or PolyBayes version 3.0 (Marth et al. 1999), with a p-prior of 0.01 and a probability that a position represents one SNP (PSNP) ≥ 0.99 (Pavy et al. 2006). A total of 727 SNPs among the in silico SNPs identified were selected and added to the 809 already discovered to complete the panel of 1536 SNPs.

SNP genotyping:

The minimum distance between any pair of SNPs selected for genotyping was 200 bp, whether considering genomic DNA or cDNA sequences. This distance was based on the information on low LD in Norway spruce, black spruce, and white spruce genes (Huertz et al. 2006; Namroud et al. 2010). About 20% of the SNPs that had been identified by resequencing were retained for genotyping.

Genotyping was conducted at the Genome Quebec Innovation Centre (McGill University, Montreal), using the highly multiplexed Illumina GoldenGate assay (Shen et al. 2005) with 250 ng of template DNA (at a rate of 50 ng/μl) for each of the 492 trees. Positive and negative controls were added to each 96-well sample plate. The intensity data for each SNP were normalized and assigned a cluster position (and resulting genotype) with the BeadStudio software (Illumina), and a quality score for each genotype was generated. A GenCall score cutoff of 0.25 was used to determine valid SNP genotypes and each retained SNP was required to have a minimum GenTrain score of 0.25 (Fan et al. 2003; Pavy et al. 2008b). Gentrain scores measure the reliability of SNP detection on the basis of the distribution of data points for each genotypic class.

Of the 1536 SNPs arrayed, 1021 (66.5%) were found to be polymorphic with a minimum GenTrain score of 0.25, which is similar to the success rate observed with a previous array in the same species (Pavy et al. 2008b). On the basis of positive controls, the repeatability of the assay for these SNPs was 99.98%. By selecting 472 trees that could be genotyped efficiently for at least 90% of their SNPs, the number of orthologous valid SNPs was reduced to 969. Moreover, rare SNPs as well as SNPs showing abnormal distributions in this population subset were excluded, such as those with very low (<1%) heterozygosity. Hence, the final data set consisted of 472 trees genotyped for 944 high-quality SNPs representing 549 different genes. Three of these 472 individuals could not be used in the association study because their wood traits could not be assessed (File S2).

Linkage disequilibrium:

Unphased linkage disequilibrium between pairs of SNPs within each candidate gene was evaluated using the squared allelic correlation coefficient (r2). The decay of LD with distance of base pairs between sites within the same candidate gene was estimated with both a linear model and a nonlinear regression, using the nls function in the nmle package of the R software v2.12.2 (Hill and Weir 1988; R Development Core Team 2010).

Population structure:

Associations between markers and traits are susceptible to errors and biases caused by covariance between the marker and polygenic effects (Kennedy et al. 1992). Although population structure is rare in outbred forest species (Neale and Savolainen 2004), it is considered the most serious systematic bias producing false-positive associations (Marchini et al. 2004). When population structure is present, bias occurs because subgroups of relatives tend to share more markers and gene alleles genome-wide (Breseghello and Sorrells 2006) than a pair of individuals drawn at random from the population.

White spruce has previously been shown to harbor little or no population genetic structure for neutral markers in the region sampled (Jaramillo-Correa et al. 2001; Namroud et al. 2008, 2010). The presence of population structure was further checked for our data set using the model-based clustering algorithm implemented in the software STRUCTURE, version 2.2 (Pritchard et al. 2000, 2007). The entire population of 469 trees was analyzed, assuming that the 944 SNP loci were unlinked. Contrary to previous studies (see, for instance, González-Martínez et al. 2007), the discovery population was not genotyped using another class of neutral genetic markers (e.g., microsatellites) to test for population structure, as the large number of available SNPs was considered sufficient. The correlated allele frequency model was used as it is likely the most appropriate as white spruce populations are almost always expected to be connected by gene flow. This model implies that allele frequencies in each putative cluster are likely to be similar (Falush et al. 2003). Preliminary tests showed that STRUCTURE was detecting substructures partially confounded with the half-sib structure. Furthermore, as K values were increasing, log probability of data Pr(X |K) (Equation 12 in Pritchard et al. 2000) never reached a maximum value, which is needed to determine the uppermost level of structure. To avoid the detection of substructure confounded with half-sib structure, we randomly selected one individual in each of the 165 open-pollinated families making up the discovery population and analyzed this subset with STRUCTURE to determine the uppermost hierarchical level of structure in our discovery population. Two other subsets containing one tree per open-pollinated family were also composed, using the remaining trees. Then, 20 replicate runs of STRUCTURE analyses at the uppermost level of structure found with the first subset were performed for each of the three subsets, and the software CLUMPP (Jakobsson and Rosenberg 2007) was used to rank the clusters and estimate Q coefficients’ averages over replicate runs.

Test statistics for association:

Single-marker–based tests were performed for each of the wood traits assessed for earlywood, latewood, and total wood, as well as for the three principal components retained. Moreover, two genetic models were tested for each trait, i.e., the additive (three genotype classes) and dominant (two genotype classes) effects models. A mixed linear model (MLM) was fit for each pair of phenotypic traits and SNP markers (Yu et al. 2006; Malosetti et al. 2007). This model accounts for relatedness among individuals through a pairwise kinship matrix estimated from marker genotypes. As kinship coefficients do not always reflect with accuracy the known relatedness among individuals, such as the half-sib structure present in our discovery population, we first built a matrix with known half-sib relationships and used it in a first series of analyses (File S3). We also built a relatedness matrix with Nason's kinship coefficients (Loiselle et al. 1995) to check whether the SNPs found to be significantly associated with wood traits using the half-sib structure of relatedness would still be significant when using Nason's kinship coefficients (File S4). Estimates were obtained using the 944 SNPs for the 469 trees with the software SPAGeDi, version 1.2 (Hardy and Vekemans 2002). As suggested by Yu et al. (2006), negative kinship values were set to zero and a 469-square matrix was fit using the SAS MIXED procedure (Littell et al. 2006; Sas Institute 2008) in the following “K” (per notation in Yu et al. 2006) MLM,y=Sα+Qν+Zu+e,(1)where y is a vector of preadjusted wood quality measurements, α is a vector of fixed SNP effects, u is the vector of random breeding values or “background” polygenic effects, e is the vector of random residual error effects, and S, Q, and Z are incidence matrices of 0’s and 1’s relating records (y) to fixed (α and ν) and random (u) effects. Expectation of the model is E(u) = E(e) = 0 and thus E(y) = Sα + , where random additive and residual variances are assumed to be Var(u)=G=2Kσa2 and Var(e)=R=Iσe2, but are estimated with REML and K and I equal to kinship and identity coefficient matrices, respectively.

Failure to appropriately adjust for multiple testing may produce excessive false positives or overlook true positive signals in association studies when using large numbers of SNP markers. To correct for multiple testing (944 analyses per wood trait), the positive false discovery rate (FDR) method (Storey and Tibshirani 2003) was used to identify significant SNPs after correcting for multiple testing using the QVALUE software, version 1.0 (Storey 2002; Storey et al. 2004). Normality of studentized residuals was checked using the Kolmogorov D and the Shapiro–Wilk W statistics obtained with the “UNIVARIATE” procedure in SAS (Sas Institute 2008). Phenotypic data were not transformed as they were normally, or approximately normally, distributed.

Microarray transcript profiling:

Given current limitations and costs of high-throughput genotyping technologies, a candidate gene approach should represent a suitable strategy for association studies with complex traits in species haboring large genomes and little LD. One criterion to consider for selecting genes is the preferential accumulation of transcripts in tissues directly linked with the traits of interest but this had not been directly tested in forest trees. This report aimed to test whether the 549 candidate genes tested in the association study were preferentially expressed in wood tissues and whether the candidate genes harboring significant SNPs were more likely to show a tissue-specific expression pattern than genes not harboring significant SNPs. A transcript profiling experiment compared secondary xylem (Xy), secondary phloem (Ph), and needles (Ne) using a large-scale custom microarray that included the 549 candidate genes.

Tissue samples were obtained from 3-year-old in vitro clonally propagated white spruce seedlings of one genotype (see Pavy et al. 2008a for details). All the tissues were frozen in liquid nitrogen. Total RNA was isolated following the method of Chang et al. (1993) with modifications described in Pavy et al. (2008a). The custom microarray was composed of a total of 32,640 oligonucleotides, 25,094 of which represented unique P. glauca gene sequences. A description of the sequence data is available at http://www.arborea.ca. Microarray hybridizations used five independent samples of each tissue type; two samples labeled with different dyes were hybridized to each slide (dye swaps) and image analyses proceeded following Pavy et al. (2008a). Data analysis was carried out using BIOCONDUCTOR (http://www.bioconductor.org) in the R statistical environment and following the method described in Pavy et al. (2008a). Briefly, data were normalized and the fold changes in spot intensity were estimated by fitting a linear model, using the limma package (release 2.0.7, http://www.bioconductor.org/packages/bioc/1.8/vignettes/limma/inst/doc/usersguide.pdf). Contrasts between the three types of tissues were estimated for each of the 549 genes, and significance of the differential expression was tested with a moderated t-statistic with P-values adjusted for multiple testing with the method of Benjamini and Hochberg (1995). The candidate genes were classified according to their expression profiles and placed into 1 of 10 tissue preference classes (see Figure 1). A likelihood-ratio test (Sokal and Rohlf 1981) allowed testing for the independence of the tissue preference class distributions between the genes that carried SNPs that were associated with each of the wood traits at α ≤ 0.05 before correction for multiple testing, as well as the genes that had nonsignificant SNPs.

Figure 1.—

Distribution of genes into tissue-preferential transcript accumulation classes. (A–F) Each chart represents a different set of gene sequences (Xy, xylem; Ph, phloem; Ne, needles). (A) Whole-trancriptome microarray (25,094 genes). (B) Genes containing a SNP significantly associated with wood property traits (13 genes). (C) Genes containing a SNP submitted to association analysis (549 genes). (D–F) Genes containing a SNP significantly associated with wood traits at P ≤ 0.05: (D) MFA (60 genes), (E) tracheid cell diameter in radial direction (RCD) (57 genes), and (F) DEN (48 genes).

Transcript analysis of SNP genotypes:

RNA transcript levels were determined for genes underlying the SNPs significantly associated with wood traits to test whether RNA transcript accumulation varied between the different genotypic classes of the SNPs. Because of the large number of determinations, only genes with significant SNPs after correction for multiple testing were targeted (see RESULTS). Transcript levels were determined by RT–qPCR with gene-specific primers. The latter could not be designed for either of the Tubulin genes because of the high level of sequence similarity among family members. For each SNP genotypic class, 8–25 trees were individually sampled by removing three small plugs of ∼1 cm2 each of actively growing secondary xylem tissue at 1.3 m from the ground; these represented the entire annual growth produced near the end of the earlywood phase. Tissue handling and RNA extractions were performed as described for white spruce (Bedon et al. 2007); complementary DNA synthesis was prepared from 500 ng of total RNA (Boyle et al. 2009). PCR mixtures were assembled using LightCycler 480 SYBR Green I Master (Roche, Basel, Switzerland) with an epMotion 5075 pipetting robot (Eppendorf, Hamburg, Germany) and amplifications used a LightCycler 480 (Roche) with conditions and melting curve analyses as described (Boyle et al. 2009). The number of RNA transcript molecules was determined using the LRE methodology (Rutledge and Stewart 2008) adapted for Excel (Boyle et al. 2009). Normalization was performed using GeNorm (Vandesompele et al. 2002) and the reference genes ef1a (BT102965), cdc2 (BT106071), and rpl3A (BT115036). The measured number of transcript molecules was transformed to a log scale, and differential expression across three or two genotypic classes was tested by ANOVA or t-tests, respectively.

RESULTS AND DISCUSSION

Linkage disequilibrium:

The pattern of the squared allelic correlation coefficient (r2) with base pair distance within candidate genes (Figure 2) illustrated the rapid LD decay in our white spruce population, with r2 values dropping to <0.25 within <100 bp. These results are congruent with those of Namroud et al. (2010) from complete or nearly complete gene sequences in the same species.

Figure 2.—

Patterns of intralocus linkage disequilibrium (LD) for white spruce estimated using the squared allelic correlation coefficient between each pair of SNPs within candidate genes.

Population structure:

A weak but nonsignificant population structure with the presence of two subpopulations was detected using the model-based clustering method, as the patterns obtained were as expected for structured populations (Pritchard et al. 2007); i.e., the log probability of Pr(X |K) values reached a maximum for K = 2 (Figures 3 and 4). The presence of a weak population structure or even its absence in white spruce at the regional level was expected on the basis of results of studies using various classes of neutral genetic markers (Alden and Loopstra 1987; Furnier and Stine 1995; Jaramillo-Correa et al. 2001; Namroud et al. 2008, 2010). Despite the nonsignificance of this weak population structure, it was nevertheless considered in the mixed model analyses to minimize the possibility of any false positive association due to the weak population structure signal.

Figure 3.—

Plot of the mean of L(K) over 20 runs for each of K clusters ranging from 1 to 7 as a function of the number of clusters. A constant of 15,500 was removed on the y-axis.

Figure 4.—

Genetic clustering of the 469 white spruce individuals of the discovery population using the program STRUCTURE. Vertical bars represent individuals, and shaded and solid areas represent proportional membership of each individual in each of the two clusters. Individuals are ranked by latitude and longitude.

Genetic association:

Results of single-marker associations with each of the 25 wood traits as well as with each of the significant principal components are presented in File S5 and File S6. For each of the three principal components, we found between 11 and 15 significant associations (P ≤ 0.01) with SNPs when using the additive effects model and between 10 and 14 SNPs for the dominant effects model. However, after correction for multiple testing with a significance level of Q ≤ 0.10, none of these associations were still significant. For the individual wood traits, a total of 23,600 (944 SNPs × 25 traits) and 22,400 (896 SNPs × 25 traits) association tests were performed for the additive and dominant effects models, respectively. Globally, 6.4% and 5.8% of these tests were significant at a level of P ≤ 0.05, respectively, and 405 (43%) SNPs and 425 (47%) SNPs were found associated with at least one of the wood traits. These proportions are of the same order of magnitude as those found for cold hardiness-related traits in Douglas fir [Pseudotsuga menziesii var. menziesii (Mirb.) Franco], for which 7.8% of the association tests performed were found to be significant at this nominal threshold level (Eckert et al. 2009). They are also close to the 8.3% of significant associations (106 of 1276 tests) found by González-Martínez et al. (2007) for wood traits in P. taeda.

The number of significant associations was slightly reduced when using a nominal threshold of P ≤ 0.01. Hence, for each of the wood traits assessed, we found 8–19 significant SNPs for earlywood, 11–18 SNPs for latewood, and 7–18 SNPs for total wood, depending on the genetic effects model (Figure 5); this represents an overall 1.5% of significant association tests and the numbers of SNPs associated with at least one wood trait were 136 (14%) and 148 (17%), for the additive and dominant effects models, respectively. The phenotypic variation explained by each of these significant SNPs varied between 1.5% and 5.4%. The total variation in wood traits that could be explained by considering all these SNPs simultaneously if we were assuming an additive effects model and the absence of interaction between significant SNPs could sum up to 60%. However, the exact percentage could be obtained only when considering simultaneously all these SNPs in a MLM analysis.

Figure 5.—

Number of significant (P ≤ 0.01) SNP associations before correction for multiple testing for eight early- (EW), late- (LW), and total wood (TW) traits (see Table 1) in white spruce. Associations were tested using an additive (ADD) or dominant (DOM) effects mixed linear model.

After correction for multiple testing, the number of significant associations was reduced to 25, with 8, 4, and 3 significant associations for early-, late-, and total wood traits, respectively, with an additive effects model, and with 4, 4, and 2 associations, respectively, with a dominant effects model (Table 2). Some SNPs were found to be significantly associated with more than one trait whereas others were positive with both the additive and the dominant effects models. Thus, globally, 13 different SNPs (discovery rate of 1.4%) harbored by as many genes (2.4% of the sampled gene set) were found to be significantly associated with 10 different characters. Among these 13 SNPs, 5 represented synonymous substitutions, 3 were nonsynonymous, and others were located in untranslated regions. Silent SNPs should not be considered a priori as potential false positives because they could affect transcript levels and codon usage (e.g., Kimchi-Sarfati et al. 2007; Chamary and Hurst 2009). While any of these SNPs could be quantitative trait nucleotides (QTNs), they could also be linked to adjacent causal but untested SNPs.

View this table:
TABLE 2

Significant associations between gene SNPs and wood traits in white spruce after correction for multiple testing

The number of significant SNPs is consistent with previous reports of association studies in forest trees. For example, González-Martínez et al. (2007) reported 7 significant associations with wood traits in P. taeda at a FDR of Q ≤ 0.10. Recently, using the same statistical threshold, Dillon et al. (2010) reported 10 SNPs significant with a variety of wood traits in P. radiata. Thumma et al. (2005) reported 2 SNPs associated with microfibril angle in E. nitens. For climate-related traits such as cold hardiness and phenological characters, there seems to be a trend toward larger numbers of significant associations. Besides Ingvarsson et al. (2008) who found 2 SNPs associated with bud set in Populus tremula L., Eckert et al. (2009) reported 30 significant associations related to adaptation in Pseudotsuga menziesii var. menziesii, whereas Holliday et al. (2010) detected 35 significant SNPs for bud set and frost hardiness in P. sitchensis (Bong.) Carr.

The phenotypic variation in wood traits captured by SNPs deemed significant by FDR in our study varied between 2.6 and 5.4% (Table 2). Small SNP effects are consistent with other association studies in which between 1.5 and 6.5% of the total phenotypic variation was accounted for by SNPs (Thumma et al. 2005; González-Martínez et al. 2007; Ingvarsson et al. 2008; Dillon et al. 2010). These small SNP effects are in accordance with polygenic quantitative models of wood traits (Zobel and Jett 1995; Brown et al. 2003; Pot et al. 2006).

When analyzed simultaneously and without taking marker interactions into consideration, the cumulative effect of multiple significant SNPs for the same trait could explain a higher proportion of the total phenotypic variation in some cases (e.g., 9.7%, 8.7%, and 11.1% for average ring width in latewood, cell wall thickness in earlywood, and percentage of earlywood, respectively, with the additive effects model; Table 2). For some of the significant SNPs in the present discovery population, one of the homozygous genotype classes had very few observations (<1%) in both the additive and the dominant effects models (Table 2). In these cases, the corresponding genotypic effects could be biased and partly drive the significance of SNP loci under the dominant effects model. However, when examining phenotypic data, we observed that these individuals possessing rare homozygous genotypes did not necessarily have the most extreme phenotypes. Thus, the effects of these loci may in fact be real. To confirm this, it would require (1) a larger discovery population to capture a greater frequency of low-frequency genotypes and/or (2) an analytical method that accounts for infrequent observations resulting in less bias.

Genes carrying significant SNPs belonged to 11 known gene families and 1 family of unknown function. These gene families are listed in Table 2. The fact that significant associations were found with both growth (e.g., average ring width, percentage of earlywood) and wood quality (e.g., tracheid coarseness, microfibril angle) traits is encouraging, as we are usually interested in selecting for both traits simultaneously. In this study, we could not find any SNP significantly associated with wood density after correction for multiple testing. Wood density is determined by several wood characteristics such as cell size and wall thickness, the ratio of earlywood to latewood, and other factors (Zobel and Van Buijtenen 1989). This is likely the main reason why no significant associations were observed between wood density and SNPs once a correction for multiple testing was performed. This trend might also be specific to the present study because significant associations with wood density were reported for P. taeda (González-Martínez et al. 2007).

Tissue profiling and genotypic transcript accumulation:

Transcript accumulation levels were determined by using a large-scale custom oligonucleotide microarray that was composed of 25,094 probes, each estimated to represent a different white spruce gene (see materials and methods). Total RNA hybridizations against these probes were for comparisons between secondary xylem, secondary phloem, and young needles (Figure 1). Xylem-preferential RNA accumulation compared to that in needles or phloem was found for 10 of the 13 genes harboring SNPs significantly associated with one or more of the wood traits, and the majority of genes (7 of 13) were xylem preferential compared to both needles and phloem (Table 3 and Figure 1B). Transcripts of the expansin sequence (GQ0172_O22), a tubulin (GQ03005_C12), and a sequence of unknown function (GQ03006_P17) were the most xylem-preferential sequences. Two of the genes had low transcript levels in all tissues, and only the receptor kinase-like transcripts accumulated to lower levels in xylem compared to either needles or phloem.

View this table:
TABLE 3

Tissue and organ preferential transcript accumulation of gene sequences carrying the gene SNPs significantly associated with wood traits

These transcript accumulation profiles were compared to microarray data for the 549 genes tested for association with wood traits and for the transcriptome represented by the microarray (Figure 1). A total of 4133 (16% of the total) of the sequences on the microarray accumulated preferentially in xylem compared to both phloem and needles, and 672 (3%) sequences were xylem preferential in reference to phloem or needles (Figure 1A). Of the 549 genes tested for association, 29% accumulated preferentially in xylem compared to both phloem and needles (Figure 1C), as was expected considering the criteria for selecting the candidate genes. Many of these 549 genes were also preferential to phloem (29%) and only 7% were preferential to needles (compared to 20% in the whole transcriptome). Considering the genes and SNPs putatively associated with wood traits before correction for multiple testing, the distributions of transcript profiles varied slightly from trait to trait but they were not significantly different from the overall set of 549 candidate genes tested (not shown). For example, the proportion of xylem-preferential transcripts varied from 36% for cellulose MFA (Figure 1D) to 24% for tracheid cell diameter (TCD) (Figure 1E), while density (DEN) was slightly enriched for phloem preferential transcripts (Figure 1F).

The association study reported here considered many more candidate genes than previous investigations related to conifer wood properties, which targeted small sets of genes known to be expressed in secondary xylem tissues (e.g., González-Martínez et al. 2007; Dillon et al. 2010). Our analysis selected xylem-expressed transcripts and included sequences that were not preferentially expressed in xylem. For example, the 24 MYB transcription factor sequences and the 13 glycosyl hydrolase protein sequences gave diverse expression profiles (not shown). The observation that most of the SNPs significantly associated with wood property traits were in genes with xylem-preferential transcript accumulation indicates that expression data are relevant for selecting candidate genes. Interestingly, a few genes with low transcript levels in xylem were also associated with phenotypic variation in wood traits, but one of these genes was found to be strongly expressed in differentiating xylem of mature trees (Table 4), whereas tissue preferential patterns were determined in whole xylem of young trees (Figure 1).

View this table:
TABLE 4

Transcript accumulation comparing genotypic classes for selected gene SNPs significantly associated with wood traits

Transcript levels could be compared among the different genotypic classes for six of the significantly associated SNPs (after FDR correction for multiple testing) by using RT–qPCR with gene-specific primers as described in Boyle et al. (2009). The assays used differentiating secondary xylem from the 30-year-old trees of the association population. Transcript accumulation varied significantly between genotypic classes of the SNPs for the LACS9, β-TIP, and ACC oxidase genes (Table 4), which represents 50% of the SNPs analyzed. The MYB and the expansin transcripts did not differ significantly between genotypes, and in determinations of the glycoside hydrolase family 28 sequence transcripts were too variable among trees to draw any conclusions. These differences in transcript accumulation among genotypic classes are discussed below with their putative function and relationship with variation in phenotypic trait. To our knowledge, only Thumma et al. (2009) previously reported a link between SNPs, transcript levels, and phenotypic variation in a forest tree population. These authors identified one SNP in a COBRA-like gene of E. nitens (EniCOBL4A) that was associated with allelic expression imbalance and with the amount of cellulose deposited in the wood.

Potential role of genes containing significant association SNPs:

Long-chain acyl-CoA synthetase 9 (LACS9):

The LACS enzymes provide the acyl-CoA pools used in the synthesis of phospholipids and triacylglycerols in vegetative and seed tissues (Pongdontri and Hills 2001; Shockey et al. 2002). They play essential roles in development via fatty acid metabolism, including the provision of energy. The nine LACS characterized in Arabidopsis were shown to activate 14–20 carbon fatty acids (Shockey et al. 2002) and are expressed in all tissues. Extended analyses of the LACS gene family have begun to reveal the roles of individual members in the peroxisome (Fulda et al. 1997) and chloroplast (Schnurr et al. 2002). Interestingly, the LACS9 transcript levels were low in the young trees used for the tissue profiling (not shown), but were among the highest in the secondary xylem of the mature trees (Table 4).

Individuals of only one of the two homozygous classes were present in the discovery population and, on average, they showed larger earlywood and total wood average ring width than heterozygotes (Table 2 and Figure 6a). The SNP (PGWD1-1094) was located in the 3′-UTR, which could indicate a role in the regulation of transcription, RNA stability, or post-transcriptional RNA processing. The transcript levels differed significantly between genotypes (P = 0.02) (Table 4). The estimated allelic effects of the SNP (PWD1-1094) on the wood trait (Figure 6a) corresponded well with estimates of transcript levels (Figure 6b); i.e., heterozygotes had both lower average width of the latewood ring (ARW) and lower transcript levels than the homozygotes.

Figure 6.—

Relationship between SNP genotypic classes and transcript accumulation. (A) Estimated genotypic effects of a significant SNP (PGWD1-1094) for earlywood average ring width (and TW-ARW, see Table 2) in the additive effects model. (B) Transcript levels of the corresponding gene determined by qPCR.

Glycoside hydrolase family 28 [polygalacturonase (pectinase) family protein] (GQ02829_F04):

The glycoside hydrolase protein family 28 (GH28) is a large family of plant proteins that comprises enzymes with several known activities. Pectin-modifying enzymes have well-established roles in the primary cell wall, in particular during cellular growth and in biotic interactions; however, data are lacking about wood formation. Seven putative family members were represented on the genotyping chip and all but one of them accumulated transcripts preferentially in phloem (four) or xylem (two). The SNP (PGWD1-1396) is a synonymous substitution (TGC–TGT) coding for a cysteine. It was represented by only one class of homozygous individuals in the discovery population. On average, the homozygotes had thinner cell walls than the heterozygotes (1.77 μm vs. 1.86 μm); the transcript levels did not differ between genotypes.

Glycoside hydrolase family 10:

The glycoside hydrolase protein family 10 (GH10), formerly known as cellulase family F, comprises several enzymes with known activities. They play important roles in primary cell wall metabolism, particularly in the reassembly and degradation process (Girke et al. 2004) and in secondary cell wall metabolism (Ichinose et al. 2010). The Populus genome encodes seven GH10 xylanases, but only PttXyn10A had detectable ESTs in the wood-forming tissue (Geisler-Lee et al. 2006) and was highly upregulated in different stages of wood formation (Aspeborg et al. 2005; Takahashi-Schmidt 2008). Only one family member was included in our genotyping chip. The SNP (PGWD1-0107) consisted of a nonsynonymous substitution (CCA–GCA) coding for proline and alanine that was found to be associated with latewood MFA in a dominance effects model.

Pectinesterase:

Pectin methylesterases (PME) are known to catalyze the demethylesterification of cell wall polygalacturonans, one of the major polysaccharides in the primary wall and middle lamella. In dicotyledonous plants, these ubiquitous cell wall enzymes are involved in important developmental processes including cellular adhesion and stem elongation (Micheli 2001). It has been suggested that the prevalence of neutral PMEs in active cells would allow cell wall expansion (Guglielmino et al. 1996). Pectinesterase transcripts were among the differentially regulated and most strongly expressed genes in tension wood of Populus (Andersson-Gunneras et al. 2006) and in the juvenile cambial region of E. grandis (Gallo De Carvalho et al. 2008). A pectin methyl esterase from wood-forming tissues of hybrid aspen (P. tremula x tremuloides Michx.) was found to act as a negative regulator of both symplastic and intrusive growth of developing wood cells (Siedlecka et al. 2008). Fourteen different pectinesterase sequences were present on our genotyping chip and all but 2 of them accumulated more transcripts in phloem (7) or xylem (5). The SNP (PGWD1-1096) represents a nonsynonymous substitution (TGC–GGC) coding for cysteine and glycine that was associated with both earlywood and total wood cell wall thickness.

MYB4 (R2R3-MYB domain DNA-binding protein):

The MYB4 gene has been linked to secondary cell wall formation, including lignification in the secondary xylem in both P. taeda and P. glauca. In P. glauca, transcripts levels were highest in the secondary xylem of the main stem and of large roots and were upregulated during the formation of compression wood (Bedon et al. 2007). Closely related MYBs of P. taeda (Patzlaff et al. 2003a,b) and of E. gunnii Hook f. (Goicoechea et al. 2005) have been shown to activate the transcription of genes encoding lignin synthesis enzymes and to regulate lignin accumulation and cell wall formation when overexpressed in angiosperm model plants. Overexpression of the P. taeda MYB8 gene in transgenic spruce resulted in ectopic secondary cell wall thickening (Bomal et al. 2008). The SNP (PGWD1-1282) represents a synonymous substitution (GCA–GCG) coding for alanine. It was associated with cell wall thickness and fiber coarseness, both of which are related to wood density and radial and tangential cell diameter. Heterozygous and dominant homozygous trees showed significantly higher fiber coarseness for the whole wood ring (279 μg/m) and for earlywood (269 μg/m) than the alternate homozygous type (266 μg/m and 255 μg/m, respectively, Table 2). A similar trend was observed for cell wall thickness (Table 2). The transcript levels did not differ significantly (P = 0.16) between genotypes, despite the large number of individuals tested (Table 4).

Galactosyl-transferase XT2 (UDP-xylosyltransferase 2):

The hemicellulose xyloglucan (XyG) is believed to play an important role in cell wall structure and function (Eckardt 2008). A number of enzymes participate in XyG biosynthesis, including b-(1/4)-glucan synthase, a-fucosyltransferases, b-galactosyltransferases, and α-xylosyltransferases (Girke et al. 2004). A lack of detectable xyloglucan in the xxt1/xxt2 double mutant in Arabidopsis caused a significant reduction in the stiffness and strength parameters of the mutants (Cavalier et al. 2008). The xyloglucan acts as a tether in the extracellular matrix and plays a key role in the loosening and tightening of cellulose microfibrils (Hayashi and Kaida 2011). The SNP (PGWD1-0511) represents a synonymous substitution (TCA–TCG) coding for a serine. It was associated with earlywood percentage.

Receptor-like kinase (ATP binding/kinase/protein serine/threonine kinase) (GQ03211_O01):

Receptor-like kinases (RLKs) control a wide range of processes, including development, disease resistance, hormone perception, and self-incompatibility. The Arabidopsis genome contains >600 RLK homologs, representing nearly 2.5% of the annotated protein-coding genes (Shiu and Bleecker 2001). In Arabidopsis, the PXY receptor-like kinase plays a role in the maintaining of cell polarity in the meristem and is required for the differentiation of specialized and spatially separated xylem and phloem cells (Fisher and Turner 2007). In poplar, transcripts of three receptor-like kinases accumulate a different stage of secondary xylem differentiation (Schrader et al. 2004). In P. glauca, transcript accumulation was clearly preferential to phloem, suggesting that communication between xylem and phloem may affect wood properties. The SNP (PGWD1-1313) is a synonymous substitution (ACA–ACG) coding for threonine in the exon. It was associated with both the ARW and the proportion of earlywood (PCT), two closely related traits.

ACC oxidase (1-aminocyclopropane-1-carboxylate oxidase):

This enzyme forms ethylene from its precursor, ACC. In plants, ethylene has signaling functions such as stimulation of fruit ripening and is required for CO2 activity. The closest ACC oxidase homolog in Arabidopsis was regulated during in vitro xylem vessel element differentiation (Kubo et al. 2005). Endogenous ACC was shown to accumulate in the vascular cambium specifically on the lower side of branches and was linked to compression wood differentiation in P. contorta Dougl. ex Loud. ssp. latifolia (Savidge et al. 1983). Evidence also suggested that ethylene plays a role in the control of xylem differentiation by inducing the activity of enzymes involved in lignification (Miller et al. 1984). In P. glauca, ACC oxidase transcripts were moderately preferential to secondary xylem (Table 3). ACC oxidase protein accumulated significantly more in secondary xylem compared to needles and in response to gravitational stress in P. pinaster Ait. (Plomion et al. 2000). The SNP (PGWD1-0560) represents a nonsynonymous substitution (GGG–GAG) coding for glycine and glutamine. It was also associated with both latewood average ring width and earlywood percentage. The transcripts varied between the two genotypic classes that could be compared (P = 0.03) (Table 4).

β-Expansin [AtEXPB3 (Arabidobsis thaliana Expansin B3)]:

Expansins are a large family of proteins that play major roles in plant growth and development by acting on cell walls. β-Expansins act through cell wall loosening in a variety of developmental functions. They have been mostly described in grasses; no studies on gene expression or function in secondary cell wall development or wood formation have been reported. The SNP was located in an intron (A*–C*) and was associated with the PCT. The homozygous trees had significantly more earlywood (84%) than the two other genotypes (81%). Our data suggested that the transcripts may accumulate to lower levels in the homozygous genotype with higher PCT, but the difference was not statistically significant. The β-expansin gene might be somehow related to the differences in stiffness and cell wall thickness between earlywood and latewood, but this remains untested.

Tubulins (TUA2 and TUB3):

Several lines of evidence suggest that the xylem predominant TUA and TUB isoforms are specifically associated with cellulose synthesis during secondary cell wall formation (Oakley et al. 2007). Both of the P. glauca tubulins were expressed preferentially in xylem and shared sequenced homologies with tubulins from P. taeda that are highly abundant among ESTs from xylem cDNA libraries (Yang et al. 2004) and SAGE sequencing data (Lorenz and Dean 2002).

The P. glauca α-tubulin gene (TUA2) is highly similar to sequences that have been linked to the formation of cortical microtubules in P. taeda (González-Martínez et al. 2007), developing xylem-preferential transcript accumulation in P. radiata (Cato et al. 2006), ESTs that are unique to xylem librairies in Populus (Pilate et al. 2004), and transcripts levels that correlated with auxin concentration in developing secondary xylem in P. tremula x tremuloides (Nilsson et al. 2008). A β-tubulin gene was shown to be preferentially expressed in xylem in Eucalyptus (Paux et al. 2004). β-Tubulin (EgrTUB1) was differentially regulated in Eucalyptus branches and was reported to be likely involved in cell wall biosynthesis (Qiu et al. 2008). The TUA2 SNP (PGWD1-0863) represents a synonymous substitution (GGA–GGC) coding for glycine in a tubulin, whereas the TUB3 SNP (PGWD1-0354) is located in the 3′-UTR (T–C). Both of the SNPs were associated with latewood ARW; the SNP in TUB3 was also associated with earlywood percentage (PCT) (Table 2).

β-TIP (β-TONOPLAST INTRINSIC PROTEIN):

Membrane intrinsic proteins form a family of 35 genes in the Arabidopsis genome, and 13 of them encode tonoplast intrinsic proteins (TIPs) that localize in the vacuolar membrane (tonoplast) (Johanson et al. 2001). Most TIPs are classified as probable or confirmed aquaporins. In Arabidopsis, both the β-TIP (also known as TIP3,2) and the highly similar α-TIP (TIP3,1) are expressed preferentially in germinating seeds, while other classes of TIPs are strongly expressed during seed maturation and dehydration (Vander Willigen et al. 2006). In P. taeda, ESTs of a TIP and an aquaporin were overrepresented in a xylem cDNA library made from compression wood compared to normal wood (Whetten et al. 2001). Compression wood produces tracheids with altered cell shape and thicker and more lignified secondary cell walls. We found that β-TIP accumulated preferentially in xylem and phloem compared to needles. The SNP (PGWD1-1070) was located in an intron (T*–C*) and was associated with the tracheid diameter in the latewood (dominance model). The homozygote (homo2) had an average diameter of 23.6 μm, whereas the other genotypic class (heterozygotes and homo1) had a slightly lower average diameter of 23 μm (Table 2). The transcript levels also differed significantly between genotypic classes but did not appear to be positively correlated with the wood phenotype. The transition from earlywood to latewood is strongly influenced by water availability, and osmotic stress-related transcripts including aquaporins have been found to accumulate in latewood (Egertsdotter et al. 2004). The final latewood morphology appears to involve β-TIP but its functional role remains to be tested.

Methodological and statistical challenges:

Association studies in conifer populations have many attractive features because discovery populations with weak or nonexistent population structure and high nucleotide diversity can be assembled (Neale and Savolainen 2004); however, several challenges remain to be overcome. This study shows that a larger number of less frequent genotypes should be captured by increasing the sample size to better estimate their effects. A larger sample size may also offer the opportunity to include experimental design features in the model (e.g., blocks) so that estimates of genetic effects are adjusted accordingly. Second, additional information on marker covariance structure can be accommodated in the model only if significant LD is present among markers and QTL. In natural conifer populations, LD is expected to sink to very low levels over distances that are spanned of average-sized genes (e.g., r20.2∼ 100 bp in our data set), making the discovery of additional SNPs unlikely to yield strong LD across candidate genes. More recently, LD has been shown to also have substantial gene-to-gene variation, in natural populations of spruces (Namroud et al. 2010). As more SNPs are discovered, LD information within genes may be used to infer haplotypes and predict their effects, which has been found to yield more accurate estimates of genetic effects (Calus et al. 2008). Although LD falls to low levels in natural populations of conifers, smaller populations contain greater potential LD partly as a function of population size (Neale and Savolainen 2004). Thus, if genomic resources are fixed, working with smaller population sizes, such as elite breeding material, would be another way of utilizing LD information in an association analysis. In such situations, genomic selection (Wong and Bernardo 2008; Grattapaglia and Resende 2011) or multiple-gene selection could also become attractive.

Finally, there is a wealth of information on longitudinal phenotypic data that at first appears underutilized by taking overall means or by discretizing the data into cross-sectional “age” traits. From one standpoint, markers found to be significant over developmental ages, year environments, and locations can be interpreted as stable and thus possibly useful for MAS (González-Martínez et al. 2007). However, if trait (co)variances change over time, more complex, functional models may better serve both researchers whose interests lie in the prediction of genetic effects (e.g., breeders) and those interested in elucidating gene functions (e.g., geneticists and molecular biologists). In this study, the application of more complex functional models, which require more parameters to be estimated, was limited by the relatively small size of the data set. A forthcoming study is planned in an expanded sampling of our discovery population with individuals genotyped at more SNPs per gene and more candidate genes, taking these considerations into account.

Conclusion:

Given the large number of SNPs that explained significant amounts of phenotypic variation before corrections for multiple testing, it stands to reason that several of these SNPs could indeed not be false positives and be declared significant if the size of the discovery population was expanded. This is a line of research that we are pursuing actively. This study represents one of the first well-documented cases showing significant differential transcript accumulation among the genotypes of a SNP associated with a trait of interest in a population of conifer trees. It also shows that transcript accumulation profiles are appropriate for selecting candidate genes and finding SNPs associated with wood properties, because most of the significant SNPs were found in genes that had higher transcripts levels in secondary xylem. While the number of significant associations and the proportion of variation explained were small, they were in line with expectations derived from statistical and quantitative genetics theory. These results obtained using large-scale gene sampling and high-throughput genotyping approaches combined with gene expression represent significant progress toward the identification of genes that control wood traits and the development of successful marker-aided selection in conifer trees. This trend will continue given the current rapid advances in high-throughput sequencing and genotyping.

Acknowledgments

This study is part of the Arborea project and all the authors are members of this research project in forest genomics. We thank E. Pouliot, D. Plourde, É. Dussault, P. Labrie, and numerous summer students for technical assistance. We also thank M. Lamothe, J. Fillon and J.-P. Dionne for ideas and support regarding the bioinformatics aspects. We are indebted to P. Watson of Canfor Pulp Limited Partnership who was involved in the early stages of the project and to the team of A. Montpetit at the Illumina genotyping platform of the McGill University and Genome Quebec Innovation Centre. We thank S. Lebihan, C. Collins, and C. Nelson of the British Columbia Microarray Facility at the Vancouver Prostate Centre for printing the microarray slides. This research was supported by funding to the Arborea project from Genome Canada, Genome Quebec, the Canadian Forest Service of Natural Resources Canada (Canadian Wood Fibre Centre and Genomics Research and Development Initiative), and FPInnovations.

Footnotes

  • Received December 9, 2010.
  • Accepted February 22, 2011.

LITERATURE CITED

View Abstract