Linking genotype to phenotype to identify genetic variation relating to host susceptibility in the mountain pine beetle system

Abstract Identifying genetic variants responsible for phenotypic variation under selective pressure has the potential to enable productive gains in natural resource conservation and management. Despite this potential, identifying adaptive candidate loci is not trivial, and linking genotype to phenotype is a major challenge in contemporary genetics. Many of the population genetic approaches commonly used to identify adaptive candidates will simultaneously detect false positives, particularly in nonmodel species, where experimental evidence is seldom provided for putative roles of the adaptive candidates identified by outlier approaches. In this study, we use outcomes from population genetics, phenotype association, and gene expression analyses as multiple lines of evidence to validate candidate genes. Using lodgepole and jack pine as our nonmodel study species, we analyzed 17 adaptive candidate loci together with 78 putatively neutral loci at 58 locations across Canada (N > 800) to determine whether relationships could be established between these candidate loci and phenotype related to mountain pine beetle susceptibility. We identified two candidate loci that were significant across all population genetic tests, and demonstrated significant changes in transcript abundance in trees subjected to wounding or inoculation with the mountain pine beetle fungal associate Grosmannia clavigera. Both candidates are involved in central physiological processes that are likely to be invoked in a trees response to stress. One of these two candidate loci showed a significant association with mountain pine beetle attack status in lodgepole pine. The spatial distribution of the attack‐associated allele further coincides with other indicators of susceptibility in lodgepole pine. These analyses, in which population genetics was combined with laboratory and field experimental validation approaches, represent first steps toward linking genetic variation to the phenotype of mountain pine beetle susceptibility in lodgepole and jack pine, and provide a roadmap for more comprehensive analyses.

The false positive discovery rate depends on how well the study system fits the demographic assumptions of the method being used (Excoffier, Hofer, & Foll, 2009;Thornton & Jensen, 2007;Whitlock & Lotterhos, 2015). But, the number of false positives can be reduced by considering loci that have their outlier status verified across multiple statistical approaches (De Mita et al., 2013;Gosset & Bierne, 2012;Narum & Hess, 2011;Nunes, Beaumont, Butlin, & Paulo, 2011). Loci verified in this manner can then be validated by testing for signatures among the same set of loci in different sets of individuals/populations (Bonin, Taberlet, Miaud, & Pompanon, 2006;Stinchcombe & Hoekstra, 2008). Ultimately, however, identification of true positives requires experimental or functional validation (Salvi & Tuberosa, 2005).
Verifying potential candidates by linking genetic variation to phenotype involves manipulating the study system using common gardens, field or laboratory experimentation (Kingsolver, 1996;Vignieri, Larson, & Hoekstra, 2010;Wikelski, Spinney, Schelsky, Scheuerlein, & Gwinner, 2003), quantitative trait locus (QTL) analysis of pedigreed material (Clutton-Brock & Sheldon, 2010;Slate et al., 2010), or association mapping (Parchman et al., 2012). Many of these options are logistically challenging for nonmodel systems, particularly long-lived species. An alternative approach is to use multiple lines of investigation, for example, using both gene expression data and population genomics to provide a first indication as to whether genes identified through population genomics have the potential to contribute to adaptive variation (Levy & Borenstein, 2012;Li, Costello, Holloway, & Hahn, 2008;Vasemägi & Primmer, 2005). Positive results obtained with such an analysis can be followed up with more in-depth functional and sequence analyses of candidate genes. We will term this "complementary validation." This approach can validate candidates without requiring experimental manipulation of the locus, that is, through targeted mutation or transgenesis of the gene, which is difficult to achieve for many nonmodel species.
Here, we apply a complementary validation approach to study host factors associated with tree defense against Dendroctonus ponderosae (Hopkins, mountain pine beetle [MPB]) attack in Pinus contorta (Dougl. ex Loud. Var. latifolia, lodgepole pine) and P. banksiana (Dougl., jack pine). These are sister species that hybridize in north-central Alberta and the Northwest Territories (Burns et al., submitted;Cullingham, James, Cooke, & Coltman, 2012). Both species are important for forest economy (Law & Valade, 1994;Lee, 1971) and provide habitat for an array of species (Kirk & Hobson, 2001;Martin, Norris, & Drever, 2006), thus making them important study organisms. Since the early 2000s, these tree species have received critical attention due to the ongoing outbreak of MPB. While lodgepole pine is a native host of this opportunistic bark beetle and its fungal associates (Safranyik & Carroll, 2006), the unprecedented scale of this most recent outbreak (Aukema, McKee, Wytrykush, & Carroll, 2016;de la Giroday, Carroll, & Aukema, 2012) has resulted in considerable range expansion of MPB and their microbiome. This range expansion includes spread into naïve lodgepole pine forests (Burke, Bohlmann, & Carroll, 2017;Cudmore, Bjorklund, Carroll, & Lindgren, 2010) and host expansion into the jack pine of Canada's boreal forest (Cullingham et al., 2011). Given the potential for this devastating forest insect pest to establish endemic populations in these novel habitats, understanding how pine susceptibility varies on the landscape is a key component to predicting risk of further spread during this and future outbreaks. Identifying loci that correlate with susceptibility provides one avenue to shed light on this important question.
Despite their ecological and economic importance, lodgepole pine and jack pine are nonmodel species that present considerable challenges for experimental manipulation typically used in more tractable systems to validate candidate genes. Both common gardens and QTL analyses are feasible, but it can take a number of years for trees to reach reproductive maturity that is required for seed collection, and also to reach a sufficient age at which reliable phenotypic measurements can be taken. Given the rapid decay of linkage in Pinus genomes, robust QTL analysis requires numerous crosses and dense genetic markers (Neale & Salvolainen, 2004). In addition, genetic transformation of lodgepole and jack pine has never been reported. Therefore, techniques that are often employed to validate candidate genes in model systems like Arabidopsis thaliana (Bergelson & Roux, 2010;Brachi et al., 2010) are infeasible for lodgepole and jack pine.

K E Y W O R D S
adaptive variation, complementary validation, gene expression, Grosmannia clavigera, jack pine, lodgepole pine, mountain pine beetle, outlier validation Accordingly, in an effort to identify pine genetic variation that may be of adaptive importance in mediating the host's response to MPB, we have taken advantage of genetic and genomic resources that have been assembled for lodgepole and jack pine (Arango-Velez et al., 2014;Cullingham, Cooke, & Coltman, 2013;Cullingham et al., 2011;Hall et al., 2013), combining population genomics with gene expression analysis of the host species. We characterized over 800 lodgepole pine, jack pine, and their hybrids across 58 locations at a subset of 96 SNP loci that were previously analyzed at 17 locations (Cullingham et al., 2014). Seventeen of these loci were identified as having statistically significant signatures of local adaptation (Cullingham et al., 2014). Building on these findings, our first objective was to validate these candidates using outlier and environmental correlation analysis similar to Cullingham et al. (2014). Our second objective was to link the validated candidates to MPB attack status of adult trees by comparing genotype frequencies among freshly attacked and unattacked trees sampled within the same stands. Our third objective was to carry out in silico annotation of validated candidate locus sequences to determine whether SNPs were located within protein-coding or untranslated regions (UTRs), and whether SNPs had the potential to impact protein function. Our final objective was to determine whether candidate loci from the above analyses constituted part of the tree's transcriptomic response to attack by using quantitative reverse transcription-PCR (qRT-PCR) in response to wounding, and inoculation with the MPB fungal associate Grosmannia clavigera ([Robinson-Jeffrey and Davidson] Zipfel, de Beer and Wingfield; Whitney, 1971). These data and analyses, when combined, provide independent evidence for the putative roles of candidate genes in adaptive traits.

| Population samples
Samples used in the population genetics component of this investigation were selected from a larger set of samples used for previous microsatellite (Cullingham et al., 2011(Cullingham et al., , 2012 and SNP (Cullingham et al., 2013)  To obtain seedlings of a sufficient size for analysis, seeds were F I G U R E 1 Distribution of sample sites used for this study ("New sites"), together with the sites used in a previous study of lodgepole and jack pine using population genetic methods to identify adaptive candidates ("Original sites," Cullingham et al., 2014). The inset includes the sampling location of all trees used in a paired analysis to examine whether any of the identified adaptive candidates showed an associated with attack status. Trees with >40 attacks were sampled, and an age-matched health tree at the same location was also sampled. Predicted distributions of lodgepole, jack pine, and their interspecific hybrids are included for reference (Burns et al., submitted) germinated following the protocol outlined in Cullingham et al. (2012). Germinated seedlings were flash-frozen in liquid nitrogen and stored at −20°C prior to DNA extraction. GPS locations were collected for field trees, and centroids of locations were used for site-level analyses ( Figure 1).
Twelve of the locations in British Columbia and Alberta were situated within the active MPB outbreak during the 2007-2008 sampling period. At these locations, foliage was sampled from both unattacked and naturally MPB attacked pines. Attacked pines typically had sustained >40 attacks, as determined by counting pitch tubes.
In Alberta, attacked trees had been identified the prior autumn via aerial surveys by Alberta Agriculture and Forestry (formerly Alberta Sustainable Resources Development). In all cases, MPB attack status of sampled trees was confirmed by inspecting gallery structure and identification of MPB within galleries. Following identification of a suitable MPB attacked pine for sampling, an unattacked, healthy tree of similar diameter by breast height (DBH) was identified for sampling no more than 250 m from the MPB attacked tree. We refer to these as paired attacked/unattacked trees. Between five and 22 attack/unattacked pairs of trees were collected per location. All field samples were stored on ice or at subzero conditions (in the case of winter sampling) until transported to the laboratory, and were kept at −20°C prior to DNA extraction. Genomic DNA was isolated according to Cullingham et al. (2011). Cullingham et al. (2013) previously typed 546 lodgepole pine, jack pine, and lodgepole × jack pine hybrids samples from 17 locations at 472 SNP loci using the Golden Gate assay (Illumina, San Diego, CA). From these data, 28 outlying SNP loci were considered adaptive candidates (Cullingham et al., 2014). To validate these findings at a broader scale of locations, we selected 96 of those 472 SNP loci for Sequenom SNP typing (Agena Biosciences, San Diego, CA) including 17 of the adaptive candidates and 79 putatively neutral loci. The adaptive loci we selected as strong candidates were statistically robust (i.e., identified consistently across outlier detection methods) and/or showed changes in transcript abundance in response to G. clavigera inoculation in a microarray expression dataset E. K. Cooke, unpublished). We included 44 samples from our previous genotyping effort to assess congruence between the two genotyping technologies. The same flanking sequence data used to develop the Illumina SNP chip assay (Supporting Information Table   S1) were used to design the Sequenom MassARRAY assay by Agena Given that a number of our sample sites occur within the area of hybridization between lodgepole and jack pine, we ran structure (Falush, Stephens, & Pritchard, 2003;Pritchard, Stephens, & Donnelly, 2000) to identify ancestry of all individuals. We used the admixture model, setting the number of clusters (K) to two; 500,000 iterations were run, discarding the first 100,000 as burn-in. We ran the model five times and averaged the q-values across runs to assign ancestry.

| Outlier detection
For consistency, we applied three of the four detection methods that we used in the previous study (Cullingham et al., 2014) together with OutFlank (Whitlock & Lotterhos, 2015), which was not publicly available at the time of the original publication. We chose not to use matsam (Joost et al., 2007), because it does not incorporate population structure information and results in a high proportion of false positives (Cullingham et al., 2014). For all analyses, we completed outlier detection on the pure species separately.
We first used the hierarchical model in arlequin (Excoffier et al., 2009;Excoffier & Lischer, 2010), considering loci to be significant outliers when p ≤ 0.01. Second, we used bayescan (Foll & Gaggiotti, 2008), using a burn-in of 10,000 iterations, thinning interval of 50, sample size of 10,000, and prior odds of 1,000 to limit the number of false positives (Lotterhos & Whitlock, 2014). Third, we used bay-env2 (Coop, Witonsky, Rienzo, & Pritchard, 2010;Gunther & Coop, 2013) to detect loci with significant correlations to the environment. Climate data were obtained from ClimateNA v5.10 (http:// tinyurl.com/ClimateNA), which is based on methodology described by Wang, Hamann, Spittlehouse, and Carroll (2016). We used the following variables averaged over a 30-yr period : temperature (annual temperature, annual maximum temperature, annual minimum temperature, mean temperature wettest quarter, mean temperature driest quarter, mean diurnal range, seasonal temperature, maximum temperature warmest period, minimum temperature coldest period, temperature annual range), precipitation (annual precipitation, precipitation driest period, precipitation wettest period), end of growing season, and start of growing season. We also included longitude, latitude, and elevation. We estimated three covariance matrices to correct for underlying genetic population structure, and ran three iterations with each covariance matrix for every SNP against environmental variables using 50,000 MCMC steps to estimate BayesFactors for each SNP × environment association. We used Jeffrey's scale of evidence (Jeffreys, 1961) to determine the significance of the BayesFactors (3-10 = substantial support, 10-30 = strong support, 30-100 = very strong, >100 = decisive support). Only those loci with consistent values across a majority of iterations were considered. The final method we used was OutFlank (Whitlock & Lotterhos, 2015). This method estimates the distribution of F ST values across loci to better identify statistical outliers by limiting the number of false positives. We verified significant association of loci across the detection methods. Loci that were significant across the majority of methods (≥3 out of the 4 methods), and were previously identified as significant (Cullingham et al., 2014) were considered to be validated candidate loci, and were the focus of further analyses.

| Phenotype association
For candidate loci validated using the multiple outlier analysis approaches described above, we compared genotype frequencies between attacked and unattacked trees using a chi-squared test to identify whether any of the genotypes were associated with attack status. For these analyses, we pooled samples typed in this study (N = 47) with previously genotyped trees (N = 309, Cullingham et al., 2014), totaling 178 pairs of attacked/unattacked trees across 12 different locations (Figure 1 inset).

| In silico sequence analysis
tBLASTx queries of sequence databases (NCBI nr, Altschul, Gish, Miller, Myers, & Lipan, 1990) provided cursory sequence annotation for all of the loci we genotyped (Cullingham et al., 2014). To improve the veracity of the annotation, and to identify whether SNPs in these loci result in a predicted amino acid change, we obtained the coding sequence of the protein of interest from available plant species.
These coding sequences were obtained from NCBI GenBank using blastx (Altschul et al., 1990). We used these sequences to generate a more complete alignment and to perform phylogenetic analysis. Alignment was completed using ClustalW (Thompson, Higgins, & Gibson, 1994). The phylogenetic analysis was completed on the aligned deduced amino acid sequences in MEGA6 (Tamura, Stecher, Peerson, Filipski, & Kumar, 2013) using maximum likelihood (ML) with the default settings (bootstrap test of phylogeny, Jones-Taylor-Thorton substitution model, and tree inference using nearest neighbor interchange). Using the alignment, we identified whether the SNP in the validated candidate locus resulted in a predicted amino acid change. If the SNP resulted in a predicted amino acid change, we then generated an in silico predicted three-dimensional structure of the protein using Phyre 2 (Kelley, Mezulis, Yates, Wass, & Sternberg, 2015) and SuSPect (Yates, Filippis, Kelley, & Sternberg, 2014) to predict potential impacts of the amino acid change on protein function.

| Gene expression analysis samples
To assess whether validated candidates identified may be involved in the physiological response to MPB attack, we analyzed changes in transcript abundance of the candidates in lodgepole and jack pine seedlings inoculated with G. clavigera. Seedlings were used for these analyses rather than mature trees to better control for the betweenindividual variance in transcript abundance that is typical of experiments that are conducted with nonclonal material. Bark for analysis was harvested from second-year seedlings of lodgepole pine (Alberta control. Seedlings (n = 5-11 replicates per treatment combination) were harvested at 1, 7, and 14 days postinoculation (dpi). Samples were derived from the same randomized complete block design experiment described in detail in Arango-Velez et al. (2015).

| Transcript abundance profiling
Total RNA was extracted from ∼100 mg tissue using the cetyltrimethylammonium bromide (CTAB) protocol of Chang, Puryear, and Cairney (1993) with modifications by Pavy et al. (2008) First, we analyzed the C t data for stability of the reference genes, and to detect outliers using BestKeeper v1.0. Following verification of the reference genes, we used the modeling approach developed by Matz, Wright, and Scott (2013) to analyze the expression data. This method builds on the generalized linear mixed model approach introduced by Steibel, Poletto, Coussens, and Rosa (2009) allowing the incorporation of both fixed and random effects. We used the normalization model in the MCMC.qpcr package (Matz et al., 2013), which is the most powerful model using the information from the reference genes. We included time × treatment as a fixed effect, and sample, genotype (see below), block, and run as random effects. We used a burn-in period of 10,000 iterations, and an additional 30,000 iterations for data collection, the number of iterations was determined to ensure the model reached consistency.

| Allelic resequencing
All lodgepole individuals from the transcript abundance experiment were genotyped at the validated candidates (Lodgc1087 and TA B L E 1 SNP loci identified as a statistical outliers in jack pine and lodgepole pine

| Genotyping
Of the 96 loci included in the Sequenom genotyping assay, 87 were retained following quality control analysis of the genotypes,

| Outlier detection
In lodgepole pine, we identified six significant outlier loci. Of these, only two were verified across three or more of the detection methods and were also detected previously (Cullingham et al., 2014), Lodgc1087 and Lodgc2304 (Table 1). We consider these two loci to be validated candidates. In jack pine, we detected five outlier loci, but these were detected with only one method (arlequin). The validated candidate loci detected in lodgepole pine were monomorphic in jack pine.

| Phenotype association
We found an association with successful MPB attack at one of the two validated candidate loci. For Lodgc1087, genotype was significantly associated with attack status ( 2 df = 2 = 6.23, p = 0.044), where the homozygous recessive genotype was more likely to be attacked than expected (Figure 2). We found no significant association for Lodgc2304 ( 2 df = 2 = 0.089, p = 0.956). Association in jack pine was not tested as these two loci were monomorphic in jack pine. To understand the spatial distribution of allelic variation at the two validated candidates, we created interpolated maps of the frequency of the minor allele using all genotype data in Alberta and British Columbia (this study and Cullingham et al., 2014). Interpolation was completed using the kriging option in ArcMap 10.5 (Figure 3).
The candidate SNP was within the coding region and conferred an amino acid change from glutamic acid to aspartic acid, both belong F I G U R E 2 Distribution of genotypes (A-major allele, a-minor allele) that were observed for lodgepole pine individuals recently attacked by mountain pine beetle ("A"), paired with equal aged, location matched unattacked ("U") lodgepole pine at two validated candidate loci, Lodgc1087 and Lodgc2304 (N = 178 pairs). Trees were obtained from 12 different locations in Alberta and British Columbia (Figures 1 and 3) to the acidic amino acid group. This amino acid change is relatively common across species and does not likely result in a change in protein function based on the prediction analysis in Phyre 2 and SuSPect.
Contig Lodgc2304 showed highest sequence similarity to the plastid triosephosphate isomerase (pTPI), which is encoded in the nuclear genome. We confirmed this was the plastid and not cytosolic version F I G U R E 3 Interpolated map of the minor allele for Lodgc1087 and Lodgc2304 across sample sites corrected for pine distribution (Yemshanov, McKenney, and Pedlar (2012). Sites that were used to assess relationship of loci with attack status are indicated as "paired sites." Included in the map is an approximation of the historic range of mountain pine beetle based on attack data from 1959 to 1999 (British Columbia Ministry of Forests, Lands, and Natural Resources [http//www.for.gov.bc.ca/ftp/HFP/external/!publish/]) F I G U R E 4 Normalized transcript abundance levels for proteasome subunit α7 ("Proteasome subunit") and plastid triosephosphate isomerase ("Triosephosphate isomerase") under three conditions (c, control; f, fungal inoculation with Grosmannia clavigera; and w, wound), across time (1, 7, and 14 dpi). Significant fold changes are indicated by an "*," where p ≤ 0.05. Normalized abundance was estimated by dividing molecule abundance for the target gene by the geometric mean of the reference genes Tif5α, and EF1α of this protein by aligning with both isoforms from a number of species (Supporting Information Figure S2). Based on the amino acid alignment, the outlier SNP was within the coding region and causes an amino acid change from glycine to alanine. These are functionally different amino acids, and the prediction/mutation analysis in Phyre 2 and SuSPect suggested a moderate probability of a change in protein function.

| Transcript abundance profiling
Quantitative RT-PCR with transcript-specific primers was used to measure transcript abundance corresponding to proteasome subunit α7 and pTPI in lodgepole and jack pine seedlings challenged with G. clavigera. We identified four outlier samples in lodgepole pine and three in jack pine using the sample-based fold-change analysis in BestKeeper. Following outlier removal, we statistically compared expression levels across treatments for each species using the nor-

| D ISCUSS I ON
While studies using population genomics to identify candidate loci are on the rise, we are not yet seeing a concomitant increase in validation studies (Shafer et al., 2015). Many studies use functional annotation to validate candidates, which has been shown to be highly biased (Pavlidis, Jensen, Stephan, & Stamatakis, 2012). Examining previously identified candidate loci in pine (Cullingham et al., 2014), across a larger set of samples representing a greater geographic distribution, we have validated two loci within protein-coding genes that may relate to host response to MPB. We have used multiple lines of evidence including independent natural, and experimental populations as part of this validation. We have called this approach complementary validation. Through this approach, we have identified a locus in lodgepole pine that could be used, in conjunction with other markers, to predict host susceptibility to MPB, and defined an analytical pipeline to identify additional genetic loci that contribute to host susceptibility.
Finding links between genetic variation and phenotype is a challenging problem in evolutionary biology (Naish & Hard, 2008). We have taken a unique approach to addressing this question, considering attack status within a natural outbreak to represent a phenotype of MPB susceptibility. Through the analysis of paired attacked and unattacked lodgepole pine trees, we have found a significant association with genotype at one of the validated candidates, Lodgc1087 with MPB attack status. To our knowledge, this is the first report of a pine locus associated with MPB host susceptibility in lodgepole pine.
Interestingly, jack pine populations do not carry the risk-associated polymorphism for Lodgc1087.
The landscape-level spatial distribution of the attack-associated allele for Lodgc1087 in lodgepole pine corresponds to the historic MPB climatic suitability classifications (Burke et al., 2017). The allele is at a high frequency where MPB climatic suitability has historically been low or very low (Figure 3), while the frequency of the allele is lowest in regions with very high MPB suitability. Forests classified as having extreme historic climatic suitability exhibit the highest likelihood for four determinants of MPB success: completion of a 1-year life cycle, overwintering temperatures favorable for larval survival, advantageous emergence and dispersal conditions, and spring precipitation as a proxy for tree defense capacity (Burke et al., 2017;Safranyik, Shrimpton, & Whitney, 1975).
Based on population genetics studies of MPB and lodgepole pine (Bentz et al., 2010;Cwynar & MacDonald, 1987;Mock et al., 2007;Wheeler & Guries, 1982) Burke et al. (2017) suggested that the northward post-glacial recolonization of western North America from southern refugia by MPB and lodgepole pine has been asynchronous, meaning that while southern populations of lodgepole pine are likely to have shared a long co-evolutionary history with MPB, northern populations of lodgepole pine are likely to have not. Burke et al. (2017) proposed that the spatial gradient of historic climatic suitability can serve as a proxy for the degree of evolutionary association between MPB and lodgepole pine.
Under this assumption, we postulate that the southern lodgepole pine populations located in high or extreme historic climatic conditions have experienced the greatest degree of evolutionary selective pressure from MPB, resulting in reduced frequencies of the attack-associated minor allele. Conversely, northern populations in which this attack-associated minor allele is present at relatively high frequencies have experienced no selection pressure from MPB.
In silico annotations suggest that both candidates participate in pathways involved in tree responses to environment, and therefore conceivably play direct or indirect roles in tree defense. Lodgc1087, the locus associated with MPB attack status, is within a gene encoding one of the proteasome alpha subunits of the 26S proteasome (Fu et al., 1999). The ubiquitin-26S proteasome machinery functions in targeted protein breakdown, and is highly conserved across eukaryotes (Dreher & Callis, 2007). The protein-coding sequence containing Lodgc1087 exhibited the highest similarity to proteasome subunit α 7 , which acts as a gate to the catalytic chamber of the 20S core protease (Book et al., 2010). Because plants are sessile, they have evolved complex proteolytic systems to manage responses to biotic stresses such as herbivore and pathogen attack, and the proteasome-ubiquitin pathway is considered an important component of these response systems (Dielen, Badaoui, Candresse, & German-Retana, 2010;Small & Vierstra, 2004;Vierstra, 2009). The 20S core protease has been implicated in plant immunity to pathogens (Üstün et al., 2016) and also plays an important role in defense-associated hormone signaling (Vierstra, 2009).
Lodgc2304 is within pTPI, which codes for a key plastid enzyme that is encoded in the nuclear genome (Henze, Schnarrenberger, Kellermann, & Martin, 1994). Plastid TPI interconverts dihydroxyacetone phosphate to glyceraldehyde 3-phosphate, controlling the triose phosphate pool. Plastid TPI and the triose phosphate pool are integral to the Calvin cycle of photosynthesis (Raines, 2003) and to glycolysis (Plaxton, 1996). Glycolysis occurs in both the cytosol and plastids in plants (Plaxton, 1996). Triose phosphates generated by pTPI are central to carbon metabolism in plants, serving as the photosynthesis-derived carbohydrate skeletons for anabolic synthesis of more complex carbohydrates such as the hexoses, disaccharides, and starch, as well as the glycolytic intermediates for catabolic energy production (Plaxton, 1996;Raines, 2003). These processes are critical for the carbon-and energy-intensive demands of plant defenses. Changes in pTPI gene expression have been found in response to fungal inoculation in cabbage (Brassica carinata; Subramanian, Bansal, & Kay, 2005), and beech trees (Fagus grandifolia Ehrh; Mason, Koch, Krasowski, & Loo, 2013). We next examined whether these genes were implicated in lodgepole and jack pine defense responses by comparing transcript abundance under control and induced conditions. Given the logistical challenges associated with conducting gene expression analyses on MPB attacked pines, we used G. clavigera inoculation as a proxy to elicit an MPB defense response. We used a novel Bayesian mixed-effect model (Matz et al., 2013) to analyze the transcript abundance data. These analyses revealed that transcript abundance corresponding to proteasome subunit α 7 significantly decreased in lodgepole pine-but not in jack pine-shortly after fungal inoculation. These data suggest that this gene, harboring the Lodgc1087 attack-associated SNP, may play some role in host susceptibility. In contrast, transcript abundance corresponding to pTPI showed a significant increase in response to both wounding and fungal inoculation at early time points in jack pine, but not lodgepole pine. As this experiment was not explicitly designed to test gene expression differences in Lodgc1087 genetic variants, we are unable to establish a link between gene expression and the mutation. Ongoing experiments will address this question in both seedling and mature tree material.

| CON CLUS IONS
In this study, we have used complementary validation appropriate for nonmodel organisms by incorporating evidence from multiple approaches to identify genetic variation potentially contributing to . Given we validated only two candidates that were previously identified using stringent criteria, should serve as a precaution for studies attempting to identify adaptive candidates without proper validation as a consideration (Meirmans, 2015).
Of the two validated candidates, one shows promise as a marker for increased susceptibility to MPB attack. We have identified a clear association between attack status and genotype for the Lodgc1087containing gene encoding proteasome subunit α 7 in lodgepole pine, and provided experimental evidence that this gene is involved in the defense response. This is the first report of a locus associated with pine host susceptibility, of which genetic resistance is one component. Genetic resistance to MPB has been reported to be moderately heritable in lodgepole pine (Yanchuk, Murphy, & Wallin, 2008), but resistance to this herbivore likely comprises hundreds-if not thousands-of loci with potentially small effects. Genotype-by-environment effects are also likely to be substantive. Nevertheless, our finding of an association of Lodg1087 with attack status demonstrates the utility of a population genomics approach to defining the genetic component of pine host suitability for MPB, and is a meaningful first step toward delineating a susceptible vs. resilient host genotype. This approach has been adopted as part of a comprehensive, large-scale effort to identify additional loci associated with attack status that are also implicated in pine defense responses.
Our goal is to develop models that will better define the linkages between host genotype and attack probability. Such models can be used for applications such as refining stand susceptibility indices used in decision support systems (Cooke & Carroll, 2017), as well as in tree improvement and ex situ conservation programs to develop more resilient forests.

ACK N OWLED G EM ENTS
The authors would like to recognize the important contributions of Sophie Dang, Stephanie Boychuk, Dominik Royko, Joël Fillion, and Bill Clark (University of Alberta) for their support and assistance in data cataloguing, collection, and manipulation, Adriana Arango-Velez for the glasshouse experiment data and samples, and the staff members at Delta Genomics for SNP genotyping and support.
Sample collection was thanks to many helping hands on the TRIA field crews. We would also like to thank Sebastian Lackey (University of Alberta) for logistics. We acknowledge funding for this research to J.E.K.C and D.W.C. from Genome Canada, the Government of Alberta through Genome Alberta, Genome British Columbia in support of the Tria I and Tria II Projects (http://www.thetriaproject.  Table S1).