QTL mapping analyses

Over two years and two climatically different sites (sites #1 and #2, Figure 1), budset timing and height growth were measured in a backcross family that had been previously used to map the genome of black spruce [30, 31]. Given that several adaptive traits were not normally distributed, both parametric and non-parametric QTL mapping approaches were used to delineate genomic regions involved in adaptive trait variation. Over all methods and traits, 35 QTLs that were related to growth and budset timing were detected (QTLs results for each trait are tabulated in Table 1). Among these QTLs, 29 were related to budset variation and 8 were related to growth variation (see Figure 2 as an example of R/qtl results for non-parametric interval mapping). Using the non-parametric interval mapping method (nIM), only one QTL was significant at the 5% genome-wide scale for budset variation. In addition, 4 and 2 QTLs were significant at the 5% chromosome-wide scale for budset and growth variation, respectively. Using interval mapping (IM) with transformed data to fit a normal distribution, 3 and 2 additional QTLs were disclosed as significantly involved in budset and growth variation, respectively, at the 5% chromosome-wide scale. The proportion of variance that was explained (PVE) by a single QTL was low, averaging 8.3% and never exceeding 12%. Assuming an additive model among QTL effects for each trait, the total PVE for one trait varied from 5.2 to 47.9%. However, given the potential complexity of QTL interactions, these values might represent overestimates of the QTL combined effects. Among QTLs, 27 were repeatedly involved in quantitative trait variation over sites, years, or methods (Table 1). Though not implicating a direct role in QTL control, 39 mapped genes were located in QTLs that were found to be significantly involved in quantitative trait variation among the progeny (24 mapped genes located in “budset” QTLs, 5 located in “height” QTLs, and 10 involved in both traits).

Genetic association tests

Relative kinship and association results

Using the genotypes that were obtained for control SNPs, pairwise relative kinship coefficients (K) among individuals were estimated using the method of Loiselle et al.[32]. Average kinship among all individuals was K =0.0004, indicating a very low degree of relatedness in the sampled populations. In addition, possible hierarchical population structure was tested using the software STRUCTURE [33]. The analysis revealed no population genetic structure in the sampled area, since the highest probability that was obtained was for k = 1 group of populations. Such low average kinship (among individuals) and the absence of population structure at the regional scale were not surprising, given that all sampled populations were part of the same glacial lineage with respect to nuclear DNA [12] and cpDNA [34]; also, black spruce is characterised by wind pollination, an outcrossing mating system, and extensive gene flow [35, 36].

For the genetic association tests, we applied the method developed by Yu et al.[37] that is based on a linear mixed-model and takes into account the degree of genetic relatedness among individuals. To correct for multiple testing, a q-value threshold of 10% (FDR [38]) was used to declare significant associations. Using this threshold, none of the control SNPs were significantly associated with trait variations (results not shown). Among the set of candidate SNPs, a total of 34 gene SNPs (representing as many distinct genes) were significantly associated with variation in adaptive traits in black spruce natural populations (Table 2). Twenty-two SNPs were significantly associated with budset variation, while 20 SNPs were significantly associated with height variation. Eight SNPs were involved with both traits, although rarely in the same year, and usually in the following growing seasons (Table 2). The total of 34 significant gene SNPs included 20 and 22 SNPs that were involved in trait variation for site #3 and site #4, respectively, and which were used for association genetic testing; 8 SNPs were detected for the same trait in both sites. In further comparisons of the two sites, no difference could be discerned regarding the average P-values or q-values of significant SNPs between sites (t-test, P = 0.45), but sites differed significantly regarding the average estimated values of PVE (proportion of explained variance) (t-test, P < 0.0001). Indeed, significant SNPs that were detected for site #3, which was situated in milder climatic conditions (Additional file 1: Table S1), rarely explained more than 1% of trait variation, while PVE values for significant SNPs that were detected for site #4, with harsher climatic conditions, frequently exceeded 1%, with several values above 2% (Table 2). When comparing SNPs among years, 11 significant SNPs were corroborated over several years, but fewer SNPs were significantly associated with tree height in the last year of measurement (age 12). Interestingly, data from the final measurement year were not used to conduct bulk segregant analysis that lead to the identification of candidate gene SNPs to be tested (see Methods), suggesting that partially different gene sets are acting at different ages for the same adaptive trait.

SNP	Gene IDa	Gene family	LGb	Site c	Trait	Age	p-value	PVEd	q-value	Correlated climatic factorse
c06511e	GQ02808_G21	AUX/IAA	1	4	budset	7	0.01	1.40	0.06	TWM
c11176n	GQ03816_N20	Unknown	1	3	budset	7	0.08	0.33	0.08	-
c16364e	GQ04002_G01	C2H2 zinc finger	1	3	height	12	0.02	2.23	0.10	PWM*
c02239c	GQ03118_J01	BAM1	3	3	budset	7	0.04	0.11	0.08	PWM***, AAP*
c03870a	GQ0168_K17	C3HC4 Ring finger	3	4	budset	8	0.02	1.95	0.08	MAT*, TWM***
c09889a	GQ02801_L07	NAC	3	4	height	7	0.07	0.74	0.10	MAT*, TWM*, TCM*
c04073m	GQ03313_I07	Zinc ion binding	4	4	budset	7	0.00	2.26	0.03	-
c05720c	GQ03009_F08	C3HC4 Ring finger	4	4	height	6	0.06	1.29	0.06	-
				3	budset	7	0.06	0.41	0.10
				3	height	12	0.01	2.36	0.1
c10254m	GQ0197_M23	Beta-glucosidase	4	3	budset	7	0.01	0.48	0.05	-
c02343e	GQ03102_D15	AP2	5	3	budset	7	0.06	0.33	0.08	TWM*, DGD*
c04662f	GQ0011_L04	AP2	5	4	height	7	0.08	0.73	0.10	-
c09406f	GQ03511_H09	LEA	5	3	budset	7	0.06	0.41	0.10	TCM*, TWM**
c09898f	GQ02813_D01	RAP2	5	4	height	6	0.01	2.01	0.03	MAT***, TWM***, TCM*
				4	height	7	0.00	2.53	0.03
c10140e	GQ04111_I24	AP2	5	4	budset	7	0.00	2.85	0.01	-
				4	budset	8	0.03	1.70	0.08
				3	budset	7	0.02	0.24	0.05
c12722e	GQ03309_B18	U-box superfamily	5	3	budset	7	0.07	0.11	0.08	-
				3	height	12	0.03	2.16	0.10
c14320f	GQ0178_D07	transcription factor	5	4	height	6	0.02	1.70	0.03	TWM
				4	height	7	0.03	1.67	0.04
				3	height	12	0.01	2.40	0.10
c04312b	GQ02754_G07	C2H2 zinc finger	6	4	height	6	0.02	1.85	0.03	-
				4	height	7	0.02	1.72	0.04
c04554n	GQ0201_E19	MYB (10)	7	4	height	6	0.03	1.69	0.03	MAT*, TCM*
				4	height	7	0.05	1.41	0.05
c04671m	GQ03310_H09	PLATZ	7	3	budset	7	0.06	0.65	0.08	TCM*
c07142s	GQ0171_A23	GASA	7	4	budset	8	0.05	1.42	0.08	-
c08080a	GQ0082_F08	MYB (5)	7	4	budset	8	0.09	1.16	0.10	TWM*, DGD*
				3	budset	7	0.05	0.12	0.08
c08085m	GQ03515_N01	Scarecrow-like	7	4	height	7	0.06	0.83	0.07	-
				3	budset	7	0.09	0.17	0.08
				4	height	6	0.10	1.07	0.09
c13572n	GQ0178_D07	Homeobox	7	4	height	6	0.02	1.86	0.03	MAT*, TWM*, TCM*
				4	height	7	0.03	1.67	0.04
				3	budset	7	0.08	0.34	0.10
c15051g	GQ0195_G17	NAM	8	3	budset	7	0.06	0.19	0.10	-
c03650g	GQ03503_G05	B-box zinc family	9	4	budset	8	0.04	1.50	0.08	TCM*, MAT**
				4	height	6	0.04	1.47	0.05
				3	budset	7	0.04	0.36	0.08
c05937s	GQ0207_N03	Chaperone DnaJ Heat	9	3	height	6	0.01	1.29	0.06	DGD*, TWM*, PWM*
		Shock Prot.		3	height	7	0.02	1.12	0.03
				3	height	12	0.04	0.99	0.01
c09063f	GQ02825_P07	PATATIN	9	4	budset	8	0.08	1.21	0.10	TCM*
				4	height	7	0.10	0.63	0.10
c13855n	GQ0255_K05	MYB (11)	9	4	budset	8	0.03	1.65	0.08	MAT***, TCM*, TWM***
c02623e	GQ03208_E14	HLH	10	4	budset	8	0.06	1.38	0.08	-
c09573a	GQ02904_H04	HTH	11	4	height	7	0.06	0.84	0.10	-
c05535s	GQ0206_M13	WRKY	12	4	height	7	0.05	0.90	0.07	-
c07602a	GQ0013_G12	MYB domain	12	3	height	6	0.01	0.85	0.06	-
c15328f	GQ0071_N15	MYB	12	4	budset	7	0.03	1.71	0.08	-
				4	budset	8	0.00	2.77	0.03

^a Gene ID follows the nomenclature of the white spruce gene catalogue GCAT [39].

^b LG: Linkage group following the nomenclature of [31].

^c Site numbers as indicated on Figure 1.

^d PVE: Proportion of variance explained.

^e Climatic factors significantly correlated with the allele frequency variation of significant gene SNPs: PWM: precipitation of the wettest month; AAP: annual amount of precipitation; TWM: temperature of the warmest month; DGD: annual number of degree-days above 5°C; TCM: temperature of the coldest month; MAT: mean annual temperature. *: P ≤ 0.05; **: P ≤ 0.01; ***:P ≤ 0.005.

Gene nomenclature conformed to the most recent white spruce gene catalogue (GCAT, [39], Table 2), and the gene SNPs that had been significantly associated with trait variations were spread over the 12 chromosomes of black spruce [31], except for chromosome # 2. Comparison between association test and QTL results showed that 17 SNPs that had been identified using the association tests were also located in QTLs that were delineated in the present work, including seven SNPs related to the same character. The associated SNPs that had been detected were representative of a large array of gene families, including more frequently occurring families such as MYB and AP2, and more interestingly, families such as Chaperone DNAj proteins, which are already known for their role in reactions to abiotic stress (Table 2). A total of 17 SNPs that were associated with adaptive trait variations displayed significant correlations between allele frequency variation and climatic factors for the population origins of the tested material (Table 2, see Figure 3), of which four were highly significant (P < 0.005).

Quantitative trait loci

A total of 22 regions of the P. mariana genome were detected as being potentially involved in budset and growth variation, which is lower than the 52 QTLs that had been found in the congeneric white spruce for each of these traits by Pelgas et al.[13] which used two pedigrees and larger sets of progeny. A relatively small number of progeny (283) was surveyed in the present study and only one biparental family was used to map QTLs. Given that the progeny were not clonally replicated, in contrast to the study conducted by Pelgas et al.[13], it was not expected to identify as many QTLs with as much precision. However, the main goal of the present QTL survey was to obtain an exploratory overview of the genomic architecture of adaptive traits in black spruce for comparative purposes.

Trait variations were not normally distributed, so that the guidelines of Asins et al.[40] were followed and, consequently, both parametric and non-parametric approaches were used to identify QTLs. The three different methods yielded overlapping complementary results. The non-parametric methods (nIM and K-W tests) were notably more efficient in QTL detection than the parametric approach (IM), despite their expected lower power to detect QTLs. However, the application of IM still yielded valuable information, such as estimates of the proportions of explained variance. Despite the limitations that were imposed by the available data, four QTLs were repeatedly identified over years or sites, thereby more likely representing genomic regions encompassing polymorphisms that were more central to trait variation. Such a restricted set of QTLs that is shared across years or environments has also been observed in the congeneric P. glauca[13].

On average, 8% of the observed phenotypic variance could be explained by a single QTL. This average proportion of explained variance (PVE) seemed low compared to several published values in conifer QTL mapping studies of traits related to growth or adaptive traits [41–43]. Yet, many QTL studies showing high PVE for single QTLs have been performed using interspecific crosses [44–47]. Given that phenotypic variation in such crosses may represent the sum of trait variations originating from both species, QTLs having effects on these adaptive traits would harbour higher PVE estimates accordingly to the average trait differences between the two species. In contrast, the average amplitude of PVE values that were estimated in the present study is consistent with results obtained in other QTL studies on phenology and growth trait variation which had been assessed with intraspecific crosses [13, 43, 48, 49]. Even though such PVE estimates may have been biased upward due to the Beavis effect [50], the congruence with results from other QTL studies relying on larger progeny [13, 48], which should be less sensitive to such effect [50], suggests that PVE values were quite accurately estimated in the present study. These results suggested that many QTLs are truly involved in adaptive variation. We noted that three QTL regions were involved in both budset and growth variations, suggesting pleiotropic effects. A similar pattern was also observed among P. glauca QTLs for phenology and growth [13]. Such patterns are expected, given that significant positive correlations between growth and budset are generally observed in trees as reported previously in both P. glauca and P. mariana[25, 28]. In other words, the later the budset, the longer the growing season and, consequently, the more growth is generally accumulated.

SNPs associated with adaptive trait variation

A total of 34 SNPs from as many distinct genes were significantly associated with adaptive trait variations in genetic association tests, representing 65% of the tested candidate SNPs. Such a rate of positives is high, exceeding by more than an order of magnitude the results obtained in other genetic association studies that are based on candidate genes by more than an order of magnitude (e.g.[18, 20, 51]), and where there had been no a priori screening of the tested gene polymorphisms. Given the much smaller number of tests that were performed with a limited set of candidate polymorphisms, the issue of multiple testing and potentially high rates of false positives that are associated with intensive testing becomes less of a concern. Although “next-gen” sequencing technologies will allow detection of many more polymorphisms at an affordable cost in the near future, with approaches such as genotyping-by-sequencing [52], statistical issues that are related to high rates of false positives when large numbers of tests are conducted will remain a recurrent theme. As shown here, a candidate SNP approach should contribute to reducing these concerns and the number of tests to be performed, especially when quantitative characters implicating large numbers of genes are considered, or for organisms characterised by very large genomes such as conifers [53]. The outlier analysis and bulk segregant analysis of sequence signatures that were used to identify candidate SNPs were quite efficient in reducing the number of candidate SNPs to be tested by more than an order of magnitude, in both cases (see Methods). The several hundred genes that were involved in the sequencing and preliminary screening procedure for segregating SNPs between opposite pools of extreme individuals represented quite a large number of screened candidate genes, given that we relied on traditional PCR and Sanger sequencing to discover segregating SNPs. New methods relying on massive sequencing of the gene space that use gene sequence capture approaches and next-gen sequencing should make it easier to discover much larger numbers of candidate polymorphisms segregating for a variety of phenotypic characters. In various spruce species, such strategies are currently employed, which involve gene capture and sequencing for a variety of pools using oligo arrays that represent most of the transcriptome of the spruce genome [39].

While screening for candidate gene SNPs was conducted using bulk segregant analysis when plants were 6- and 7-years old, more SNPs were found to be significant in association testing around these ages (6-, 7- or 8-years old) than at a later age of 12 years (Table 2). This trend suggests that different genes and SNPs (and physiological processes) could be involved at different ages, as the transition from juvenile to mature adult is taking place. This scenario is additionally supported by the fact that data from 12-year-old plants were not used in the screening procedure to identify candidate gene SNPs from pool segregant analysis, and that the number of significant genetic association tests was the lowest for field data at that age. A similar trend was not observed in a QTL study of adaptive characters that had been assessed over several years in P. glauca, probably because the trees were still at the juvenile stage of development [13]. In the present study, genetic correlations between characters at age 12 and ages 6, 7 and 8 years were generally modest for both sites (data not shown), further pointing to partially different sets of genes and polymorphisms involved. These modest correlations are consistent with previously published estimates from a clonal test in P. mariana[54].

In our study, each SNP that was associated with trait variation explained a low proportion of variance (PVE ~ 1.25%), which significantly differed between the two study sites. Average PVE was generally higher for site #4 than for site #3. Site #4 was characterised by harsher climatic conditions at higher altitude (Additional file 1: Table S1), with freezing temperatures occurring earlier in the season than at site #3, the latter site being located at lower elevation and lower latitude. Since budset timing is also significantly related to temperature conditions [23], it was likely that quantitative trait differences between adapted and maladapted trees would have been amplified on the site with harsher climatic conditions, leading to greater ability in detecting significant associations. With respect to this likely trend, the variance in quantitative traits was generally greater on the site with harsher climatic variation (data not shown).

The low PVE values that were observed with significant SNPs suggested that a large number of adaptive polymorphisms from various genes and gene families are involved in the genetic control of adaptation. These low PVE values were not surprising given that similar values have been observed in genetic association studies involving other tree species (3% in [15]; <1% in [55]; 3% in [17]; 0.7-5.4% in [18]; 3-5% in [20]). This trend was consistent with that found for QTLs of the same characters (although lower PVE values were observed in association testing), and in agreement with quantitative genetics theory, which predicts that genetic variation in these quantitative traits is controlled by many genes, each exerting small effects [56, 57]. Because linkage disequilibrium usually decays rapidly within gene limits in natural populations of boreal spruces [10, 11], the co-occurrence of several physically linked adaptive polymorphisms that would explain a higher proportion of the variance was unlikely.

Several lines of evidence supported the involvement of these significant SNPs in adaptation. First, half of these polymorphisms displayed allele frequency variation significantly correlated with climatic variations (Table 2, Figure 3), hence corroborating, with an extensive population sampling, the trends that had been previously observed in black spruce molecular adaptation studies [26]. In the absence of population structure, these correlations suggested that the alleles are under divergent selection. However, the remaining SNPs should not be considered as false positives given that variation in SNP frequencies may not follow a linear distribution [12]. Although the distribution of an adaptive trait is expected to follow the selective environmental variable, adaptive allele frequency variation along the environmental gradient may follow various and sometimes barely predictable distributions given the complexity of interacting genes and alleles that likely govern adaptive trait variation [12, 42].

This idea of complexity of gene and allelic interactions is also reinforced by the large diversity of families and functions observed among the genes that were found to harbour significant SNPs. We noted that 53% of significant SNPs belonged to gene families previously found to harbour genetic polymorphisms and expression patterns that were related to stress resistance and adaptation in other conifers, trees or plant model organisms. For instance, genes of the MYB family have already been reported in several studies as transcription factors that are involved in growth control, response to cold stress and organ differentiation, and related to budset in P. glauca[20, 58–60]. In addition, MYB members have also been discovered in genetic association studies of cold hardiness in Pseudotsuga menziensii var. menziesii[17] and Picea sitchensis[18]. A SNP in a BAM gene was also associated with ring growth variation in a previous association study in P. glauca[20]. As budset timing and cold resistance are preceded by growth cessation [23], we expected to find genes regulating growth development. As was the case for genes from the AP2 family, which have been identified previously in the gene network activating plant freezing tolerance in Arabidopsis[61] and Triticum aestivum[62]. Similarly, variation in the expression of genes from the C3HC4 RING family has been observed under drought and low temperature stress in Arabidopsis and rice [63, 64]. In the latter species (Oryza sativa L.), a protein of the WRKY family was also found to be activated by a Heat Shock protein under drought stress [65]. Dnaj Heat Shock proteins are produced under diverse abiotic stresses [66] and both WRKY and Heat Shock families were represented by a gene SNP that was associated with height variation in the present study. Among other significant gene SNPs that were found in the present study, the NAC, LEA, PATATIN and RAP2 gene families were also represented. NAC proteins are over-expressed in response to drought and salinity stresses in rice [67], as well as in drought responses of Arabidopsis[68]. In addition, several LEA genes are known to improve resistance to dehydration and freezing temperatures in plants [69] and were reported in a drought adaptation study in Pinus pinaster Aiton [70]. It should be also noted that a PATATIN gene has been found to be stimulated by drought stress in Vigna unguiculata (L.) Walp. and a RAP2 gene was found to be up-regulated by drought stress in Arabidopsis. Hence, a wide variety of gene families was found, indicative of the large number of physiological processes involved in budset and growth.

Of the 34 SNPs significantly associated with adaptive trait variation, seven were also located in a QTL for the same character. Given the known relationships between budset timing, the length of the growing season and growth (see above), it is plausible that several SNPs that were associated to one trait were located in a QTL for the other trait. Considering such cases, 17 significant SNPs (50%) were located in QTLs. Such corroboration would support the identified polymorphisms as true positives. Furthermore, the two analytical approaches should not be expected to be congruent but complementary. QTL mapping permits the identification of polymorphisms involved in trait variation in one or few pedigree populations, thereby reducing background genetic variations as well as the number and complexity of interacting alleles, while genetic association testing permits the identification of polymorphisms of more general value over a wider range of genetic backgrounds.

In total, 65% of the significant SNPs that had been identified by genetic association tests were corroborated by additional indirect evidence of various sources, as indicated above (relationships with climatic factors, adaptation QTLs and key gene families). This evidence highlights their relevance for possible in situ monitoring of genetic variation for adaptation or for inclusion in marker-aided breeding selection schemes. The remaining gene SNPs could still represent polymorphisms of high adaptive value, even though additional research would be required to identify their specific role in adaptation.

QTL mapping analyses

The QTL analysis was conducted using a backcross family (#9920002: ♀11307-03 [♀83 x ♂425] × ♂425) of 283 individuals already used to map the genome of black spruce [30, 31]. The genome of black spruce contains 12 chromosomes corresponding to as many different linkage groups, as is the case for most other Pinaceae and all spruces that have been mapped so far [10, 30, 31, 77]. The controlled cross used for mapping was performed in 1999. The progeny was seeded in greenhouse in 2000 and planted in 2002 on site #1 (46°38′N, 72°00′W; 35 m elevation; Figure 1). Clones (rooted cuttings in 2002) of these individuals were subsequently planted in 2004 on site #2 (47°45′N, 70°50′W; 792 m elevation; Figure 1), which was characterised by harsher climatic conditions at higher altitude (Additional file 1: Table S1). However, no clonal replicates were available on each site. Budset timing was assessed weekly between mid-July and October (~13 measurements) and scored following a three-stage morphological index as in [25]: (0) absence of terminal bud, (1) presence of a white terminal bud, (2) brownish terminal bud, and (3) terminal bud with scales ready to endure winter conditions. For each tree, these scores were averaged over the measurement period. Growth was assessed by measuring terminal shoot elongation of each tree at the end of the growing season. These measurements were collected on site #1 for 254 and 244 progeny in 2006 and 2007, respectively and on site #2 for 155 and 142 progeny for the same years (differences in sample size between years was due to mortality).

As suggested for traits not normally distributed [40], the QTL analysis was conducted following parametric and non-parametric approaches. First, a traditional interval mapping approach was applied to the transformed data to fit the normal distribution required for this method. It was implemented in MapQTL 5.0 [78], and QTLs having LOD scores greater than a threshold value that had been determined by a permutation test were retained (1000 permutations were applied at the genome-wide level or each linkage group separately). Second, non-parametric interval mapping was applied to untransformed phenotypic data using the R/qtl package [79] for R (R Development Core Team, 2011). As was the case for the parametric interval mapping approach, the minimum LOD score threshold was determined using a permutation test that allowed the identification of significant QTLs. Finally, a Kruskal-Wallis test was applied to phenotypic data to test each marker, keeping in mind that this represented a multiple testing procedure. Loci repeatedly identified over years, sites and methods were considered as QTLs of putative interest.

Genetic association testing