Genome-Wide Association Study of Wood Anatomical and Morphological Traits in Populus trichocarpa

To understand the genetic mechanisms underlying wood anatomical and morphological traits in Populus trichocarpa, we used 869 unrelated genotypes from a common garden in Clatskanie, Oregon that were previously collected from across the distribution range in western North America. Using GEMMA mixed model analysis, we tested for the association of 25 phenotypic traits and nine multitrait combinations with 6.741 million SNPs covering the entire genome. Broad-sense trait heritabilities ranged from 0.117 to 0.477. Most traits were significantly correlated with geoclimatic variables suggesting a role of climate and geography in shaping the variation of this species. Fifty-seven SNPs from single trait GWAS and 11 SNPs from multitrait GWAS passed an FDR threshold of 0.05, leading to the identification of eight and seven nearby candidate genes, respectively. The percentage of phenotypic variance explained (PVE) by the significant SNPs for both single and multitrait GWAS ranged from 0.01% to 6.18%. To further evaluate the potential roles of candidate genes, we used a multi-omic network containing five additional data sets, including leaf and wood metabolite GWAS layers and coexpression and comethylation networks. We also performed a functional enrichment analysis on coexpression nearest neighbors for each gene model identified by the wood anatomical and morphological trait GWAS analyses. Genes affecting cell wall composition and transport related genes were enriched in wood anatomy and stomatal density trait networks. Signaling and metabolism related genes were also common in networks for stomatal density. For leaf morphology traits (leaf dry and wet weight) the networks were significantly enriched for GO terms related to photosynthetic processes as well as cellular homeostasis. The identified genes provide further insights into the genetic control of these traits, which are important determinants of the suitability and sustainability of improved genotypes for lignocellulosic biofuel production.

Introduction

It is of increasing interest to identify the molecular variants underlying adaptive and morphological trait variation in plant populations. Loci highlighted by such analyses have great potential for optimizing the trait of interest through genetic engineering or breeding, thereby producing trees with increased productivity, enhanced abiotic stress tolerance, and/or improved quality of end products. Because of their wide geographic distribution and climatic gradients, large effective population sizes, and high genetic variation, forest trees are excellent model systems for understanding local adaptation and the genetic architecture of complex traits (Neale and Savolainen, 2004; González-Martínez et al., 2006; Neale and Kremer, 2011; Street and Ingvarsson, 2011; Ingvarsson et al., 2016). In this regard, efforts have been made to optimize the ecologically and economically important tree Populus for lignocellulosic biofuel production. Vast amounts of genomic and phenotypic resources are available for the genus. Several large-scale genome-wide association studies have identified the underlying genetic architecture related to morphological, physiological, wood property and chemistry, salinity tolerance, and disease resistance traits (Ma et al., 2013; McKown et al., 2014b; Muchero et al., 2015; Zhang et al., 2018; Bdeir et al., 2019; Quan et al., 2019; Jia et al., 2020; Lu et al., 2020). Furthermore, the biology of wood formation, cell wall ultrastructure and composition, and cell wall recalcitrance are fairly well studied (Groover et al., 2010; Wegrzyn et al., 2010; Studer et al., 2011; Du et al., 2013; Porth et al., 2013a; Porth et al., 2013b; Muchero et al., 2015; Porth et al., 2015; Allwright et al., 2016; Du et al., 2016; Escamez et al., 2017; Fahrenkrog et al., 2017; Johnson et al., 2017; Xi et al., 2017; Du et al., 2018; Gandla et al., 2018; Du et al., 2019) and a few studies have also recognized the role of microRNA in controlling tree growth and wood property traits in Populus (Quan et al., 2016; Chen B. et al., 2018). However, the genetic architecture underlying wood anatomical traits such as vessel size and density are relatively uncharacterized, despite the importance of these traits for cell wall composition and the overall performance of the tree.

Wood anatomy not only contributes to the structural integrity of the tree, but it is also critical for transport and storage processes (Sperry, 2003; Hietz et al., 2017). Anatomical structures like vessel size and density are related to cell wall structure and composition and wood density. These traits together affect long-distance axial transport of nutrients and hydraulic conductivity. Lignified cell walls and fibers add strength and living parenchyma cells provide radial transport and storage. There are often trade-offs among vessel properties, wood density and hydraulic conductivity (Preston et al., 2006). Wood traits are generally heritable (Carlquist, 2012) and serve as useful traits in phylogenetic analyses (Hietz et al., 2017). Radial variation in wood anatomical properties in the stem affects wood functional traits such as hydraulic conductivity and often scales with tree size and leaf characteristics (Lachenbruch et al., 2011). Genome-wide association studies of wood anatomy traits undoubtedly complement the current understanding of genetic architecture of other structural and functional traits in Populus and will serve as an important step for optimizing traits suitable for lignocellulosic biofuels.

Nevertheless, the limitation of power due to sample size in detecting the genetic variants associated with complex traits is a major hurdle in trees and other non-model organisms. The recent technological revolution in genomics has made comprehensive genotyping routine, but challenges remain due to the cost and logistical difficulties inherent in the establishment, management, and intensive phenotyping of large common gardens of trees. Furthermore, SNP loci and candidate genes identified thus far explain only a small proportion of the genetic variation in complex traits in general (Visscher et al., 2017). Genes controlling complex traits do not work in isolation, but instead are interconnected in networks of hundreds to thousands of genes, each of which may contribute incrementally to the variation in complex traits (Boyle et al., 2017). Recently, methods such as multitrait GWAS and meta analyses using summary statistics (e.g., principal components) have become increasingly popular due to their role in enhancing the power of GWAS and the identification of potentially pleiotropic loci (Porter and O’Reilly, 2017). Performing multitrait GWAS directly on raw phenotypic data can enhance power and lead to discovery of novel associations (Chhetri et al., 2019).

Multi-omics network methods can be used to boost the signal of true positives and reduce noise in GWAS results by demonstrating concurrent support for hypotheses from multiple lines of evidence (LOE) (Pendergrass et al., 2013; Verma et al., 2018; Weighill et al., 2019). This can be particularly helpful for studies in which sample size or power are limiting factors (Genovese et al., 2006). Furthermore, LOE methods help identify GWAS hits that are strong gene candidates for experimental validation by examining the gene across multiple contexts in parallel. By considering expression patterns in conjunction with phenotypic associations, genes that have a higher probability of being feasibly validated in lab, greenhouse, and common garden settings can be prioritized.

Here we focus on Populus trichocarpa, a targeted species for lignocellulosic biofuel production that has a distribution spanning from northern California to northern British Columbia. Tremendous resources including whole-genome resequencing data, multiple common gardens for association mapping, transcriptome and metabolite data and expression networks are available for this species (Chhetri et al., 2019; Weighill et al., 2019). We present a genome-wide association study of wood anatomical traits for the first time for this species. We also present a GWAS for important morphological traits from the same trees that together with anatomical traits affect overall plant productivity. We complemented the single trait GWAS with multitrait analyses, followed by multi-omic LOE and functional enrichment analyses. Since this study is based on data collected from a common garden in Clatskanie, Oregon this allowed for the direct comparison of GWAS genes identified for some of the same traits from other plantations (Evans et al., 2014; Chhetri et al., 2019).

Materials and Methods

Phenotypic Data Collection

Wood anatomical, and morphological trait data were collected from a field trial of 1,100 P. trichocarpa genotypes that was established in Clatskanie, Oregon in 2009 (Figure 1). These genotypes were previously collected from across the natural range of P. trichocarpa from northern California to northern British Columbia and were clonally replicated and planted in a randomized block design with three replicates of each genotype at 2 m x 3 m spacing in the field trial (Slavov et al., 2012; Evans et al., 2014). Phenotyping was performed on clonal replicates for all traits as described in Table 1. In some cases replicate observations were recorded for a subset of the genotypes, with the number of replicated genotypes ranging from 59 to 860.

FIGURE 1

Figure 1 Source locations of 869 P. trichocarpa genotypes sampled in this study (purple dots). The trees were grown in a common garden in Clatskanie, Oregon, USA (black star).

TABLE 1

Table 1 Broad-sense heritability estimates (H²) and the number of SNP-trait associations for wood anatomical and morphological traits in P. trichocarpa.

Wood Anatomy

In June 2012, 557 trees were sampled for wood anatomical traits. Wood cores of 5.5 mm diameter were extracted from bark to pith from the southern face of the main trunk of the tree using an increment borer. Free hand cross-sections were made from the previous year’s growth ring and fixed in 70% alcohol. The tissue sections were stained in 1% Safranin O solution for 30 s before preparing the slides for imaging. Images at 100x magnification were taken to sample the early, intermediate and late wood from the growth ring. These images were used for measuring vessel density, size and number (Figure 2). All images were processed using the software imageJ to extract quantitative measurements—vessel count, density and size (Schindelin et al., 2015).

FIGURE 2

Figure 2 Populus trichocarpa wood cross section from the previous year’s growth ring showing variation in early, intermediate and late wood types.

Leaf Anatomy

Leaf characteristics were measured for 676–687 trees (Table 1). The first and second fully expanded leaf (counting from the apex) were collected from a branch receiving full sunlight. One of the leaves was used for measuring petiole length and minimum and maximum petiole diameter with a caliper and then scanned using a hand-held scanner. The software imageJ (Schindelin et al., 2015) was used to estimate leaf area, leaf circularity, leaf length, leaf width, and leaf perimeter. Dry and wet weights were measured for the same leaf and leaf area and leaf dry weight were used to estimate the specific leaf area (SLA). Leaf chlorophyll content (SPAD) was assessed using a SPAD 502 Plus meter (Spectrum Technologies) with an average of three replicate measures on a leaf section. The second leaf was used for measuring the abaxial stomatal density. Clear nail polish was applied to the broadest part of the leaf close to the midrib. A clear piece of tape was then used to capture an imprint of the epidermal leaf surface and stored on microscope slides. The number of stomata in 1 mm² area in four random microscopic fields at 400x magnification were counted.

Leaf Metabolites

Leaf material was collected from 851 P. trichocarpa genotypes. Metabolites were extracted from tissues, measured by GCMS, and peaks were characterized (Weighill et al., 2018). For each sample, a Median Absolute Deviation (MAD) threshold of six was used to remove phenotypic values that were more than six MADs from the population median for that phenotype. In addition, phenotypes with greater than 20% missing data were removed prior to analysis, resulting in a total of 816 metabolite phenotypes. A Benjamini-Hochberg false-discovery rate cutoff of 0.1 was applied to the adjusted p-values.

Wood Composition

The wood cores were ground, processed and analyzed by pyrolysis molecular beam mass spectrometry (py-MBMS) as previously described (Muchero et al., 2015; Weighill et al., 2018). py-MBMS data were TIC normalized to account for variation in sample mass. Phenotypes with > 20% missing data or zeroes were removed prior to analysis (one phenotype, “m/z 76”). For each sample, a Median Absolute Deviation (MAD) threshold of five was used to remove phenotype values that were more than five MADs from the population median for that phenotype, resulting in the removal of an average of 1.49 py-MBMS phenotype values per sample (in total, 625 phenotype values).

Statistical Analyses

To estimate the genetic control of quantitative traits, broad-sense heritability (H²) was estimated for all traits using the genotypes with replicate clonal measurements using the following formula:

$$H^2=\frac{{\sigma }_G^2}{{\sigma }_G^2+{\sigma }_E^2},$$

where, ${\sigma }_G^2$ is genotypic variance due to clonal differences and ${\sigma }_E^2$ is environmental variance.

Outliers were removed, and the data were evaluated for normality. Variance components were estimated employing the linear regression model with the lmer and ranef functions of the lme4 package implemented in R. Genotype and the position of the tree (i.e. row and column) in the garden was used as a random effect in the model. Error was estimated from the residuals of the model. Genetic correlation between the traits was estimated using the Best Linear Unbiased Predictors (BLUPs) from the same model (Supplementary Figure S1). Phenotypic BLUPs were also used for estimating the correlation between the phenotypic traits and 26 geoclimate variables and their first four principal components of the source location of the trees. Genetic correlations were also calculated between phenotypic traits collected at this site to those collected in a clonally replicated trial with the same genotypes in Corvallis, Oregon (Chhetri et al., 2019).

Genotypic Data

Methods for obtaining genotypic data were as described previously (Weighill et al., 2018; Chhetri et al., 2019). Briefly, whole genome re-sequencing data was obtained from 1,053 trees using Illumina genetic analyzers at the DOE Joint Genome Institute. After removing trees related more closely than first cousins and highly differentiated California trees, 869 trees were left, which were used for all analyses. A genetic relationship matrix was estimated using GEMMA and used as a covariate in the GWAS analyses. Furthermore, SNPs with minor allele frequency ≤0.05 were removed.

Association Analysis

The tests for statistical association for all phenotypic traits were conducted using the Genome-wide Efficient Mixed Model Association package (GEMMA, Zhou & Stephens, 2012; Zhou & Stephens, 2014), except for the py-MBMS data, which were analyzed using the efficient mixed-model association eXpedited (EMMAX) package (Kang et al., 2010). Phenotypic BLUPs, a genetic relationship matrix and 6,741,160 SNPs were used for the association test. Single trait GWAS was run for 25 phenotypes (Table 1). The tested model was:

$$y=x\beta +u+є,$$

where y is an n-vector of phenotypic BLUP values, where n is the number of individuals tested; x is an n-vector of marker genotypes, β is the effect size of the marker, u is an n-vector of random effects that includes a relatedness matrix and є is an n-vector of errors. We also tested Principal Components as covariates in the model, as described previously (Evans et al., 2014), but we determined that inclusion of PCs resulted in overcorrection, as determined by a preponderance of points below the 1:1 line in Q-Q plots (not shown).

Multitrait GWAS was run for 9 sets of phenotype combinations (Table 2). The same procedure was used for selecting trait sets as described in Chhetri et al. (2019). Briefly, traits were combined if there was reason to believe that they could be functionally related. For example, leaf traits were combined with wood anatomy traits due to the expected interactive roles of both sets of traits in water stress tolerance (Sperry, 2000). Multitrait association was conducted with GEMMA using the same model as for single trait associations, except y is an n×d matrix of d phenotypes for n individuals.

TABLE 2

Table 2 List of traits used for multitrait associations and significant SNPs identified in P. trichocarpa.

We used an FDR cutoff of 0.05 and a more liberal FDR cutoff of 0.1 to identify suggestive associations (Storey and Tibshirani, 2003). For the purpose of summarizing the results, significant SNPs within 10 kb of one another were merged into peaks. Gene models that were closest to significant SNPs were identified based on v3 of the P. trichocarpa genome. Annotation information including gene expression level in different plant tissues and annotation of putative gene function was obtained from Phytozome (Goodstein et al., 2012). Percentage of variance explained (PVE) by SNPs was estimated using the formula in Shim et al. (2015).

Multiple Lines of Evidence (LOE) Analysis

A network was created using five data layers generated from P. trichocarpa genetic and phenotypic analyses: pairwise gene coexpression, pairwise gene comethylation, metabolite GWAS, rare variant regional metabolite GWAS, and py-MBMS GWAS. The first four of these layers were previously described (Furches et al., 2019). The py-MBMS GWAS analysis was performed for this study. For each layer, a network was constructed that consisted of gene pairs (nodes) linked by a correlation value or association p-value (edge).

Coexpression Network

P. trichocarpa (Nisqually-1 reference) RNA-seq data was obtained from the DOE Joint Genome Institute Plant Gene Atlas (available from Phytozome at: https://phytozome.jgi.doe.gov/pz/portal.html). For each sample, expression was measured in leaf, stem, root and bud tissues at multiple developmental stages and four different conditions (three nitrogen treatments and an untreated control). Plant cultivation, tissue collection, and RNA extraction and sequencing were previously described in Weighill et al. (2018), as were RNA-seq read trimming, read alignment, TPM calculations, Spearman correlation analysis, and network construction. A correlation value threshold of 0.85 was set as criteria for retention of gene pairs in the coexpression network. As previously demonstrated, this threshold minimized false positives, and the resulting network was significantly different from random (Furches et al., 2019).

Comethylation Network

P. trichocarpa MeDIP-seq data (Vining et al., 2012) was obtained from Phytozome. This data set contained MeDIP-seq reads from ten tissues: bud, callus, male and female catkins, internode, regenerated internode, leaf, phloem, xylem, and root. Read mapping, calculation of transcripts per million reads, Spearman correlation analysis, and network construction were conducted as previously described (Weighill et al., 2018). Gene pairs with correlation values greater than or equal to 0.95 were retained in the comethylation network (Furches et al., 2019).

Rare Variant Regional Metabolite GWAS Networks

Leaf metabolite data were analyzed using a rare variant regional GWAS analysis as described elsewhere (Weighill et al., 2018). This analysis jointly tested rare SNPs (MAF ≤ 0.01) within each gene’s boundary and the 2-kb up- and downstream flanking regions of that gene using a linear mixed model as implemented by SKAT in RVtest (Zhan et al., 2016). An FDR cutoff of 0.1 was applied to the resulting p-values.

LOE Scoring of Candidate Genes

The purpose of this LOE method was to explore the networks of genes underlying the traits analyzed in the primary leaf and wood anatomy GWAS analyses to assess support for the functional hypotheses suggested by the gene annotations. For each of the significant (FDR<0.1) gene models identified by the single trait and multitrait GWAS analyses of wood and leaf anatomy traits (hereafter referred to as “anchor genes”), a merged network of anchor genes and their nearest one-hop neighbors in six data layers (including the gene-trait associations detected in this study) was created. Network statistics were calculated, including breadth (the number of layers to which a gene is connected), depth per layer (the number of nodes to which the gene is connected), and total depth. Anchor genes were directly scored based on the number of shared edges in each layer rather than scoring the nearest neighbors of anchor genes. These scores were then used to prioritize genes for further analysis.

Functional Enrichment of Coexpression Nearest Neighbors

The purpose of this analysis is to test whether functional categories represented in the one-hop networks occur more frequently than expected from a random draw of the same number of genes from the genome. For each anchor gene, two separate Gene Ontology (GO) functional enrichment analyses were performed on positively and negatively coexpressed one-hop neighbors, respectively using the software package Fuento (Weichselbaum et al., 2017). Fisher’s Exact tests were performed for all available namespaces (molecular function, cellular compartment, biological process, aberrant functions) using 100 tries and an FDR adjusted p-value cutoff of p ≤ 0.1.

Results

Heritabilities of Wood Anatomy and Morphology Traits

Overall the broad-sense heritabilities estimated for wood anatomy and morphology related traits in this study were low to moderate (Table 1). Broad-sense heritabilities for wood anatomy traits ranged from 0.114 for late wood vessel area to 0.392 for early wood vessel size. Late wood anatomy (late wood vessel count, vessel area and vessel size) had low heritabilities compared to other wood traits. For morphology related traits, broad-sense heritabilities ranged from 0.122 for specific leaf area to 0.477 for stomatal density.

Heritability estimates derived from the SNP data (“chip heritabilities”) for wood anatomy traits ranged from 0.199 for early wood vessel size to 0.611 for early wood vessel count whereas for morphology related traits it ranged from 0.038 for specific leaf area to 0.965 for tree height (Table 1). Although there was no correlation between broad-sense and chip heritabilities, low chip heritability was generally reflective of low broad-sense heritability and there was a significant correlation between the sample size (number of genotypes used for GWAS) and chip heritability values (r=0.441, P=0.021) for the traits measured in this study (Table 1). GWAS results for traits with low chip heritabilities should be treated with caution.

Genetic Correlation of Phenotypic Traits Within and Between the Common Gardens

As expected, most morphological traits measured in the Clatskanie common garden were correlated to each other—leaf traits such as leaf area, leaf length, leaf dry weight, leaf wet weight, petiole diameter and length, leaf perimeter, and leaf aspect ratio were highly significantly correlated to each other and SPAD and stomatal density had low, but significant correlations with most leaf traits (Supplementary Figure S2). Similarly, most wood traits were correlated to each other. Wood anatomy traits such as vessel area and size and vessel count within each of the growth ring regions (early, intermediate or late wood indicating a different maturation stages) were correlated to each other. Wood anatomy traits were also significantly correlated to each other across growth ring areas (Supplementary Figure S2). Furthermore, some morphological traits such as tree height and diameter were significantly correlated with wood traits such as, early wood vessel size and area, intermediate wood vessel size and count, and late wood vessel size and count (Supplementary Figure S2).

Genetic correlation of the same phenotypic traits measured between the Clatskanie and Corvallis common gardens showed that most traits were weak to moderately, but significantly correlated to each other (Supplementary Figure S3). Pairwise genetic correlations of tree height (r=0.436, p<0.001) and stomatal density (0.345, p<0.001) between the sites showed the strongest correlations among all traits compared (Supplementary Figure S3).

Phenotypic Trait Correlations With Climate Variables

As seen for the phenotypic traits in the Corvallis common garden (Chhetri et al., 2019), most phenotypic traits had significant correlations with latitude, and therefore the correlations between the phenotypic traits and other geoclimate variables cannot easily be disentangled. Nevertheless, most morphological traits had significant (although weak) correlations with most geoclimatic variables (Supplementary Figure S4). Similarly, wood traits such as late wood vessel area had significant (although weak) correlations with most geoclimatic variables (Supplementary Figure S4).

Genes Identified From Single Trait and Multitrait GWAS

We performed single trait GWAS with 6.741 million SNPs for 25 morphological and wood anatomical traits. Fifty-seven SNPs passed the FDR threshold of p<0.05 (Table 1). However, we identified a total of 112 SNPs that passed suggestive association FDR P-value cutoff of 0.1 (Table 1, Figures 3 and 4, Supplementary Figures S5 and S6, Supplementary Table S1). These associated SNPs belonged to 12 separate SNP peaks, and 21 P. trichocarpa gene models were identified as nearest neighbors to the significant SNPs (Table 3, Supplementary Table S2). Percentage of variance explained (PVE) for significant SNPs ranged from 3.30% to 5.65% (Supplementary Table S1).

FIGURE 3

Figure 3 Single trait GWAS – Manhattan (left) and QQ plots (right). Numbers 1 to 19 represent chromosomes; scaffolds are the reads that did not align to any of the 19 chromosomes. (A, B) Leaf dry weight; (C, D) Leaf wet weight. SNPs above gray line have p<1×10^-7, which is roughly equivalent to FDR ≤ 0.1 and SNPs above black line (the top line) have p<7.417×10^-9 (the Bonferroni correction threshold), which is roughly equivalent to FDR ≤ 0.05 in this study.

FIGURE 4

Figure 4 Single trait GWAS – Manhattan (left) and QQ plots (right). Numbers 1 to 19 represent chromosomes; scaffolds are the reads that did not align to any of the 19 chromosomes. (A, B) Stomatal density; (C, D) Vessel size (intermediate wood). SNPs above gray line have p<1×10^-7, which is roughly equivalent to FDR ≤ 0.1 and SNPs above black line (the top line) have p<7.417×10^-9 (the Bonferroni correction threshold), which is roughly equivalent to FDR ≤ 0.05 in this study.

TABLE 3

Table 3 Genes identified from P. trichocarpa single trait GWAS.

Multitrait GWAS performed on 9 sets of traits identified 11 SNPs that passed the FDR p-value cutoff of 0.05 and 19 SNPs that passed the suggestive association FDR p-value cutoff of 0.1. (Table 2, Figure 5, Supplementary Figure S7, Supplementary Table S3). These SNPs belonged to 13 separate SNP peaks and were within or close to 13 P. trichocarpa gene models (Table 4, Supplementary Table S4). PVE of these SNPs ranged from 0.01% to 6.18% for the individual traits comprising the multitrait set (Supplementary Table S3).

FIGURE 5

Figure 5 Multitrait GWAS—Manhattan (left) and QQ plots (right). Numbers 1 to 19 represent chromosomes; scaffolds are the reads that did not align to any of the 19 chromosomes. (A, B) Leaf area, leaf dry weight, leaf length and leaf width; (C, D) Vessel area and vessel count (late wood). SNPs above gray line have P<1×10^-7, which is roughly equivalent to FDR ≤ 0.1 and SNPs above black line (the top line) have P<7.417×10^-9 (the Bonferroni correction threshold), which is roughly equivalent to FDR ≤ 0.05 in this study.

TABLE 4

Table 4 Genes identified from P. trichocarpa multitrait GWAS.

Out of the single trait GWAS for nine wood anatomical traits, we identified only one gene model belonging to intermediate wood vessel size (Figures 4C, D). A total of three gene models for seven wood anatomy multitrait sets (Figures 5C, D) and three gene models for one morphology and wood anatomy multitrait set (Figures S6A, B) were identified.

LOE Scores and Rankings

Breadth scores ranged from one to four (Supplementary Table S5, Supplementary Figure S8). Two anchor genes had breadth scores of four: Potri.019G067300 and Potri.010G031900. Both of these genes were associated with the multitrait GWAS for leaf morphology (LA_LD_LL_LW). They have connections to coexpressed and comethylated genes, metabolite, and py-MBMS phenotypes (Supplementary Tables S6 and S7). Two anchor genes had breadth scores of three, both of which were identified through single trait GWAS analysis: Potri.019G037900 (wood anatomy), and Potri.014G159800 (morphology). Total depth scores ranged from one to 2,014. The ranges of layer specific depth scores were: zero to 1,115 for coexpression, zero to 1,914 for comethylation, zero to one for metabolite GWAS, zero to one for RV metabolite GWAS, and zero to one for py-MBMS GWAS. A total of 7630 unique coexpressed and comethylated genes, four metabolites and three py-MBMS GWAS were found to be associated with 33 GWAS (defined as anchor) genes for single and multitrait GWAS (Figures 6–8, Supplementary Tables S6 and S7, Supplementary Figures S9–S13).

FIGURE 6

Figure 6 One-hop multi-omic network for the gene Potri.019G042600, which was associated with leaf wet weight.

FIGURE 7

Figure 7 One-hop multi-omic network for the gene Potri.015G117300, which was associated with leaf abaxial stomatal density.

FIGURE 8

Figure 8 One-hop multi-omic network for the gene Potri.019G067300, which was associated with the multitrait phenotype leaf area, leaf dry weight, leaf length, and leaf width (LA_LD_LL_LW).

Functional Enrichment Results

Thirteen anchor genes had significantly enriched functions among positively coexpressed one-hop neighbors. However, after FDR correction, only five anchor genes had enriched neighborhoods (Supplementary Table S8). Seven anchor genes had significantly enriched functions among negatively coexpressed neighbors, but only a single enrichment for one anchor gene remained after FDR correction. Overall there was an abundance of signaling, transport, growth and development, RNA regulation, hormone metabolism, cell wall related and stress related genes for the coexpressed and comethylated genes in the network for all anchor genes (Supplementary Figure S14).

Shared LOE for Single Trait and Multi-Trait Associations

Only one gene, Potri.015G113800, was significant for both a multitrait set (early wood vessel area-leaf area-stomatal density) and for one of the traits making up the set (stomatal density). We also searched for overlaps between the one-hop networks for multi-traits and the individual traits, but again this was the only case where there were shared genes in the one-hop networks (Supplemental Figure S15).

Discussion

Patterns of Genetic Variation

Broad-sense heritabilities for most traits in this study ranged from 0.2 to 0.4, suggesting that the traits were under moderate genetic control. Heritability estimates for morphological traits were comparable to similar studies (Evans et al., 2014; Mckown et al., 2014; Chhetri et al., 2019). However, some wood related traits and SLA had low heritabilities (<0.2), which may be due to relatively small sample size compared to morphological traits. Low heritability may also be due to micro-environmental variation in the common garden, which can have a strong effect on the growth ring patterns and the anatomical traits studied here. Furthermore, some imprecision in wood trait phenotypes may result from the high throughput preparation method that we were forced to adopt due to the large sample numbers. Wood samples were prepared via free-hand sectioning, which results in less uniform images than can be obtained with embedded sections prepared with a microtome. In any case, the results for these low heritability traits should be interpreted with caution.

As reported in previous studies geography plays a major role in shaping adaptive trait variation in P. trichocarpa. Correlation of these adaptive traits with latitude can obscure the actual between-trait and trait-climate relationships (Chhetri et al., 2019). Most morphological traits such as tree height and leaf characteristics had strong significant correlations with latitude compared to previous studies (Mckown et al., 2014; Chhetri et al., 2019). In contrast to Chhetri et al. (2019), tree height including diameter at breast height had a very low or no correlation with leaf traits as in McKown et al. (2014c). However, tree height and diameter had significant negative correlations with intermediate and late wood vessel counts and significant positive correlations with vessel size of all maturation stages (early, intermediate and late wood). Interestingly, vessel density was positively correlated with latitude and vessel size was negatively correlated with latitude. It is therefore difficult to deconvolute the indirect effects of latitudinal variation from the direct effects of wood anatomy on productivity in this study.

Further corroboration for the effects of environment shaping these traits were evident when comparing data from the two common gardens. Genetic correlation of the directly comparable traits between the Clatskanie and Corvallis common garden showed very low or no correlations except for stomatal density and tree height (Supplementary Figure S3). Furthermore, different loci seem to be controlling phenotypic traits in the two gardens, based on a complete lack of overlap in significant loci discovered by GWAS in the Corvallis, Oregon common garden (Chhetri et al., 2019). This might be due to differential effects of local environments on phenotypes in each of the respective common gardens, or may simply indicate that our analysis is underpowered.

As temperature decreases, and moisture and precipitation increase with latitude, relationships of the phenotypic traits with climate variables corroborates the trait relationships explained above. Leaf characteristics including leaf dry weight and petiole length and diameter had significant positive correlation with latitude, mean annual precipitation and relative humidity, but significant negative correlation with mean annual temperature. Similarly, abaxial stomatal density was negatively correlated with precipitation (Supplementary Figure S4). Taken together, southern trees had larger but less dense vessels, smaller leaves with higher abaxial stomatal density and low chlorophyll content, and larger tree height and diameter. Late wood vessel area (similar to the relationship of height and diameter with latitude and temperature) decreased with latitude and increases with temperature.

GWAS Genes

The GWAS analysis for wood anatomical traits we reported here is the first such study in Populus. While we found none to only a few significantly associated or suggestive SNPs for most traits we studied here, we think that this is mainly due to lack of power due to sample size although this is one of the most comprehensive GWAS analyses to date for P. trichocarpa in terms of the number of SNPs used. We did not find any gene model controlling more than one trait from our single trait GWAS analyses, but we identified 13 gene models with potential pleiotropic effects controlling sets of morphological and wood anatomical traits (Table 4).

Single Trait GWAS Genes

For single trait GWAS, we identified eight genes belonging to four traits based on 57 SNPs that passed the FDR threshold of 0.05 (Table 3). The top five genes identified were for two morphological traits, leaf wet weight (three genes) and stomatal density (two genes).

The gene model identified based on the top GWAS hit for both leaf wet weight and leaf dry weight (Figure 3) was Potri.019G042600 which is a mitochondrial outer membrane protein based on PANTHER annotation (https://phytozome.jgi.doe.gov/pz/portal.html). The gene is 6.1 kb downstream of the GWAS SNP (FDR = 0.015, p-value = 7.18 × 10^-9). Outer membrane proteins in mitochondria (and chloroplasts) play a role in important cellular functions such as protein import and ion exchange (Lee et al., 2014). The Arabidopsis homolog of this gene, AT5G17170 encodes rubredoxin family protein located in the plant chloroplast and is involved in electron transport and metal ion binding activity in Arabidopsis (https://www.arabidopsis.org/). Expression of rubredoxin family protein in a grass, Puccinellia tenuiflora (PutRUB) was shown to be induced by abiotic stresses due to NaCl, NaHCO₃, CuCl₂, and H₂O₂ suggesting a role of PutRUB in stress tolerance. Overexpression of PutRUB in Arabidopsis plants under stressed conditions showed an increase in seedling weight in overexpressed plants compared to wild type (Li et al., 2016). The GWAS gene has a breadth score of 2 and a total depth score of 1,116 (3^rd highest among all single and multitrait anchor genes, Figure 6, Supplementary Table S5). It is highly expressed in late winter apical bud and fully open bud, stem, young and first fully expanded leaf and stem treated with ammonia, nitrate and urea (Supplementary Table S7). The anchor gene is coexpressed with 1,115 other genes (Figure 6) of which 915 are positively and 200 are negatively coexpressed. The coexpressed genes have functions related to photosystem I and II in photosynthesis (155 genes), oxidation reduction reaction (26), RNA regulation and binding (62), signaling (53, including light, G-protein, calcium, receptor kinases and sugar and nutrient physiology), abiotic and biotic stresses (25), transport (36), glycolysis (11), and hormone metabolism (25). All of the top 10 coexpressed genes are related to the photosynthesis electron transport chain including Photosystem I and Photosystem II complexes (Supplementary Table S6). Similarly, Potri.019G042600 showed significant functional enrichment for GO term Biological Processes (photosynthesis and respiration related) with functional annotations including ATP synthesis, oxidation reduction, glycolysis and photosynthesis. The gene also showed enrichment for GO term Molecular Function that included hydrolase, oxidoreductase, calcium ion binding, antioxidant activity and proton transport (Supplementary Table S8). In sum, network and functional enrichment analyses suggest that the function of this GWAS gene is related to photosynthetic and respiration processes and responses to biotic and abiotic stresses. One of the possible mechanisms by which plants respond to abiotic stresses, especially nutrient and water stress, is by changing morphology (e.g, leaf area and weight) and physiological processes (e.g., respiration) (Rymaszewski et al., 2017), thereby raising the expectation for coordinated responses of genes related to photosynthesis, oxidation reduction reactions, glycolysis, signaling, hormone metabolism, and stress responses, along with those controlling leaf morphology.

Potri.015G117500 was the top gene associated with stomatal density (p-value = 2.77 × 10^-9, Figures 4A, B). The top GWAS SNP is a synonymous substitution harbored within the gene. Potri.015G117500 encodes peptide-O-fucosyltransferase, which belongs to a family of genes involved in post-translational modification (Lira-Navarrete et al., 2011). The gene had a breadth score of 1 and the total depth score of 17 (comethylated with 17 other genes, Supplementary Table S5). Notable comethylated genes include a HSP20-like chaperone super family protein (Potri.009G153200), photosystem II oxygen evolving PsbQ complex subunit protein (Potri.011G031300) and farnesylcysteine lyase (Potri.007G101700) that indirectly has a negative regulatory effect in the ABA mediated signaling pathway in Arabidopsis (Huizinga et al., 2010). Another coexpressed gene, Potri.019G006900 encodes kunitz trypsin inhibitor and is involved in apoptosis during plant-pathogen interactions and responding to salt stress and salicylic acid. Kunitz trypsin inhibitors are known to have a function in pest resistance in poplar (Major and Peter Constabel, 2008). Similarly, another comethylated gene, Potri.001G045500, encodes cinnamoyl CoA reductase which is involved in lignin biosynthesis (Lacombe et al., 1997; Kawasaki et al., 2006). Lignin biosynthetic genes are known to affect growth and defense mechanisms in plants (Xie et al., 2018).

In Arabidopsis the O-fucosyltransferase SPINDLY (SPY, similar to the protein encoded by Potri.015G117500) affects DELLA and associated regulators such as brassinazole-resistant1 (BZR1) phytochrome-interacting-factor3 (PIF3) and PIF4 in the brassinosteroid- and light-signaling pathways (Zentella et al., 2017). DELLA proteins are important integrators of multiple signaling pathways in flowering plants. DELLAs are repressors of phytohormone (GA) signaling, master growth repressors that restrict plant growth by affecting genes in cell division, expansion and differentiation (Zentella et al., 2007). In Arabidopsis O-fucosyltransferase modifies the DELLA protein RGA (repressor of ga1-3). The GA-deficient mutant ga1-3 is a nongerminating, extreme dwarf that flowers late and produces male-sterile flowers. This protein is also known to be involved in cell-to-cell adhesion (deficiency in pectin biosynthesis pathway leading to loss of cell adhesion) together with GAUT8 (QUASIMODO1, see below) in Arabidopsis (Verger et al., 2016). McKown et al. (2014a) identified a similar gene, Potri.004G059000, encoding BRASSINOSTEROID-INSENSITIVE 2 for stomatal density based on the GWAS using a 34K SNP array data for 464 P. trichocarpa trees. BRASSINOSTEROID-INSENSITIVE 2 potentially affects the brassinosteroid signaling pathway and the regulation of stomatal development (Gudesblat et al., 2012; Khan et al., 2013). Brassinosteroids also play important role in plant growth, development and stress responses (Nolan et al., 2020).

The other top GWAS hit for stomatal density was within a gene model Potri.015G117300, a gene of unknown function (Figures 4A, B). The anchor gene has a breadth score of 2 and a total depth score of 2014 (highest breadth score among all single and multitrait anchor genes, Figure 7, Supplemental Table S5). The gene is comethylated and coexpressed with 1914 and 100 other genes, respectively, that encode proteins related to hormone signaling, stress response and transport functions. The anchor gene is coexpressed or comethylated with the genes related to diverse functions including photosynthesis (22), cell wall (44), development (84), glycolysis (6), hormone metabolism (56), UDP glucosyl and glucoronyl transferases (6), redox (12), signaling (153, including receptor kinase, MAP kinase, light, sugar, G-proteins and calcium), stress (64, including abiotic and biotic, PR-proteins, respiratory burst), transport (123). For example, Potri.015G117300 is comethylated with two other genes, Potri.001G139700 and Potri.005G109300, both of which encode for plasma membrane proteins that belong to the subfamily of plant aquaporins (AQPs) that play important roles in plant physiological processes (Chaumont and Tyerman, 2014). Membrane proteins have been shown to have an involvement in regulating plant-water homeostasis and abiotic stress responses (Sade et al., 2009). Similarly, the anchor gene is comethylated with Potri.005G123200, that encodes for Sec14p-like phosphatidylinositol transfer family protein that has important role in lipid signaling pathways in plants (Huang et al., 2016). Another comethylated gene, Potri.005G134900, encodes a LATERAL ORGAN BOUNDARIES (LOB) domain protein. Proteins in this family play a role in organ development and stress responses (Grimplet et al., 2017), and also in regulating the expression of VND7, a master regulator of tracheary elements in plants (Ohashi-Ito et al., 2018). Overexpression of VND7 represses KNOX gene activity, which in turn is responsible for meristem and organ development including the development of leaf primordia in plants (Hake et al., 2004) The anchor gene is coexpressed with germin-like (GLP) protein 1 (Potri.013G141900). GLPs are involved in a wide range of functions including growth and development, and biotic and abiotic stress responses (Dunwell et al., 2008). Another coexpressed gene, Potri.009G140900, encodes Small Auxin-Up RNA (SAUR)-like protein is involved in plant growth and development (cell division in particular) and the expression of the gene is regulated by auxin (Markakis et al., 2013). The coexpressed gene, Potri.010G112800 encodes auxin efflux carrier family protein, PIN-FORMED 3 (PIN3) which plays a role in seedling elongation and shade avoidance in Arabidopsis (Keuskamp et al., 2010). Another coexpressed gene, Potri.011G067500, encodes CBL-interacting protein kinase which is involved in drought tolerance mechanisms and the regulation of stomatal aperture (Wang et al., 2016).

Stomatal density in leaves facilitates gas exchange, and variation in the density of stomata can have a significant effect on overall plant fitness as it can change the dynamics of photosynthetic rate, respiration and growth (Franks and Farquhar, 2007; Franks and Beerling, 2009). Our GWAS and network analyses highlight the importance of genes involved in hormonal signaling and stress responses. Other recent studies in Populus have suggested a relationship between stomatal patterning and the trade-off between growth and disease susceptibility (and presence of defense compounds) (McKown et al., 2014a; McKown et al., 2019).

Multitrait GWAS Genes

We identified a total of seven genes based on the 11 SNPs that passed an FDR threshold of 0.05 in the multitrait GWAS, of which 6 genes were detected for the leaf morphology multitrait, (leaf area, leaf dry weight, leaf length and leaf width) and 1 gene was detected for wood anatomy multitrait (late wood vessel area and vessel count) (Figure 5)

The first anchor gene for the leaf morphology multitrait set was Potri.019G067300 (breadth score 4, total depth score 395, FDR=0.021, p-value=1.9 x 10^-8; Figures 5A, B and 8, Supplementary Table S5). The gene encodes a protein related to galactosyltransferase-8 (GAUT8). The gene is comethylated with 392 other genes, positively coexpressed with one other gene, and related to one metabolite and linked to a py-MBMS (wood chemistry) GWAS. The anchor gene is comethylated with genes with various functions including cell wall (11 genes), growth and development related (11), glycolysis and Krebs cycle related (6), photosystem I and II related (4), oxidation reduction related (6), RNA regulation & binding (37), signaling (31, includes receptor kinase, light, photosynthesis, G-proteins and calcium signaling), stress (16, abiotic such as drought/salt and heat, biotic such as PR-proteins, respiratory burst), and transport (15) (Figure 8). The coexpressed gene Potri.006G049600 encodes concanavalin A-like lectin protein kinase family protein that is involved in amino acid phosphorylation. Lectin receptor-like kinases (LecRLK) are known to regulate abiotic and biotic stresses, and development in plants (Vaid et al., 2013; Yang et al., 2016; Labbe et al., 2019). A comethylated gene with a similar function as the anchor genes is Potri.006G137100, which encodes a pectin lyase-like superfamily protein that functions in pectin methylesterase (PME) activity. PMEs are involved in regulation of plant growth and development via structural modeling of the cell wall and also play an important role in heat stress tolerance (Wu et al., 2018). Another gene with a similar function, Potri.004G111000 (GAUT9) was detected for the same multitrait GWAS from the Corvallis common garden (Chhetri et al., 2019). GAUT8 and GAUT9 belonged to the same clade B1 in a phylogenetic tree based on the GAUT protein family of Arabidopsis thaliana and P. trichocarpa (Biswal et al., 2018). GAUT9 affects leaf size in Populus deltoides (Chhetri et al., 2019) whereas GAUT8 (also known as QUASIMODO1) has high expression in stems in Arabidopsis (Caffall et al., 2009) and is thought to be involved in cell wall pectic homogalacturonan (HG) and xylan biosynthesis (Orfila et al., 2005; Verger et al., 2016). Genes encoding glycosyltransferases, a large family of enzymes mainly involved in biosynthesis of polysaccharides and glycoproteins in the plant cell wall (Hansen et al., 2012), were also detected in single trait GWAS for stomatal density in this study (Table 3).

Another gene associated with the leaf morphology multi-trait set was Potri.010G031900 (FDR = 0.043, P-value = 5.23 x 10^-8) is coexpressed and comethylated with 12 and 67 other genes, respectively (breadth score 4, total depth score 81, Figures 5A, B, Supplementary Table S5). The gene encodes E1 C-terminal protein, which is a RUB-activating enzyme that indirectly triggers auxin response by interacting with the cullin AtCUL1 in Arabidopsis (Del Pozo and Estelle, 1999). The anchor gene is coexpressed or comethylated with genes with diverse functions including signaling (9, including G-proteins, MAP kinase, receptor kinase and calcium), RNA regulation (9), transport (3), stress (5, abiotic- drought/salt and heat, biotic—PR-proteins), hormone metabolism (4, auxin 2, ABA 2), and lipid metabolism (2) (Supplementary Figure S13).

A few of the top coexpressed genes in the network are related to the RING/U-box superfamily protein (e.g., Potri.018G064400). U-box (PUB) proteins have diverse functions in plants, including hormonal, defense and abiotic stress responses (Yee and Goring, 2009). RNA binding protein (Potri.019G017900) and mitogen-activated protein kinase kinase 3 (MKK3, Potri.001G345500)—MKK3 belong to a bigger group of mitogen-activated protein kinases (MAPKs) that play a role in activating downstream hormonal signaling including auxin, jasmonic acid, salicylic acid, brassinosteroid, ethylene and ABA signaling targets in plants (Jagodzik et al., 2018). The gene is also comethylated with genes with similar roles such as Potri.009G015800 encoding IQ-domain 5 (IQD5) that functions in calmodulin binding. Calmodulin and calmodulin-like proteins are involved in responding to abiotic stresses in plants (Zeng et al., 2015; Badmi et al., 2018). There are other comethylated genes in the network with functions in stress response, and growth and development. Potri.009G001200 encodes RING-H2 group F2A (RHF2A) protein. RING-H2 belong to a group of E3 ubiquitin protein complexes and are involved in defense response, regulation of growth and development, and apoptosis (Guzmán, 2012). Potri.008G010800 encodes ABA-responsive element binding protein 3. Potri.002G152100 encodes phospholipase D beta 1 (PLDBETA1). PLD functions in regulating plant growth and development, hormonal signaling and stress responses (Li et al., 2007). The anchor gene is also associated with a metabolite, caffeoyl-quercetin glycoside and one py-MBMS GWAS, m/z 414). Quercetin glycosides are known to play a role in osmotic adjustment during drought stress (Tschaplinski et al., 2019). Caffeoyl-quercetin glycoside was found to be associated with Potri.001G411800 (EF-hand Calcium-Binding Domain protein) in a multitrait GWAS for leaf area, stomatal density and carbon isotope composition in a previous study (Chhetri et al., 2019). Ethylene and ABA play important roles in stress responses with prominent effects in plant tissues (Cramer et al., 2011). Reduction in Leaf water potential and thus leaf wet weight is possible due to abiotic (e.g. drought) and biotic stress (e.g. pathogen) (Pandey et al., 2017). Drost et al. (2015) provided evidence for the effect of PtARF1 on leaf morphology in a hybrid population of P. deltoides and P. trichocarpa that is mediated by auxin signaling such that the PIN polarization in leaf cells possibly plays a role in differential expression of PtARF1, ultimately causing the expansion of leaf in the length direction with the decrease in leaf width in P. trichocarpa and vice versa in P. deltoides (Kalluri et al., 2007). Various signaling genes including G-protein signaling (Lease et al., 2001), and hormonal signaling genes (McSteen and Zhao, 2008) are known to have an effect on plant growth and development. Calcium signaling regulates the function of stress related genes ultimately affecting the plant phenotype (Tuteja and Mahajan, 2007).

Multitrait GWAS Genes: Wood Anatomy

The gene model Potri.001G058100 was associated with the wood anatomy multitrait (Figure 5). This gene has no known function but is coexpressed and comethylated with 9 and 972 other genes (breadth score 2, total depth score 981, 4^th highest depth score among all single and multitrait GWAS genes), respectively (Supplementary Tables S5 and S9). The gene is comethylated with 3 RNA binding genes (Potri.009G152900, Potri.013G094400, Potri.017G063700) that are involved in biogenesis of small RNA in Arabidopsis (Eamens et al., 2009; Eamens et al., 2012). Small RNAs are known to regulate the expression of genes related to growth and development, and wood formation in plants (Tang et al., 2016; Chen C. et al., 2018). The anchor gene is coexpressed or comethylated with diverse genes that have functions in transport (61 genes), biotic and abiotic stress tolerance (35), signaling mechanism (84, hormonal, receptor kinase, light and calcium, sugar and G-protein signaling), cell wall related (20), developmental mechanisms (39), hormone metabolism (25) and RNA regulation (106) (Figures 4C, D; Supplementary Table S9). Notable top 10 coexpressed genes include Potri.001G043100 that encodes Chaperone DnaJ-domain superfamily protein which is a heat shock protein with a function in stress tolerance, Potri.002G169400 (encodes RCAR1, regulatory components of ABA receptor), Potri.009G033300 (encodes SCARECROW-like 14, SCL14) and Potri.017G134800 (encoding zinc finger protein 7, ZFP7). SCL14 is a member of GRAS family of transcription factors that interacts with the TGA2 transcription factors and affects the transcription of stress-responsive genes (Fode et al., 2008). ZFP7 is a member of C2H2 zinc finger protein related to ZFP3. ZFP3 and related ZFP proteins including ZFP7 acts as a negative regulator of ABA signaling during seedling development, affect vegetative growth and regulate light signaling in Arabidopsis (Joseph et al., 2014). ABA receptors are involved in transmitting ABA signals that allow plants to respond to drought and other abiotic stresses (Fan et al., 2016; Duarte et al., 2019). The anchor gene is comethylated with Leucine-rich repeat protein kinase family protein (Potri.019G001800) that functions in ATP binding activity and regulating hormonal signaling and stress response in plants (Máthé et al., 2019). Taken together, this suggests a role of this gene in growth and development including wood formation by maintaining cellular homeostasis, signaling and stress response.

Comparison With Previous Studies

We compared our GWAS results with previous GWAS studies for similar traits in P. trichocarpa (Wegrzyn et al., 2010; Porth et al., 2013a; Porth et al., 2013b; Evans et al., 2014; McKown et al., 2014b; Chhetri et al., 2019) and other Populus species including P. deltoides (Fahrenkrog et al., 2017), P. euphratica (Ma et al., 2013; Jia et al., 2020) and P. tomentosa (Quan et al., 2016; Chen B. et al., 2018; Quan et al., 2019; Lu et al., 2020). We found very minimal overlap in the genes and/or gene families detected. Ma et al. (2013) identified multiple gene families that control adaptation in saline environments. Similar to their results, heat shock protein families and some ABA signaling proteins were found in the co-expression networks in our study. Jia et al. (2020) identified 82 genes related to seed salinity tolerance in P. euphratica, including a zinc finger protein similar to one that we detected for leaf abaxial density that regulates gas exchange in plants. Lu et al. (2020) identified 23 genes involved in three growth and six wood property traits and Quan et al. (2019) identified 203 lignin biosynthetic genes, 81 TF genes, 36 microRNAs and 71 long non-coding RNA genes for 10 growth and wood property traits in the same P. tomentosa population. However, these studies have no apparent overlap in gene functions for similar morphology or wood property traits in our study.

Comparing our findings with previous P. trichocarpa studies we identified some non-overlapping gene models with similar functions in shared gene families. For example, we identified a gene Potri.019G067300 in a multi-trait GWAS for leaf morphology that encodes the GAUT8 protein. In a previous multi-trait GWAS study of leaf morphology at a different site, we identified GAUT9 (Potri.004G111000), which affects leaf size and is involved in pectin and xylan biosynthesis (Chhetri et al., 2019). Likewise, we identified a gene (Potri.005G247700) encoding “SKU5 similar 5” (sks5) that functions in oxidoreductase activity and copper ion binding in this study. Porth et al. (2013a) also identified a gene with similar function for soluble lignin (Potri.001G000500) encoding “SKU5 similar 12” (sks12) based on 34K SNP genotyping array data. McKown et al. (2019) identified several genes related to plant growth, physiological processes, biotic and abiotic stresses, and defense and immunity for abaxial and adaxial stomata related traits (stomatal density and pore length). We did not find any of the genes they identified in their study, despite analyzing similar traits here. The discrepancies in the gene models identified across the studies may be due to a) differential environmental effects on the traits due to plantation sites, b) traits measured in different seasons of the year and different years were subjected to different environmental conditions, c) variability in the microenvironment within the sites as reflected by the different broad-sense heritability values, d) differences in the developmental stages of the traits across the studies, e) difference in the genotypes and the sequence data types (e.g. SNP array vs whole genome sequence data) and f) the lack of power in this and previous studies. Furthermore, we were not able to replicate any significant SNP hits or the gene models for the same traits from another P. trichocarpa common garden in Corvallis, Oregon in our study, although several single traits such as stomatal density, petiole diameter and tree height and the multitrait set containing leaf area, leaf dry weight, leaf length and leaf width were common between the two common gardens, and measurements were performed using the same methods. This might very well reflect the differential influence of the environment on genotypes in the two populations. Nevertheless, we believe that our comprehensive GWAS study with genome-wide sequencing data highlights the underlying genes controlling adaptive traits in P. trichocarpa and complements the findings from other similar studies. Moreover, we have reported here potential genes underlying complex wood anatomy traits in Populus for the first time.

Conclusion

We presented here the first comprehensive GWAS for wood anatomical traits in Populus that provides insights into the type of genes controlling structural and functional properties important for plant development, function and stress tolerance. We complemented this with GWAS for important morphological traits that have functional relationships with wood anatomical traits. As shown in a previous P. trichocarpa GWAS study in another common garden (Chhetri et al., 2019), we identified additional genes with multitrait GWAS in this study as well. Furthermore, the multitrait sets formed based on the genetic correlations and functional relationships of the traits within and among the wood anatomical and morphological traits provided insight into pleiotropic genes controlling these traits. Some of the genes we identified in this study had no known functions. However, the LOE network analyses we performed here provided additional information not only for suggesting the role of the unknown genes, but also added an extra layer of support for the genes with known functions. Genes identified here, especially the genes controlling wood anatomical traits, can be good targets for biotechnological approaches toward optimizing wood traits for biofuel production. However, as is the case with most GWAS studies, despite the use of large scale whole genome resequencing data and much broader sampling compared to most previous studies in trees, a very small percent of the variation in the traits were explained by the significant SNPs. It is possible that we may simply be underpowered to detect SNPs with small effect sizes, given that we only detected a handful of genes or none for most of the complex traits that might potentially be controlled by hundreds or thousands of genes (Boyle et al., 2017). A more robust approach to GWAS is required, with large and homogeneous sampling from across the range of distribution. Alternatively, more populations with controlled structure can be used to take full advantage of the range of variation present while explicitly controlling for underlying structure, as is the case in nested association mapping (Yu et al., 2008). Moreover, none of the genetic variants detected in the GWAS from another common garden were validated in this study. This may be because of the differential effects of the environment between the common gardens, which is somewhat reflected by very low correlation of the traits between the two common gardens. Furthermore, additional variation due to the difference in the timing (different years) of data collection might have influenced the GWAS study. Nevertheless, by integrating data from multiple sources we were able to apply a LOE approach that provided further corroboration for relatively weak associations, while providing insights into gene function.

Data Availability Statement

Populus trichocarpa genome sequence, annotation, and Gene Atlas expression and methylation data sets are available on Phytozome (http://phytozome.jgi.doe.gov). Populus trichocarpa variant data (DOI 10.13139/OLCF/1411410) is available from https://doi.ccs.ornl.gov/ui/doi/55.

Author Contributions

HC, AF, DM-S, AW, DK, PJ, and DJ performed statistical and computational analyses. AH-W and TT performed metabolite analyses and analyzed the data. DJ, GT, and SD oversaw the study. HC, AF, and SD drafted the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Center for Bioenergy Innovation (CBI) and the Bioenergy Science Center. CBI is supported by the Office of Biological and Environmental Research in the DOE Office of Science. This manuscript has been coauthored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). The work conducted by the U.S. Department of Energy Joint Genome Institute is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research described herein was supported by an award of computer time provided by the INCITE program and used resources of the Oak Ridge Leadership Computing Facility (OLCF) at the Oak Ridge National Laboratory. This work was authored in part by Alliance for Sustainable Energy, LLC, the manager and operator of the National Renewable Energy Laboratory for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding was provided by U.S. Department of Energy Office of Energy Efficiency and Renewable Energy Bioenergy Technologies Office.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the multitude of researchers from the Bioenergy Science Center and the DOE Joint Genome Institute who provided invaluable logistical support for this work. In particular, we would like to thank Kat Haiby, Brian Stanton, Rich Shuren, Carlos Gantz, and Austin Himes of Greenwood Resources for their work in establishing and maintaining the plantation, for facilitating our work at the site, and for the many insights that they have provided about Populus biology and silviculture. We would also like to thank Crissa Doeppke and Robert Sykes for their help with py-MBMS analysis.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2020.545748/full#supplementary-material

References

Allwright, M. R., Payne, A., Emiliani, G., Milner, S., Viger, M., Rouse, F., et al. (2016). Biomass traits and candidate genes for bioenergy revealed through association genetics in coppiced European Populus nigra (L.). Biotechnol. Biofuels 9, 1–22. doi: 10.1186/s13068-016-0603-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Badmi, R., Payyavula, R. S., Bali, G., Guo, H. B., Jawdy, S. S., Gunter, L. E., et al. (2018). A new calmodulin-binding protein expresses in the context of secondary cell wall biosynthesis and impacts biomass properties in Populus. Front. Plant Sci. 871:1669. doi: 10.3389/fpls.2018.01669

CrossRef Full Text | Google Scholar

Bdeir, R., Muchero, W., Yordanov, Y., Tuskan, G. A., Busov, V., Gailing, O. (2019). Genome-wide association studies of bark texture in Populus trichocarpa. Tree Genet. Genomes 15, 14. doi: 10.1007/s11295-019-1320-2

CrossRef Full Text | Google Scholar

Biswal, A. K., Atmodjo, M. A., Pattathil, S., Amos, R. A., Yang, X., Winkeler, K., et al. (2018). Working towards recalcitrance mechanisms: Increased xylan and homogalacturonan production by overexpression of GAlactUronosylTransferase12 (GAUT12) causes increased recalcitrance and decreased growth in Populus. Biotechnol. Biofuels 11, 1–26. doi: 10.1186/s13068-017-1002-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyle, E. A., Li, Y., II, Pritchard, J. K. (2017). An Expanded View of Complex Traits: From Polygenic to Omnigenic. Cell 169, 1177–1186. doi: 10.1016/j.cell.2017.05.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Caffall, K. H., Pattathil, S., Phillips, S. E., Hahn, M. G., Mohnen, D. (2009). Arabidopsis thaliana T-DNA mutants implicate GAUT genes in the biosynthesis of pectin and xylan in cell walls and seed testa. Mol. Plant 2, 1000–1014. doi: 10.1093/mp/ssp062

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlquist, S. (2012). How wood evolves: a new synthesis. Botany 90, 901–940. doi: 10.1139/b2012-048

CrossRef Full Text | Google Scholar

Chaumont, F., Tyerman, S. D. (2014). Aquaporins: Highly regulated channels controlling plant water relations. Plant Physiol. 164, 1600–1618. doi: 10.1104/pp.113.233791

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Chen, J., Du, Q., Zhou, D., Wang, L., Xie, J., et al. (2018). Genetic variants in microRNA biogenesis genes as novel indicators for secondary growth in Populus. New Phytol. 219, 1263–1282. doi: 10.1111/nph.15262

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Zeng, Z., Liu, Z., Xia, R. (2018). Small RNAs, emerging regulators critical for the development of horticultural traits. Hortic. Res. 5, 6–8. doi: 10.1038/s41438-018-0072-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Chhetri, H. B., Macaya-sanz, D., Kainer, D., Biswal, A. K., Evans, L. M., Chen, J., et al. (2019). Multitrait genome-wide association analysis of Populus trichocarpa identifies key polymorphisms controlling morphological and physiological traits. New Phytol. 223, 293–309. doi: 10.1111/nph.15777

PubMed Abstract | CrossRef Full Text | Google Scholar

Cramer, G. R., Urano, K., Delrot, S., Pezzotti, M., Shinozaki, K. (2011). Effects of abiotic stress on plants: A systems biology perspective. BMC Plant Biol. 11, 163. doi: 10.1186/1471-2229-11-163

PubMed Abstract | CrossRef Full Text | Google Scholar

Del Pozo, J. C., Estelle, M. (1999). The Arabidopsis cullin AtCUL1 is modified by the ubiquitin-related protein RUB1. Proc. Natl. Acad. Sci. U. S. A. 96, 15342–15347. doi: 10.1073/pnas.96.26.15342

PubMed Abstract | CrossRef Full Text | Google Scholar

Drost, D. R., Puranik, S., Novaes, E., Novaes, C. R. D. B., Dervinis, C., Gailing, O., et al. (2015). Genetical genomics of Populus leaf shape variation. BMC Plant Biol. 15, 1–10. doi: 10.1186/s12870-015-0557-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Q., Pan, W., Xu, B., Li, B., Zhang, D. (2013). Polymorphic simple sequence repeat (SSR) loci within cellulose synthase (PtoCesA) genes are associated with growth and wood properties in Populus tomentosa. New Phytol. 197, 763–776. doi: 10.1111/nph.12072

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Q., Gong, C., Wang, Q., Zhou, D., Yang, H., Pan, W., et al. (2016). Genetic architecture of growth traits in Populus revealed by integrated quantitative trait locus (QTL) analysis and association studies. New Phytol. 209, 1067–1082. doi: 10.1111/nph.13695

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Q., Lu, W., Quan, M., Xiao, L., Song, F., Li, P., et al. (2018). Genome-Wide Association Studies to Improve Wood Properties: Challenges and Prospects. Front. Plant Sci. 9, 1912. doi: 10.3389/fpls.2018.01912

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, Q., Yang, X., Xie, J., Quan, M., Xiao, L., Lu, W., et al. (2019). Time-specific and pleiotropic quantitative trait loci coordinately modulate stem growth in Populus. Plant Biotechnol. J. 17, 608–624. doi: 10.1111/pbi.13002

PubMed Abstract | CrossRef Full Text | Google Scholar

Duarte, K. E., de Souza, W. R., Santiago, T. R., Sampaio, B. L., Ribeiro, A. P., Cotta, M. G., et al. (2019). Identification and characterization of core abscisic acid (ABA) signaling components and their gene expression profile in response to abiotic stresses in Setaria viridis. Sci. Rep. 9, 1–16. doi: 10.1038/s41598-019-40623-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunwell, J. M., Gibbings, J. G., Mahmood, T., Saqlan Naqvi, S. M. (2008). Germin and germin-like proteins: Evolution, structure, and function. CRC. Crit. Rev. Plant Sci. 27, 342–375. doi: 10.1080/07352680802333938

CrossRef Full Text | Google Scholar

Eamens, A. L., Smith, N. A., Curtin, S. J., Wang, M. B., Waterhouse, P. M. (2009). The Arabidopsis thaliana double-stranded RNA binding protein DRB1 directs guide strand selection from microRNA duplexes. Rna 15, 2219–2235. doi: 10.1261/rna.1646909

PubMed Abstract | CrossRef Full Text | Google Scholar

Eamens, A. L., Kim, K. W., Curtin, S. J., Waterhouse, P. M. (2012). DRB2 is required for microRNA biogenesis in Arabidopsis thaliana. PloS One 7, e35933. doi: 10.1371/journal.pone.0035933

PubMed Abstract | CrossRef Full Text | Google Scholar

Escamez, S., Latha Gandla, M., Derba-Maceluch, M., Lundqvist, S. O., Mellerowicz, E. J., Jönsson, L. J., et al. (2017). A collection of genetically engineered Populus trees reveals wood biomass traits that predict glucose yield from enzymatic hydrolysis. Sci. Rep. 7, 1–11. doi: 10.1038/s41598-017-16013-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, L. M., Slavov, G. T., Rodgers-Melnick, E., Martin, J., Ranjan, P., Muchero, W., et al. (2014). Population genomics of Populus trichocarpa identifies signatures of selection and adaptive trait associations. Nat. Genet. 46, 1089–1096. doi: 10.1038/ng.3075

PubMed Abstract | CrossRef Full Text | Google Scholar

Fahrenkrog, A. M., Neves, L. G., Resende, M. F. R., Dervinis, C., Davenport, R., Barbazuk, W. B., et al. (2017). Population genomics of the eastern cottonwood (Populus deltoides). Ecol. Evol. 7, 9426–9440. doi: 10.1002/ece3.3466

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, W., Zhao, M., Li, S., Bai, X., Li, J., Meng, H., et al. (2016). Contrasting transcriptional responses of PYR1/PYL/RCAR ABA receptors to ABA or dehydration stress between maize seedling leaves and roots. BMC Plant Biol. 16, 1–14. doi: 10.1186/s12870-016-0764-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fode, B., Siemsen, T., Thurow, C., Weigel, R., Gatz, C. (2008). The arabidopsis GRAS protein SCL14 interacts with class II TGA transcription factors and is essential for the activation of stress-inducible promoters. Plant Cell 20, 3122–3135. doi: 10.1105/tpc.108.058974

PubMed Abstract | CrossRef Full Text | Google Scholar

Franks, P. J., Beerling, D. J. (2009). Maximum leaf conductance driven by CO2 effects on stomatal size and density over geologic time. Proc. Natl. Acad. Sci. U. S. A. 106, 10343–10347. doi: 10.1073/pnas.0904209106

PubMed Abstract | CrossRef Full Text | Google Scholar

Franks, P. J., Farquhar, G. D. (2007). The mechanical diversity of stomata and its significance in gas-exchange control. Plant Physiol. 143, 78–87. doi: 10.1104/pp.106.089367

PubMed Abstract | CrossRef Full Text | Google Scholar

Furches, A., Kainer, D., Weighill, D., Large, A., Jones, P., Walker, A. M., et al. (2019). Finding New Cell Wall Regulatory Genes in Populus trichocarpa Using Multiple Lines of Evidence. Front. Plant Sci. 10, 1249. doi: 10.3389/fpls.2019.01249

PubMed Abstract | CrossRef Full Text | Google Scholar

Gandla, M. L., Martín, C., Jönsson, L. J. (2018). Analytical Enzymatic Saccharification of Lignocellulosic Biomass for Conversion to Biofuels and Bio-Based Chemicals. Energies 11, 2936. doi: 10.3390/en11112936

CrossRef Full Text | Google Scholar

Genovese, C. R., Roeder, K., Wasserman, L. (2006). False discovery control with p-value weighting. Biometrika. 93, 509–524. doi: 10.1093/biomet/93.3.509

CrossRef Full Text | Google Scholar

González-Martínez, S. C., Ersoz, E., Brown, G. R., Wheeler, N. C., Neale, D. B. (2006). DNA sequence variation and selection of tag single-nucleotide polymorphisms at candidate genes for drought-stress response in Pinus taeda L. Genetics 172, 1915–1926. doi: 10.1534/genetics.105.047126

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodstein, D. M., Shu, S., Howson, R., Neupane, R., Hayes, R. D., Fazo, J., et al. (2012). Phytozome: A comparative platform for green plant genomics. Nucleic Acids Res. 40, 1178–1186. doi: 10.1093/nar/gkr944

CrossRef Full Text | Google Scholar

Grimplet, J., Pimentel, D., Agudelo-Romero, P., Martinez-Zapater, J. M., Fortes, A. M. (2017). The LATERAL ORGAN BOUNDARIES Domain gene family in grapevine: Genome-wide characterization and expression analyses during developmental processes and stress responses. Sci. Rep. 7, 1–18. doi: 10.1038/s41598-017-16240-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Groover, A. T., Nieminen, K., Helariutta, Y., Mansfield, S. D. (2010). “Wood Formation in Populus,” in Genetics and Genomics of Populus. Eds. Jansson, S., Bhalerao, R., Groover, A. T. (New York: Springer), 201–224.

Google Scholar

Gudesblat, G. E., Schneider-Pizoń, J., Betti, C., Mayerhofer, J., Vanhoutte, I., Van Dongen, W., et al. (2012). SPEECHLESS integrates brassinosteroid and stomata signalling pathways. Nat. Cell Biol. 14, 548–554. doi: 10.1038/ncb2471

PubMed Abstract | CrossRef Full Text | Google Scholar

Guzmán, P. (2012). The prolific ATL family of RING-H2 ubiquitin ligases. Plant Signal. Behav. 7, 1014–1021. doi: 10.4161/psb.20851

PubMed Abstract | CrossRef Full Text | Google Scholar

Hake, S., Smith, H. M. S., Holtan, H., Magnani, E., Mele, G., Ramirez, J. (2004). The Role of Knox Genes in Plant Development. Annu. Rev. Cell Dev. Biol. 20, 125–151. doi: 10.1146/annurev.cellbio.20.031803.093824

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, S. F., Harholt, J., Oikawa, A., Scheller, H. V. (2012). Plant Glycosyltransferases Beyond CAZy: A Perspective on DUF Families. Front. Plant Sci. 3, 59. doi: 10.3389/fpls.2012.00059

PubMed Abstract | CrossRef Full Text | Google Scholar

Hietz, P., Rosner, S., Hietz-Seifert, U., Wright, S. J. (2017). Wood traits related to size and life history of trees in a Panamanian rainforest. New Phytol. 213, 170–180. doi: 10.1111/nph.14123

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Ghosh, R., Bankaitis, V. A. (2016). Sec14-like phosphatidylinositol transfer proteins and the biological landscape of phosphoinositide signaling in plants. Biochim. Biophys. Acta - Mol. Cell Biol. Lipids 1861, 1352–1364. doi: 10.1016/j.bbalip.2016.03.027

CrossRef Full Text | Google Scholar

Huizinga, D. H., Denton, R., Koehler, K. G., Tomasello, A., Wood, L., Sen, S. E., et al. (2010). Farnesylcysteine lyase is involved in negative regulation of abscisic acid signaling in Arabidopsis. Mol. Plant 3, 143–155. doi: 10.1093/mp/ssp091

PubMed Abstract | CrossRef Full Text | Google Scholar

Ingvarsson, P. K., Hvidsten, T. R., Street, N. R. (2016). Towards integration of population and comparative genomics in forest trees. New Phytol. 212, 338–344. doi: 10.1111/nph.14153

PubMed Abstract | CrossRef Full Text | Google Scholar

Jagodzik, P., Tajdel-Zielinska, M., Ciesla, A., Marczak, M., Ludwikow, A. (2018). Mitogen-activated protein kinase cascades in plant hormone signaling. Front. Plant Sci. 9, 1387. doi: 10.3389/fpls.2018.01387

PubMed Abstract | CrossRef Full Text | Google Scholar

Jia, H., Liu, G., Li, J., Zhang, J., Sun, P., Zhao, S., et al. (2020). Genome resequencing reveals demographic history and genetic architecture of seed salinity tolerance in Populus euphratica. J. Exp. Bot. 71, 4308–4320. doi: 10.1093/jxb/eraa172

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, A. M., Kim, H., Ralph, J., Mansfield, S. D. (2017). Natural acetylation impacts carbohydrate recovery during deconstruction of Populus trichocarpa wood. Biotechnol. Biofuels 10, 1–12. doi: 10.1186/s13068-017-0734-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Joseph, M. P., Papdi, C., Kozma-Bognár, L., Nagy, I., López-Carbonell, M., Rigó, G., et al. (2014). The Arabidopsis ZINC FINGER PROTEIN3 interferes with abscisic acid and light signaling in seed germination and plant development. Plant Physiol. 165, 1203–1220. doi: 10.1104/pp.113.234294

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalluri, U. C., Difazio, S. P., Brunner, A. M., Tuskan, G. A. (2007). Genome-wide analysis of Aux/IAA and ARF gene families in Populus trichocarpa. BMC Plant Biol. 7, 1–14. doi: 10.1186/1471-2229-7-59

PubMed Abstract | CrossRef Full Text | Google Scholar

Kang, H. M., Sul, J. H., Service, S. K., Zaitlen, N. A., Kong, S. Y., Freimer, N. B., et al. (2010). Variance component model to account for sample structure in genome-wide association studies. Nat. Genet. 42, 348–354. doi: 10.1038/ng.548

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawasaki, T., Koita, H., Nakatsubo, T., Hasegawa, K., Wakabayashi, K., Takahashi, H., et al. (2006). Cinnamoyl-CoA reductase, a key in lignin biosynthesis, is an effector of small GTPase Rac in defense signaling in rice. Proc. Natl. Acad. Sci. U. S. A. 103, 230–235. doi: 10.1073/pnas.0509875103

PubMed Abstract | CrossRef Full Text | Google Scholar

Keuskamp, D. H., Pollmann, S., Voesenek, L. A. C. J., Peeters, A. J. M., Pierik, R. (2010). Auxin transport through PIN-FORMED 3 (PIN3) controls shade avoidance and fitness during competition. Proc. Natl. Acad. Sci. U. S. A. 107, 22740–22744. doi: 10.1073/pnas.1013457108

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, M., Rozhon, W., Bigeard, J., Pflieger, D., Husar, S., Pitzschke, A., et al. (2013). Brassinosteroid-regulated GSK3/Shaggy-like kinases phosphorylate mitogen-activated protein (MAP) kinase kinases, which control stomata development in Arabidopsis thaliana. J. Biol. Chem. 288, 7519–7527. doi: 10.1074/jbc.M112.384453

PubMed Abstract | CrossRef Full Text | Google Scholar

Labbe, J., Muchero, W., Czarnecki, O., Wang, J., Wang, X., Bryan, A. C., et al. (2019). Mediation of plant-mycorrhizal interaction by a lectin receptor-like kinase. Nat. Plants 5, 676–680. doi: 10.1038/s41477-019-0469-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lachenbruch, B., Moore, J. R., Evans, R. (2011). “Radial Variation in Wood Structure and Function in Woody Plants, and Hypotheses for its Occurence,” in Size and Age Related Changes in Tree Structure and Function. Eds. Meinzer, F. C., Lachenbruch, B., Dawson, T. E. (Dordrecht: Springer), 121–164.

Google Scholar

Lacombe, E., Hawkins, S., Van Doorsselaere, J., Piquemal, J., Goffner, D., Poeydomenge, O., et al. (1997). Cinnamoyl CoA reductase, the first committed enzyme of the lignin branch biosynthetic pathway: Cloning, expression and phylogenetic relationships. Plant J. 11, 429–441. doi: 10.1046/j.1365-313X.1997.11030429.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lease, K. A., Wen, J., Li, J., Doke, J. T., Liscum, E., Walker, J. C. (2001). A mutant Arabidopsis heterotrimeric G-protein β subunit affects leaf, flower, and fruit development. Plant Cell 13, 2631–2641. doi: 10.1105/tpc.13.12.2631

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J., Kim, D. H., Hwang, I. (2014). Specific targeting of proteins to outer envelope membranes of endosymbiotic organelles, chloroplasts, and mitochondria. Front. Plant Sci. 5, 173. doi: 10.3389/fpls.2014.00173

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G., Lin, F., Xue, H. W. (2007). Genome-wide analysis of the phospholipase D family in Oryza sativa and functional characterization of PLDB1 in seed germination. Cell Res. 17, 881–894. doi: 10.1038/cr.2007.77

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Liu, P., Takano, T., Liu, S. (2016). A chloroplast-localized rubredoxin family protein gene from Puccinellia tenuiflora (PutRUB) increases NaCl and NaHCO3 tolerance by decreasing H2O2 accumulation. Int. J. Mol. Sci. 17, 804. doi: 10.3390/ijms17060804

CrossRef Full Text | Google Scholar

Lira-Navarrete, E., Valero-González, J., Villanueva, R., Martínez-Júlvez, M., Tejero, T., Merino, P., et al. (2011). Structural insights into the mechanism of protein O-fucosylation. PloS One 6, e25365. doi: 10.1371/journal.pone.0025365

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, W., Xiao, L., Quan, M., Wang, Q., El-Kassaby, Y. A., Du, Q., et al. (2020). Linkage-linkage disequilibrium dissection of the epigenetic quantitative trait loci (epiQTLs) underlying growth and wood properties in Populus. New Phytol. 225, 1218–1233. doi: 10.1111/nph.16220

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, T., Wang, J., Zhou, G., Yue, Z., Hu, Q., Chen, Y., et al. (2013). Genomic insights into salt adaptation in a desert poplar. Nat. Commun. 4, 2797. doi: 10.1038/ncomms3797

PubMed Abstract | CrossRef Full Text | Google Scholar

Major, I. T., Peter Constabel, C. (2008). Functional analysis of the kunitz trypsin inhibitor family in poplar reveals biochemical diversity and multiplicity in defense against herbivores. Plant Physiol. 146, 888–903. doi: 10.1104/pp.107.106229

PubMed Abstract | CrossRef Full Text | Google Scholar

Markakis, M. N., Boron, A. K., Van Loock, B., Saini, K., Cirera, S., Verbelen, J. P., et al. (2013). Characterization of a small auxin-up RNA (SAUR)-like gene involved in Arabidopsis thaliana development. PloS One 8, 1–13. doi: 10.1371/journal.pone.0082596

CrossRef Full Text | Google Scholar

Máthé, C., Garda, T., Freytag, C., Hamvas, M. M. (2019). The role of serine-threonine protein phosphatase pp2a in plant oxidative stress signaling—facts and hypotheses. Int. J. Mol. Sci. 20, 3028. doi: 10.3390/ijms20123028

CrossRef Full Text | Google Scholar

McKown, A. D., Guy, R. D., Quamme, L., Klápště, J., La Mantia, J., Constabel, C. P., et al. (2014a). Association genetics, geography and ecophysiology link stomatal patterning in Populus trichocarpa with carbon gain and disease resistance trade-offs. Mol. Ecol. 23, 5771–5790. doi: 10.1111/mec.12969

PubMed Abstract | CrossRef Full Text | Google Scholar

McKown, A. D., Klápště, J., Guy, R. D., Geraldes, A., Porth, I., Hannemann, J., et al. (2014b). Genome-wide association implicates numerous genes underlying ecological trait variation in natural populations of Populus trichocarpa. New Phytol. 203, 535–553. doi: 10.1111/nph.12815

PubMed Abstract | CrossRef Full Text | Google Scholar

Mckown, A. D., Guy, R. D., Klapste, J., Geraldes, A., Friedmann, M., Cronk, Q. C. B., et al. (2014c). Geographical and environmental gradients shape phenotypic trait variation and genetic structure in Populus trichocarpa. New Phytol. 201, 1263–1276. doi: 10.1111/nph.12601

PubMed Abstract | CrossRef Full Text | Google Scholar

McKown, A. D., Klápště, J., Guy, R. D., Corea, O. R. A., Fritsche, S., Ehlting, J., et al. (2019). A role for SPEECHLESS in the integration of leaf stomatal patterning with the growth vs disease trade-off in poplar. New Phytol. 223, 1888–1903. doi: 10.1111/nph.15911

PubMed Abstract | CrossRef Full Text | Google Scholar

McSteen, P., Zhao, Y. (2008). Plant Hormones and Signaling: Common Themes and New Developments. Dev. Cell 14, 467–473. doi: 10.1016/j.devcel.2008.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Muchero, W., Guo, J., DiFazio, S. P., Chen, J., Ranjan, P., Slavov, G. T., et al. (2015). High-resolution genetic mapping of allelic variants associated with cell wall chemistry in Populus. BMC Genomics 16, 24. doi: 10.1186/s12864-015-1215-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Neale, D. B., Kremer, A. (2011). Forest tree genomics: Growing resources and applications. Nat. Rev. Genet. 12, 111–122. doi: 10.1038/nrg2931

PubMed Abstract | CrossRef Full Text | Google Scholar

Neale, D. B., Savolainen, O. (2004). Association genetics of complex traits in conifers. New Phytol. 9, 325–330. doi: 10.1111/j.1469-8137.2010.03593.x

CrossRef Full Text | Google Scholar

Nolan, T. M., Vukašinović, N., Liu, D., Russinova, E., Yin, Y. (2020). Brassinosteroids: Multidimensional Regulators of Plant Growth, Development, and Stress Responses. Plant Cell 32, 295–318. doi: 10.1105/tpc.19.00335

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohashi-Ito, K., Iwamoto, K., Fukuda, H. (2018). LOB DOMAIN-CONTAINING PROTEIN 15 Positively Regulates Expression of VND7, a Master Regulator of Tracheary Elements. Plant Cell Physiol. 59, 989–996. doi: 10.1093/pcp/pcy036

PubMed Abstract | CrossRef Full Text | Google Scholar

Orfila, C., Sørensen, S. O., Harholt, J., Geshi, N., Crombie, H., Truong, H. N., et al. (2005). QUASIMODO1 is expressed in vascular tissue of Arabidopsis thaliana inflorescence stems, and affects homogalacturonan and xylan biosynthesis. Planta 222, 613–622. doi: 10.1007/s00425-005-0008-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Pandey, P., Irulappan, V., Bagavathiannan, M. V., Senthil-Kumar, M. (2017). Impact of combined abiotic and biotic stresses on plant growth and avenues for crop improvement by exploiting physio-morphological traits. Front. Plant Sci. 8, 537. doi: 10.3389/fpls.2017.00537

PubMed Abstract | CrossRef Full Text | Google Scholar

Pendergrass, S. A., Brown-Gentry, K., Dudek, S., Frase, A., Torstenson, E. S., Goodloe, R., et al. (2013). Phenome-Wide Association Study (PheWAS) for Detection of Pleiotropy within the Population Architecture using Genomics and Epidemiology (PAGE) Network. PloS Genet. 9, e1003087. doi: 10.1371/journal.pgen.1003087

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, H. F., O’Reilly, P. F. (2017). Multivariate simulation framework reveals performance of multi-trait GWAS methods. Sci. Rep. 7, 1–12. doi: 10.1038/srep38837

PubMed Abstract | CrossRef Full Text | Google Scholar

Porth, I., Klápště, J., Skyba, O., Friedmann, M. C., Hannemann, J., Ehlting, J., et al. (2013a). Network analysis reveals the relationship among wood properties, gene expression levels and genotypes of natural Populus trichocarpa accessions. New Phytol. 200, 727–742. doi: 10.1111/nph.12419

PubMed Abstract | CrossRef Full Text | Google Scholar

Porth, I., Ranjan, P., Klapšte, J., Guy, R. D., Tuskan, G. A., Ehlting, J., et al. (2013b). Genome-wide association mapping for wood characteristics in Populus identifies an array of candidate single nucleotide polymorphisms. New Phytol. 200, 710–726. doi: 10.1111/nph.12422

PubMed Abstract | CrossRef Full Text | Google Scholar

Porth, I., Klápště, J., McKown, A. D., La Mantia, J., Guy, R. D., Ingvarsson, P. K., et al. (2015). Evolutionary Quantitative Genomics of Populus trichocarpa. PloS One 10, e0142864. doi: 10.1371/journal.pone.0142864

PubMed Abstract | CrossRef Full Text | Google Scholar

Preston, K. A., Cornwell, W. K., DeNoyer, J. L. (2006). Wood density and vessel traits as distinct correlates of ecological strategy in 51 California coast range angiosperms. New Phytol. 170, 807–818. doi: 10.1111/j.1469-8137.2006.01712.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Quan, M., Wang, Q., Phangthavong, S., Yang, X., Song, Y., Du, Q., et al. (2016). Association studies in Populus tomentosa reveal the genetic interactions of PTO-MIR156c and its targets in wood formation. Front. Plant Sci. 7, 1159. doi: 10.3389/fpls.2016.01159

PubMed Abstract | CrossRef Full Text | Google Scholar

Quan, M., Du, Q., Xiao, L., Lu, W., Wang, L., Xie, J., et al. (2019). Genetic architecture underlying the lignin biosynthesis pathway involves noncoding RNAs and transcription factors for growth and wood properties in Populus. Plant Biotechnol. J. 17, 302–315. doi: 10.1111/pbi.12978

PubMed Abstract | CrossRef Full Text | Google Scholar

Rymaszewski, W., Vile, D., Bediee, A., Dauzat, M., Luchaire, N., Kamrowska, D., et al. (2017). Stress-related gene expression reflects morphophysiological responses to water deficit. Plant Physiol. 174, 1913–1930. doi: 10.1104/pp.17.00318

PubMed Abstract | CrossRef Full Text | Google Scholar

Sade, N., Vinocur, B. J., Diber, A., Shatil, A., Ronen, G., Nissan, H., et al. (2009). Improving plant stress tolerance and yield production: Is the tonoplast aquaporin SlTIP2;2 a key to isohydric to anisohydric conversion? New Phytol. 181, 651–661. doi: 10.1111/j.1469-8137.2008.02689.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Schindelin, J., Rueden, C. T., Hiner, M. C., Eliceiri, K. W. (2015). The ImageJ ecosystem: An open platform for biomedical image analysis. Mol. Reprod. Dev. 82, 518–529. doi: 10.1002/mrd.22489

PubMed Abstract | CrossRef Full Text | Google Scholar

Shim, H., Chasman, D., II, Smith, J. D., Mora, S., Ridker, P. M., Nickerson, D. A., et al. (2015). A multivariate genome-wide association analysis of 10 LDL subfractions, and their response to statin treatment, in 1868 Caucasians. PloS One 10, 1–20. doi: 10.1371/journal.pone.0120758

CrossRef Full Text | Google Scholar

Slavov, G. T., Difazio, S. P., Martin, J., Schackwitz, W., Muchero, W., Rodgers-Melnick, E., et al. (2012). Genome resequencing reveals multiscale geographic structure and extensive linkage disequilibrium in the forest tree Populus trichocarpa. New Phytol. 196, 713–725. doi: 10.1111/j.1469-8137.2012.04258.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sperry, J. S. (2000). Hydraulic constraints on plant gas exchange. Agric. For. Meteorol. 104, 13–23. doi: 10.1016/S0168-1923(00)00144-1

CrossRef Full Text | Google Scholar

Sperry, J. S. (2003). Evolution of Water Transport and Xylem Structure. Int. J. Plant Sci. 164, S115–S127. doi: 10.1086/368398

CrossRef Full Text | Google Scholar

Storey, J. D., Tibshirani, R. (2003). Statistical significance for genomewide studies. Proc. Natl. Acad. Sci. U. S. A. 100, 9440–9445. doi: 10.1073/pnas.1530509100

PubMed Abstract | CrossRef Full Text | Google Scholar

Street, N. R., Ingvarsson, P. K. (2011). Association genetics of complex traits in plants. New Phytol. 189, 909–922. doi: 10.1111/j.1469-8137.2010.03593.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Studer, M. H., DeMartini, J. D., Davis, M. F., Sykes, R. W., Davison, B., Keller, M., et al. (2011). Lignin content in natural Populus variants affects sugar release. Proc. Natl. Acad. Sci. 108, 6300–6305. doi: 10.1073/pnas.1009252108

CrossRef Full Text | Google Scholar

Tang, F., Wei, H., Zhao, S., Wang, L., Zheng, H., Lu, M. (2016). Identification of microRNAs involved in regeneration of the secondary vascular system in populus tomentosa carr. Front. Plant Sci. 7, 724. doi: 10.3389/fpls.2016.00724

PubMed Abstract | CrossRef Full Text | Google Scholar

Tschaplinski, T. J., Abraham, P. E., Jawdy, S. S., Gunter, L. E., Martin, M. Z., Engle, N. L., et al. (2019). The nature of the progression of drought stress drives differential metabolomic responses in Populus deltoides. Ann. Bot. 124, 617–626. doi: 10.1093/aob/mcz002. In press.

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuteja, N., Mahajan, S. (2007). Calcium signaling network in plants: An overview. Plant Signal. Behav. 2, 79–85. doi: 10.4161/psb.2.2.4176

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaid, N., Macovei, A., Tuteja, N. (2013). Knights in action: Lectin receptor-like kinases in plant development and stress responses. Mol. Plant 6, 1405–1418. doi: 10.1093/mp/sst033

PubMed Abstract | CrossRef Full Text | Google Scholar

Verger, S., Chabout, S., Gineau, E., Mouille, G. (2016). Cell adhesion in plants is under the control of putative O-fucosyltransferases. J. Cell Sci. 129, e1.2–e1.2. doi: 10.1242/jcs.195537

CrossRef Full Text | Google Scholar

Verma, A., Lucas, A., Verma, S. S., Zhang, Y., Josyula, N., Khan, A., et al. (2018). PheWAS and Beyond: The Landscape of Associations with Medical Diagnoses and Clinical Measures across 38,662 Individuals from Geisinger. Am. J. Hum. Genet. 102, 592–608. doi: 10.1016/j.ajhg.2018.02.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Vining, K. J., Pomraning, K. R., Wilhelm, L. J., Priest, H. D., Pellegrini, M., Mockler, T. C., et al. (2012). Dynamic DNA cytosine methylation in the Populus trichocarpa genome: Tissue-level variation and relationship to gene expression. BMC Genomics 13:27. doi: 10.1186/1471-2164-13-27

PubMed Abstract | CrossRef Full Text | Google Scholar

Visscher, P. M., Wray, N. R., Zhang, Q., Sklar, P., McCarthy, M., II, Brown, M. A., et al. (2017). 10 Years of GWAS Discovery: Biology, Function, and Translation. Am. J. Hum. Genet. 101, 5–22. doi: 10.1016/j.ajhg.2017.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Sun, T., Li, T., Wang, M., Yang, G., He, G. (2016). A CBL-interacting protein kinase TaCIPK2 confers drought tolerance in transgenic tobacco plants through regulating the stomatal movement. PloS One 11, e167962. doi: 10.1371/journal.pone.0167962

CrossRef Full Text | Google Scholar

Wegrzyn, J. L., Eckert, A. J., Choi, M., Lee, J. M., Stanton, B. J., Sykes, R., et al. (2010). Association genetics of traits controlling lignin and cellulose biosynthesis in black cottonwood (Populus trichocarpa, Salicaceae) secondary xylem. New Phytol. 188, 515–532. doi: 10.1111/j.1469-8137.2010.03415.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Weichselbaum, D., Zagrovic, B., Polyansky, A. A. (2017). Fuento: functional enrichment for bioinformatics. Bioinformatics 33, 2604–2606. doi: 10.1093/bioinformatics/btx179

PubMed Abstract | CrossRef Full Text | Google Scholar

Weighill, D., Jones, P., Shah, M., Ranjan, P., Muchero, W., Schmutz, J., et al. (2018). Pleiotropic and Epistatic Network-Based Discovery: Integrated Networks for Target Gene Discovery. Front. Energy Res. 6, 30. doi: 10.3389/fenrg.2018.00030

CrossRef Full Text | Google Scholar

Weighill, D., Jones, P., Bleker, C., Ranjan, P., Shah, M., Zhao, N., et al. (2019). Multi-Phenotype Association Decomposition: Unraveling Complex Gene-Phenotype Relationships. Front. Genet. 10, 417. doi: 10.3389/fgene.2019.00417

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H. C., Bulgakov, V. P., Jinn, T. L. (2018). Pectin methylesterases: Cell wall remodeling proteins are required for plant response to heat stress. Front. Plant Sci. 871, 1612. doi: 10.3389/fpls.2018.01612

CrossRef Full Text | Google Scholar

Xi, W., Song, D., Sun, J., Shen, J., Li, L. (2017). Formation of wood secondary cell wall may involve two type cellulose synthase complexes in Populus. Plant Mol. Biol. 93, 419–429. doi: 10.1007/s11103-016-0570-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, M., Zhang, J., Tschaplinski, T. J., Tuskan, G. A., Chen, J.-G., Muchero, W. (2018). Regulation of Lignin Biosynthesis and Its Role in Growth-Defense Tradeoffs. Front. Plant Sci. 9, 1427. doi: 10.3389/fpls.2018.01427

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Labbé, J., Muchero, W., Yang, X., Jawdy, S. S., Kennedy, M., et al. (2016). Genome-wide analysis of lectin receptor-like kinases in Populus trichocarpa. BMC Genomics 17, 1–16. doi: 10.1186/s12864-016-3026-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Yee, D., Goring, D. R. (2009). The diversity of plant U-box E3 ubiquitin ligases: From upstream activators to downstream target substrates. J. Exp. Bot. 60, 1109–1121. doi: 10.1093/jxb/ern369

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Holland, J. B., McMullen, M. D., Buckler, E. S. (2008). Genetic design and statistical power of nested association mapping in maize. Genetics. 178, 539–551. doi: 10.1534/genetics.107.074245

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, H., Xu, L., Singh, A., Wang, H., Du, L., Poovaiah, B. W. (2015). Involvement of calmodulin and calmodulin-like proteins in plant responses to abiotic stresses. Front. Plant Sci. 6, 600. doi: 10.3389/fpls.2015.00600

PubMed Abstract | CrossRef Full Text | Google Scholar

Zentella, R., Zhang, Z.-L., Park, M., Thomas, S. G., Endo, A., Murase, K., et al. (2007). Global Analysis of DELLA Direct Targets in Early Gibberellin Signaling in Arabidopsis. Plant Cell 19, 3037–3057. doi: 10.1105/tpc.107.054999

PubMed Abstract | CrossRef Full Text | Google Scholar

Zentella, R., Sui, N., Barnhill, B., Hsieh, W. P., Hu, J., Shabanowitz, J., et al. (2017). The Arabidopsis O-fucosyltransferase SPINDLY activates nuclear growth repressor della. Nat. Chem. Biol. 13, 479–485. doi: 10.1038/nchembio.2320

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhan, X., Hu, Y., Li, B., Abecasis, G. R., Liu, D. J. (2016). RVTESTS: An efficient and comprehensive tool for rare variant association analysis using sequence data. Bioinformatics 32, 1423–1426. doi: 10.1093/bioinformatics/btw079

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Yang, Y., Zheng, K., Xie, M., Feng, K., Jawdy, S. S., et al. (2018). Genome-wide association studies and expression-based quantitative trait loci analyses reveal roles of HCT2 in caffeoylquinic acid biosynthesis and its regulation by defense-responsive transcription factors in Populus. New Phytol. 220, 502–516. doi: 10.1111/nph.15297

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Stephens, M. (2012). Genome-wide efficient mixed-model analysis for association studies. Nat. Genet. 44, 821–826. doi: 10.1038/ng.2310

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Stephens, M. (2014). Efficient multivariate linear mixed model algorithms for genome- wide association studies. Nat. Genet. 11, 407–411. doi: 10.1038/nmeth.2848

CrossRef Full Text | Google Scholar