Results
We first performed a characterisation of the two heteroblastic leaf forms produced byP. tremula . Pre-formed leaves were orbicular and those produced from the extending shoot apical meristem (terminal) were cordate (Figure 1A, B), which is in agreement with previous observations in poplar (Y. Liu et al., 2015; Russin & Evert, 1984). Pre-formed leaves were notably thicker (~2x) than terminal leaves, with contrasting spatial arrangement of cell layers between the two leaf types (Figure 1C, D). Pre-formed leaves had a thicker adaxial epidermal layer followed by two rows of densely packed palisade mesophyll; spongy mesophyll cells were separated by many air spaces and a thick abaxial epidermis. The spongy mesophyll cell layer of terminal leaves had fewer air spaces than in pre-formed leaves (Figure 1C, D).

Genome wide association mapping identifies complex genetic architecture

We characterised variation in pre-formed leaf shape within the Swedish Aspen (SwAsp) collection of P. tremula genotypes, which was sampled from local populations in Sweden and established in two replicated common garden experiments located in the north (Sävar, Umeå, 62°N) and south (Ekebo, 56°N) of Sweden (Luquez et al., 2008). We measured three representative leaf physiognomy traits: leaf area, circularity and indent depth (margin/boundary serration) in both common gardens in two years. We observed extensive natural variation in pre-formed leaf shape (Figure 2A), with measured leaf shape traits having high broad sense heritability (H 2) and low sub-population differentiation (Q ST) (Table 1). Genetic correlations indicated that a substantial proportion of the heritable variation for each trait is controlled by genetic factors unique to that trait (Figure 2B; Table 1). We additionally tested whether any of the traits correlated significantly to a range of environmental and climatic factors (Table 1). There was no evidence of sub-population differentiation for any of the traits, as indicated by low Q ST (Table 1) and genotypes did not cluster by sub-population of origin on the basis of the three phenotypic traits (Figure 2D). Similarly, there were no large-effect significant correlations between the traits and environmental and climatic factors (Table 1). Leaf circularity and indent depth had higherH 2 and lower Genotype by Environment (GxE) interaction than did leaf area (indicated by comparisons of ANOVAF values; Table 1; Figure S1).
We used the trait values to perform association mapping using genome-wide single nucleotide polymorphism (SNP) markers identified from population-wide resequencing data (Wang et al., 2018). Since there were no statistically significant associations after multiple testing correction for any trait in any of the repeated datasets, we examined whether there was overlap among the top 1,000 SNP associations (rank ordered by P -value), indicative of consistency in ranking among the repeated measures. In general, there were few associations in common among the repeated measures, with the majority of the top-ranked associations being unique to each dataset (Figure 2C). In agreement with the lower GxE and higher H 2 for leaf circularity and indent depth, these traits displayed greater overlap among the associations than for leaf area. We also calculated Best Linear Unbiased Predictions (BLUPs) from the repeated measures of each trait, but again found no statistically significant associations after multiple testing correction for the BLUPs (Figure S2, Table S1). There was no evidence of substantial inflation due to population stratification (Yang et al., 2011), as indicated by Genomic Control (GC) values (λGC 1.05 for area, 1.04 for circularity, 1.04 for indent depth) while SNP-based estimates of percentage variance explained (PVE) were relatively high (0.40 ± s.d. 0.31 for area, 0.63 ± s.d. 0.31 for indent depth, 0.80 ± s.d. 0.28 for circularity), further supporting that higher plasticity of leaf area and that trait variance was under genetic control.
To select candidate genes for these traits, we again rank ordered the SNP associations by P -value and selected the top 1,000 genes for each of the traits (Table S2). The majority of associations were unique to each trait (Figure S3), in line with genetic correlations suggesting largely independent genetic control (Figure 2B). We examined the genomic context of SNPs within the top 1,000 gene sets (Figure 3), observing that the highest density of SNPs occurred in regulatory regions (UTRs and flanking regions, which contain the promoter).

Differential gene expression between genotypes with contrasting leaf shape

We used an existing resource assaying gene expression in flushing leaf buds from the SwAsp collection (Mähler et al., 2017) to examine correlations between gene expression and the leaf physiognomy traits. Although none of these correlations were significant (Figure S4A), we did find significantly Differentially Expressed Genes (DEGs) between sets of genotypes at the population distribution extremes for the three phenotypic measures (Figure 4) including 182 DEGs for area, 203 for circularity and 223 for indent depth (Table S3). As the two shape traits had higher H 2, we concentrated more specifically on these in subsequent analyses.

Extensive remodulation of the transcriptome during leaf development

To characterise the developmental role of our candidate genes, we profiled gene expression using RNA-Seq in a developmental series of terminal leaves from a single genotype. Developmental time accounted for the largest proportion of variance in the corresponding gene expression data with the apical region sample being distinctly separated from later developmental stages and with more extensive separation of later stages (Figure 5A). We used these data to perform differential expression analysis (Figure 5B) and to confirm that we observed the expected expression profiles for Gene Ontology (GO) categories (Figure 5C) and homologs of known leaf development regulators (Figure 5D). To be able to examine the relationship to network topology of our GWAS and DEG candidate genes we performed unsupervised network analysis to obtain an unbiased overview of major expression profiles and underlying processes active during leaf development and to identify the most central genes involved in these processes. We calculated a gene co-expression network (Table 2) by aggregating networks from multiple inference algorithms (Marbach et al., 2012). Our analysis involved graph partitioning to define modules (clusters) of genes and subsequent calculation of node centrality statistics (Table 3, Table S4). We identified statistical enrichment of GO and PFAM (Protein family) categories within graph clusters to annotate common processes represented by cluster members (Figure S5 Table S5). We overlaid stage-wise DE results (Table S6) onto the network nodes to visualise network regions or clusters active at each leaf development stage transition and to examine whether most significantly DE genes (DEGs) exhibited an increase or decrease in expression (Figure 6).
To aid community utilisation, we have made the gene expression data and the co-expression network available at http://aspleaf.plantgenie.org and within PopGenIE (Sundell et al., 2015) as the dataset ‘AspLeaf’.

Developmental characterisation of population candidate genes

The majority of the GWAS candidate genes were expressed within the sampled leaf ages but were not included in the developmental network, indicating low levels of connectivity (436 GWAS leaf area genes, 451 circularity genes and 424 indent depth genes were included in the network). None of the genes that were present in the network were highly ranked in terms of network connectivity metrics and none of the network clusters were enriched for these GWAS genes. Furthermore, we performed Gene Set Enrichment Analysis (GSEA) of the GWAS genes using network PageRank score as a measure of network centrality to rank-order genes within the developmental network, finding no significant enrichment (Daub et al., 2013). However, the trend was for GWAS genes being present among low connectivity genes (Figure S6). As such, genes identified by GWAS were clearly not central within the developmental network.
We performed GO enrichment analysis of the GWAS candidate genes, however none had significant enrichment, suggesting that they included genes spanning a diverse range of biological processes and that no particular process was associated with the genetic control of trait variation among these genes. As an alternative to analysing the discrete set of GWAS candidate genes, we used the SNP with the strongest association (lowestP -value) within each gene to rank order all genes for each trait. We then performed GSEA using gene sets from GO terms and expression network clusters. A number of GO terms were enriched for the three traits, including “protein phosphorylation” (GO:0006468) for circularity, “sulfur compound metabolic process” (GO:0006790) for indent depth, and “carbohydrate metabolic process” (GO:0005975) for area. Among the network clusters, five were significant at a 5% false discovery rate. Clusters [1:4], [1:5], and [1:10] were enriched in the leaf area GWAS results and clusters [1:1], [1:4], and [1:5] were enriched for circularity. No cluster was significant for indent depth.
Many of the DEGs between leaf shape extremes were actively regulated during the process of leaf development (Figure 7; Figure S7). There was no significant enrichment of GO terms among these sets of DEGs, however GSEA showed that indent depth GWAS genes were significantly enriched in indent depth DEGs (q =0.0009), that circularity GWAS genes were significantly enriched in circularity DEGs (q =0.009), but that there was no enrichment for area. Examination of the correlations between DEGs and the corresponding phenotypic traits showed that, although no single correlation was significant, correlation values for the DEGs were significantly larger than for non-DEGs (Figure S4B). Two of these DEGs were in common for all three traits (Figure S8) but had relatively low and consistent expression within the leaf development datasets: Potra196739g30199 (ATP synthase subunit C) and Potra003791g32371 (mRNA cap guanine-N7 methyltransferase). Within the terminal development network, 71 area, 83 circularity and 85 indent depth DEGs were present. We performed GSEA of the sets of DEGs within the terminal development network (Figure S6), which revealed significant under-representation for genes with low connectivity for all three traits (P =0.004, 0.011 and 0.001, respectively).
The population gene expression data has previously been utilised to perform expression Quantitative Trait Locus (eQTL) mapping and co-expression network analyses (Mähler et al., 2017). A relatively large proportion of the DEGs between leaf shape extremes were also eGenes (i.e. genes for which an eQTL was present): 80 for area, 78 for circularity, and 83 for indent depth. Out of these genes, 22 for circularity and 35 for indent depth were also part of the subset of genes that contributed to the GWAS GSEA signal. We examined the distribution of network centrality (degree) of DEGs and non-DEGs and GWAS genes within the population co-expression network presented in Mähler et al. (2017) and the developmental network (Figure S9) and performed GSEA to test for a relationship to network connectivity within the population co-expression network (Figure S10). DEGs had significantly lower centrality than non-DEGs, a pattern that we confirmed was specific only to DEGs and not to random gene sets (Figure S9). The GSEA revealed that GWAS genes were significantly enriched among low connectivity genes while DEGs were not, although the trend for DEGs was similar (Figure S10). Furthermore, DEGs displayed signatures of relaxed selective constraint (i.e. lower negative selection), as indicated by Tajima’s D values (Figure S11). This pattern was observed for DEGs but not for the top ranked 1,000 GWAS gene sets.

Discussion

Despite the vast variation in leaf shape within and between species, we still have a rather limited understanding of the gene regulatory network underlying the process of leaf development and of the genetic determinants of leaf shape variation among individuals. There is extensive variation in leaf form within an individual (Figure 1) as well as shape variation of pre-formed leaves among individuals (Figure 2A) in our current study system. We used leaf area, circularity and indent depth as representative physiognomy traits, with leaf shape being under tight genetic control (Table 1). There was no apparent link between leaf shape or size and environment, longitude or latitude (Table 1) andQ ST values were low. As such, there was no identifiable adaptive role or clear signature of directional selection for leaf shape within the SwAsp collection. In line with this, and in stark contrast to the date of autumn bud set, which is highly adaptive (Wang et al., 2018), we did not identify any significant SNP associations for the three traits. Given the highH 2 for leaf circularity and indent depth, this indicates complex genetic architecture with a large number of small effect size polymorphisms contributing to the control of leaf shape variation. As the population size of the SwAsp collection is certainly underpowered to detect such small effect associations, we reasoned that calculated P -values would still be informative for rank ordering the importance of SNP associations. We therefore used this rank ordering to identify the top 1,000 associated genes (i.e. GWAS candidate genes) and examined these further. After normalising for feature length, the highest density of SNPs within the candidate genes was observed in UTR and regulatory (up- and down-stream) regions (Figure 3), suggesting that those SNPs potentially act by influencing gene expression.
We note that these findings are highly concordant with the omnigenic model (Boyle, Li, & Pritchard, 2017; X. Liu, Li, & Pritchard, 2019), which posits that a majority of genes contributing to trait variance are not directly biologically connected with the trait of interest (X. Liu et al., 2019). The model states that there is a limited set of genes that directly affect a trait (referred to as core genes), while all other genes with expression variation in the relevant tissue are peripheral to the trait. These peripheral genes can havetrans -acting effects that propagate through the highly connected regulatory network that cause individually small changes in the expression of core genes. While each trans effect is small in isolation, the combined effect of these peripheral trans -acting effects will explain the majority of the variation in core genes and, therefore, in trait variation. It has not been resolved whether core genes defined on the basis of topology within a co-expression network represent core genes as defined by the omnigenic model. It also remains to be explored whether the ability to predict phenotype from gene expression profiles, or co-expression modules, is an effective means of identifying omnigenic core genes.
We further explored network topology of our identified GWAS candidate genes. A common characteristic of core genes (as defined by the omnigenic model) and of network hubs is that they are under strong selective constraint as natural selection acts strongly against large effect variants (Gouy, Daub, & Excoffier, 2017; X. Liu et al., 2019; Simons, Bullaughey, Hudson, & Sella, 2018). To test whether this held up in our system, we focused on a reference genotype to establish a developmental timeline for terminal leaves. As predicted by the model, our GWAS candidate genes were not hubs within the co-expression network constructed from the leaf developmental series.
Our observation that the highest density of top-ranked SNPs was within regulatory regions prompted us to make use of an existing resource of gene expression data assaying flushing buds from the same population (Mähler et al., 2017). We have previously performed eQTL mapping using these data, the results of which also highlighted the importance of SNPs located within regulatory regions in controlling natural variation in transcript abundance (Mähler et al., 2017). We did not identify any genes with significant correlation of gene expression variation among genotypes to variation in our target leaf traits and correlations were of low magnitude (Figure S4). This is, however, expected given the apparently highly polygenic nature of these traits.
In addition to the GWAS candidate genes, we also derived a set of candidate genes based on the population expression data. To this end, we identified a set of genes with the maximum signal strength between expression and phenotype by performing DE analyses between sets of genotypes with the most extreme phenotype values (Figure 4). There was significantly higher correlation for DEGs than non-DEGs and three-way overlap between GWAS, DEGs and eQTL presence provided circumstantial evidence that these genes may influence leaf shape via gene expression variation and that this effect is under genetic control. The sets of DEGs were of lower centrality within the population wide co-expression network (Figure S10) and displayed evidence of weaker negative selection (Figure S11). However, for area and indent depth there was an under-representation of low connectivity genes within the development co-expression network. As such, these candidate genes were under-represented for genes for which changes to expression would be expected to have the lowest impact to the network in general and, by extension, on phenotype. The majority of these genes have no known function in leaf development, but many were actively regulated during leaf development (Figure 7, Figure S7). These results are highly concordant with the omnigenic model, indicating that the DEG analysis potentially identified a set of core trait genes (as defined by the omnigenic model) that have higher effect size on phenotype and higher correlation of expression to phenotype alongside a larger set of genes that are peripheral to the phenotype and that are distributed throughout the developmental co-expression network. Within the GWAS candidate genes there was no such under-representation, suggesting that the DEGs were enriched for genes of higher importance within the development co-expression network and, therefore, of larger impact on phenotype, although this signal is weak. The lack of GO enrichment within the GWAS or DEG candidate gene sets is also congruent with the omnigenic model as the majority of genes contributing to trait heritability are expected to be peripheral and not associated to particular biological processes or directly connected to the focal trait. Taken together our results suggest that leaf shape is an omnigenic trait, with variation resulting from SNPs within regulatory regions that act by causing variation in gene expression. The genes affected are enriched within the periphery of the population co-expression network and are associated with signatures of relaxed selective constraint. There was a lack of evidence for directional selection acting on leaf shape, as indicated by low sub-population differentiation and a lack of any link to environmental variables. This potentially indicates that variation in leaf shape is adaptively neutral within the morphospace represented among the sampled genotypes, although there are other plausible explanations congruent with the results. However, neutrality does appear to be an harmonious interpretation as allele frequencies of SNPs amoung the GWAS candidate genes reflected those of all SNPs globally and of the 1000 lowest ranked SNPs. There was, therefore, no clear indication of balancing or stabilising selection or of extensive purifying selection. Of note, the pattern was not the same for area, for which SNPs within GWAS candidate genes showed a skew towards lower allele frequencies. This suggests that there may be contrasting selection pressures acting on leaf size and shape, with leaf size being subjected to stronger purifying selection, although the signal was weak, and this was not reflected in the distribution of Tajima’s D (Figure S11).
Similar approaches have been used to identify and prioritise candidate genes in maize, with small numbers of those candidates confirmed to influence leaf characteristics in transgenic lines (Baute et al., 2015, 2016; Schaefer et al., 2018). A distinct difference in maize is the clear separation of cell proliferation and subsequent expansion into defined zones that can easily be sampled separately. This sampling approach increases signal strength for associating gene expression to division and expansion processes, which were mixed within single samples within our current study. In maize there are reported loci of higher effect size influencing leaf development and phenotype, possibly as those traits are correlated to general yield characteristics (Baute et al., 2016) and have therefore been targets of artificial selection. A similar study in apple reported a small number of significant associations for leaf shape traits but concluded that those were false-positives (Crouch et al., 2018) and a study in sweet potato used a similar DEG approach to identify genes associated with variation in leaf form (Gupta et al., 2019).
The analysis of leaf shape has many parallels with that of face shape in humans, for which initial GWAS of basic parameters (such as length:width) yielded few to no significant associations. However, analyses that defined more specific sub-features have yielded substantially more informative associations (Crouch et al., 2018). In aspen, there is similarly considerable scope to improve phenotype GWAS results through improved software tools to deconvolute more specific features of leaf shape, such as the angle of veins or of specific regions of the leaf, such as the base and tip. Each of these features is likely under specific molecular control that exhibits genetic variation among individuals and that are not captured by the larger scale features, such as overall leaf circularity, employed in our current study. The fact that such sub-features of leaf shape are under local, spatially defined control and that growth distribution across the leaf lamina controls final leaf shape and venation patterning (Runions, Tsiantis, & Prusinkiewicz, 2017) may also explain the low correlation between leaf shape traits and gene expression in the whole-bud samples used to generate the gene expression values for the eQTL mapping results that we considered here. Furthermore, the use of a single developmental snapshot for such analysis is limiting as causal variation in expression could occur at any point during the developmental program.
We integrated developmental gene expression profiling, GWAS, population-wide gene expression data and population genetics to define the genetic architecture of leaf shape variation in aspen. We show that leaf shape variation is a highly complex trait likely determined by small-effect variations in gene expression caused by numerous small-effect size SNPs in regulatory regions. Genes with evidence of association to variation in leaf shape were peripheral within a population-wide gene co-expression network, were not hubs within the leaf developmental network and displayed signatures of relaxed selection. We therefore suggest that leaf shape is an omnigenic trait. Combined with low sub-population differentiation and a lack of correlation to climatic or environmental variables, we suggest that variation in leaf shape within the morphospace represented within the SwAsp collection may be neutral. We believe that this integrated approach is a pragmatic strategy for identifying candidate genes and disentangling the genetic basis of complex traits, such as leaf shape variation in aspen, which have no apparent evidence of being adaptive and for which no significant associations are identified using association mapping.

Acknowledgements

NRS, BT and KMR were supported by the Trees and Crops for the Future (TC4F) project. NM, BS and CM were supported by the Knut and Alice Wallenberg foundation and the VINNOVA UPSC Centre for Forest Biotechnology. This work was supported by the Trygger Foundation (#CTS 12:471) and Gunnar och Ruth Björkmans fond för Norrländsk Botanisk Forskning. The authors acknowledge support from the Umeå Plant Science Centre Bioinformatics platform (UPSCb), the Science for Life Laboratory (SciLifeLab)and the National Genomics Infrastructure (NGI). Computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX). We thank Skogforsk for hosting the SwAsp collection common gardens.

References

Aida, M., Ishida, T., Fukaki, H., Fujisawa, H., & Tasaka, M. (1997). Genes involved in organ separation in Arabidopsis: an analysis of the cup-shaped cotyledon mutant. The Plant Cell , 9 (6), 841–857. doi: 10.1105/tpc.9.6.841
Anders, S., Pyl, P. T., & Huber, W. (2014). HTSeq - A Python framework to work with high-throughput sequencing data. Bioinformatics (Oxford, England) , 31 (2), 166–169. doi: 10.1101/002824
Andres, R. J., Coneva, V., Frank, M. H., Tuttle, J. R., Samayoa, L. F., Han, S.-W., … Kuraparthy, V. (2017). Modifications to a LATE MERISTEM IDENTITY1 gene are responsible for the major leaf shapes of Upland cotton ( Gossypium hirsutum L.). Proceedings of the National Academy of Sciences , 114 (1), E57–E66. doi: 10.1073/pnas.1613593114
Andrews, S. (n.d.). FastQC. Retrieved from http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ website: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
Ashburner, M., Ball, C., Blake, J., Botstein, D., Butler, H., Cherry, J., … Sherlock, G. (2000). Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nature Genetics , 25 (1), 25–29. doi: 10.1038/75556
Barnes, B. V. (1969). Natural variation and delineation of clones of Populus tremuloides and P. grandidentata in northern lower Michigan.Silvae Genetica , 18 , 130–142.
Barnes, B. V. (1975). Phenotypic variation of Trembling Aspen in western North America. Forest Science , 21 (3), 319–328.
Baute, J., Herman, D., Coppens, F., De Block, J., Slabbinck, B., Dell’Acqua, M., … Inzé, D. (2015). Correlation analysis of the transcriptome of growing leaves with mature leaf parameters in a maize RIL population. Genome Biology , 16 (1), 168. doi: 10.1186/s13059-015-0735-9
Baute, J., Herman, D., Coppens, F., De Block, J., Slabbinck, B., Dell’Acqua, M., … Inzé, D. (2016). Combined Large-Scale Phenotyping and Transcriptomics in Maize Reveals a Robust Growth Regulatory Network. Plant Physiology , 170 (3), 1848–1867. doi: 10.1104/pp.15.01883
Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics ,30 (15), 2114–2120. doi: 10.1093/bioinformatics/btu170
Boyle, E. A., Li, Y. I., & Pritchard, J. K. (2017). An Expanded View of Complex Traits: From Polygenic to Omnigenic. Cell , 169 (7), 1177–1186. doi: 10.1016/J.CELL.2017.05.038
Bylesjö, M., Segura, V., Soolanayakanahally, R. Y., Rae, A. M., Trygg, J., Gustafsson, P., … Street, N. R. (2008). LAMINA: a tool for rapid quantification of leaf size and shape parameters. BMC Plant Biology , 8 (1), 82. doi: 10.1186/1471-2229-8-82
Cheng, X., Peng, J., Ma, J., Tang, Y., Chen, R., Mysore, K. S., & Wen, J. (2012). NO APICAL MERISTEM (MtNAM) regulates floral organ identity and lateral organ separation in Medicago truncatula. New Phytologist , 195 (1), 71–84. doi: 10.1111/j.1469-8137.2012.04147.x
Chitwood, D. H. H., & Sinha, N. R. R. (2016). Evolutionary and Environmental Forces Sculpting Leaf Development. Current Biology ,26 (7), R297–R306. doi: 10.1016/j.cub.2016.02.033
Chitwood, D. H., Kumar, R., Headland, L. R., Ranjan, A., Covington, M. F., Ichihashi, Y., … Sinha, N. R. (2013). A Quantitative Genetic Basis for Leaf Morphology in a Set of Precisely Defined Tomato Introgression Lines. The Plant Cell , 25 (7), 2465–2481. doi: 10.1105/tpc.113.112391
Chitwood, D. H., Maloof, J. N., & Sinha, N. R. (2013). Dynamic transcriptomic profiles between tomato and a wild relative reflect distinct developmental architectures. Plant Physiology ,162 (2), 537–552. doi: 10.1104/pp.112.213546
Chitwood, D. H., Ranjan, A., Kumar, R., Ichihashi, Y., Zumstein, K., Headland, L. R., … Sinha, N. R. (2014). Resolving Distinct Genetic Regulators of Tomato Leaf Shape within a Heteroblastic and Ontogenetic Context. The Plant Cell , 26 (9), 3616–3629. doi: 10.1105/tpc.114.130112
Chitwood, D. H., Ranjan, A., Martinez, C. C., Headland, L. R., Thiem, T., Kumar, R., … Sinha, N. R. (2014). A modern ampelography: a genetic basis for leaf shape and venation patterning in grape.Plant Physiology , 164 (1), 259–272. doi: 10.1104/pp.113.229708
Clark, S. E. (1997). Organ Formation at the Vegetative Shoot Meristem.The Plant Cell , 9 (7), 1067–1076. doi: 10.1105/tpc.9.7.1067
Clauw, P., Coppens, F., Korte, A., Herman, D., Slabbinck, B., Dhondt, S., … Inzé, D. (2016). Leaf Growth Response to Mild Drought: Natural Variation in Arabidopsis Sheds Light on Trait Architecture.The Plant Cell , 28 (10), 2417–2434. doi: 10.1105/tpc.16.00483
Cleland, R. E. (2001). Unlocking the mysteries of leaf primordia formation. Proceedings of the National Academy of Sciences of the United States of America , 98 (20), 10981–10982. doi: 10.1073/pnas.211443498
Coscia, M., & Neffke, F. (2017). Network Backboning with Noisy Data.ArXiv , arXiv:1701 .
Cox, S. E. (2005). Elevational gradient of neoformed shoot growth inPopulus tremuloides . Canadian Journal of Botany ,83 (10), 1340–1344. doi: 10.1139/b05-103
Critchfield, W. B. (1960). Leaf Dimorphism in Populus Trichocarpa.American Journal of Botany , 47 (8), 699. doi: 10.2307/2439521
Crouch, D. J. M., Winney, B., Koppen, W. P., Christmas, W. J., Hutnik, K., Day, T., … Bodmer, W. F. (2018). Genetics of the human face: Identification of large-effect single gene variants. Proceedings of the National Academy of Sciences of the United States of America ,115 (4), E676–E685. doi: 10.1073/pnas.1708207114
Curtis, J. D., & Lersten, N. R. (1978). Heterophylly in Populus grandidentata (Salicaceae) with Emphasis on Resin Glands and Extrafloral Nectaries. American Journal of Botany , 65 (9), 1003. doi: 10.2307/2442687
Czesnick, H., & Lenhard, M. (2015). Size control in plants–lessons from leaves and flowers. Cold Spring Harbor Perspectives in Biology , 7 (8), a019190. doi: 10.1101/cshperspect.a019190
Daub, J. T., Hofer, T., Cutivet, E., Dupanloup, I., Quintana-Murci, L., Robinson-Rechavi, M., & Excoffier, L. (2013). Evidence for polygenic adaptation to pathogens in the human genome. Molecular Biology and Evolution , 30 (7), 1544–1558. doi: 10.1093/molbev/mst080
Delhomme, N., Mähler, N., Schiffthaler, B., Sundell, D., Mannapperuma, C., Hvidsten, T. R., … Street, N. R. (2015). Guidelines for RNA-Seq data analysis. Retrieved March 13, 2016, from Epigenesis website: http://www.epigenesys.eu/en/protocols/bio-informatics/1283-guidelines-for-rna-seq-data-analysis
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., … Gingeras, T. R. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics (Oxford, England) , 29 (1), 15–21. doi: 10.1093/bioinformatics/bts635
Drost, D. R., Puranik, S., Novaes, E., Novaes, C. R. D. B., Dervinis, C., Gailing, O., & Kirst, M. (2015). Genetical genomics of Populus leaf shape variation. BMC Plant Biology , 15 , 166. doi: 10.1186/s12870-015-0557-7
Du, J., Mansfield, S. D., & Groover, A. T. (2009). The Populus Homeobox Gene ARBORKNOX2 Regulates Cell Differentiation During Secondary Growth.The Plant Journal : For Cell and Molecular Biology ,60 (6), 1000–1014. doi: 10.1111/j.1365-313X.2009.04017.x
Eddelbuettel, D., & Balamuta, J. J. (2017). Extending R with C++: A Brief Introduction to Rcpp. PeerJ Preprints , 5 , e3188v1. doi: 10.7287/peerj.preprints.3188v1
Erickson, R. O., & Michelini, F. J. (1957). The Plastochron Index.American Journal of Botany , 44 (4), 297–305.
Finn, R., Mistry, J., Tate, J., Coggill, P., Heger, A., Pollington, J., … Bateman, A. (2010). The Pfam protein families database.Nucleic Acids Research , 38 (suppl 1), D211–D222. doi: doi: 10.1093/nar/gkp985
Fulop, D., Ranjan, A., Ofner, I., Covington, M. F., Chitwood, D. H., West, D., … Sinha, N. R. (2016). A New Advanced Backcross Tomato Population Enables High Resolution Leaf QTL Mapping and Gene Identification. G3 (Bethesda, Md.) , 6 (10), 3169–3184. doi: 10.1534/g3.116.030536
Gouy, A., Daub, J. T., & Excoffier, L. (2017). Detecting gene subnetworks under selection in biological pathways. Nucleic Acids Research , 45 (16). doi: 10.1093/nar/gkx626
Grossmann, S., Bauer, S., Robinson, P. N., & Vingron, M. (2007). Improved detection of overrepresentation of Gene-Ontology annotations with parent child analysis. Bioinformatics , 23 (22), 3024–3031. doi: 10.1093/bioinformatics/btm440
Gupta, S., Rosenthal, D. M., Stinchcombe, J. R., & Baucom, R. S. (2019). The remarkable morphological diversity of leaf shape in sweet potato ( Ipomoea batatas ): the influence of genetics, environment, and G×E. New Phytologist , nph.16286. doi: 10.1111/nph.16286
Ichihashi, Y., Aguilar-Martínez, J. A., Farhi, M., Chitwood, D. H., Kumar, R., Millon, L. V, … Sinha, N. R. (2014). Evolutionary developmental transcriptomics reveals a gene network module regulating interspecific diversity in plant leaf shape. Proceedings of the National Academy of Sciences of the United States of America ,111 (25), E2616-21. doi: 10.1073/pnas.1402835111
Kopylova, E., Noé, L., & Touzet, H. (2012). SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data.Bioinformatics (Oxford, England) , 28 (24), 3211–3217. doi: 10.1093/bioinformatics/bts611
Larson, P. R., & Isebrands, J. G. (1971). The Plastochron Index as Applied to Developmental Studies of Cottonwood. Canadian Journal of Forest Research , 1 (1), 1–11.
Lin, S. M., Du, P., Huber, W., & Kibbe, W. A. (2008). Model-based variance-stabilizing transformation for Illumina microarray data.Nucleic Acids Research , 36 (2), e11–e11. doi: 10.1093/nar/gkm1075
Lindtke, D., González-Martínez, S. C., Macaya-Sanz, D., & Lexer, C. (2013). Admixture mapping of quantitative traits in Populus hybrid zones: power and limitations. Heredity , 111 (6), 474–485. doi: 10.1038/hdy.2013.69
Liu, X., Li, Y. I., & Pritchard, J. K. (2019). Trans Effects on Gene Expression Can Drive Omnigenic Inheritance. Cell , 177 (4), 1022-1034.e6. doi: 10.1016/j.cell.2019.04.014
Liu, Y., Li, X., Chen, G., Li, M., Liu, M., & Liu, D. (2015). Epidermal Micromorphology and Mesophyll Structure of Populus euphratica Heteromorphic Leaves at Different Development Stages. PloS One ,10 (9), e0137701. doi: 10.1371/journal.pone.0137701
Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology , 15 (12), 550. doi: 10.1186/s13059-014-0550-8
Luquez, V., Hall, D., Albrectsen, B. R., Karlsson, J., Ingvarsson, P., & Jansson, S. (2008). Natural phenological variation in aspen (Populus tremula): the SwAsp collection. Tree Genetics & Genomes ,4 (2), 279–292. doi: 10.1007/s11295-007-0108-y
Mähler, N., Wang, J., Terebieniec, B. K. B. K., Ingvarsson, P. K. P. K., Street, N. R. N. R., & Hvidsten, T. R. (2017). Gene co-expression network connectivity is an important determinant of selective constraint. PLOS Genetics , 13 (4), e1006402. doi: 10.1371/journal.pgen.1006402
Marbach, D., Costello, J. C., Küffner, R., Vega, N. N. M. N., Prill, R. J., Camacho, D. M., … Stolovitzky, G. (2012). Wisdom of crowds for robust gene network inference. Nature Methods , 9 (8), 796–804. doi: 10.1038/nmeth.2016
Meicenheimer, R. D. (2014). The plastochron index: Still useful after nearly six decades. American Journal of Botany , 101 (11), 1821–1835. doi: 10.3732/ajb.1400305
Moon, J., & Hake, S. (2011). How a leaf gets its shape. Current Opinion in Plant Biology , 14 (1), 24–30. doi: 10.1016/J.PBI.2010.08.012
Nystedt, B. B., Street, N. R., Wetterbom, A., Zuccolo, A., Lin, Y.-C., Scofield, D. G., … Jansson, S. (2013). The Norway spruce genome sequence and conifer genome evolution. Nature , 497 (7451), 579–584. doi: 10.1038/nature12211
Peppe, D. J., Royer, D. L., Cariglino, B., Oliver, S. Y., Newman, S., Leight, E., … Wright, I. J. (2011). Sensitivity of leaf size and shape to climate: global patterns and paleoclimatic applications.New Phytologist , 190 (3), 724–739. doi: 10.1111/j.1469-8137.2010.03615.x
R Core Team. (2018). R: A language and environment for statistical computing (R. F. for S. Computing, Ed.). Retrieved from http://www.r-project.org
Rae, A. M., Ferris, R., Tallis, M. J., & Taylor, G. (2006). Elucidating genomic regions determining enhanced leaf growth and delayed senescence in elevated CO2. Plant, Cell and Environment , 29 (9), 1730–1741. doi: 10.1111/j.1365-3040.2006.01545.x
Ranjan, A., Budke, J. M., Rowland, S. D., Chitwood, D. H., Kumar, R., Carriedo, L., … Sinha, N. R. (2016). eQTL Regulating Transcript Levels Associated with Diverse Biological Processes in Tomato.Plant Physiology , 172 (1), 328–340. doi: 10.1104/pp.16.00289
Ratke, C., Terebieniec, B. K., Winestrand, S., Derba-Maceluch, M., Grahn, T., Schiffthaler, B., … Mellerowicz, E. J. (2018). Downregulating aspen xylan biosynthetic GT43 genes in developing wood stimulates growth via reprograming of the transcriptome. New Phytologist , 219 (1), 230–245. doi: 10.1111/nph.15160
Rosvall, M., & Bergstrom, C. T. (2008). Maps of random walks on complex networks reveal community structure. Proceedings of the National Academy of Sciences of the United States of America , 105 (4), 1118–1123. doi: 10.1073/pnas.0706851105
Rottmann, W. H., Meilan, R., Sheppard, L. A., Brunner, A. M., Skinner, J. S., Ma, C., … Strauss, S. H. (2000). Diverse effects of overexpression of LEAFY and PTLF, a poplar (Populus) homolog of LEAFY/FLORICAULA, in transgenic poplar and Arabidopsis. Plant J ,22 (3), 235–245. doi: 10.1046/j.1365-313x.2000.00734.x
Royer, D. L., McElwain, J. C., Adams, J. M., & Wilf, P. (2008). Sensitivity of leaf size and shape to climate within Acer rubrumand Quercus kelloggii . New Phytologist , 179 (3), 808–817. doi: 10.1111/j.1469-8137.2008.02496.x
Royer, D. L., Wilf, P., Janesko, D. A., Kowalski, E. A., & Dilcher, D. L. (2005). Correlations of climate and plant ecology to leaf size and shape: potential proxies for the fossil record. Am. J. Bot. ,92 (7), 1141–1151. doi: 10.3732/ajb.92.7.1141
Runions, A., Tsiantis, M., & Prusinkiewicz, P. (2017). A common developmental program can produce diverse leaf shapes. New Phytologist , 216 (2), 401–418. doi: 10.1111/nph.14449
Russin, W. A., & Evert, R. F. (1984). Studies on the Leaf of Populus deltoides (Salicaceae): Morphology and Anatomy. American Journal of Botany , 71 (10), 1398. doi: 10.2307/2443707
Schaefer, R., Michno, J.-M., Jeffers, J., Hoekenga, O. A., Dilkes, B. P., Baxter, I. R., & Myers, C. (2018). Integrating co-expression networks with GWAS to prioritize causal genes in maize. The Plant Cell , tpc.00299.2018. doi: 10.1105/tpc.18.00299
Schiffthaler, B., Serrano, A., Delhomme, N., & Street, N. R. (2018). Seidr: A toolkit for calculation of crowd networks. BioRxiv , 250696. doi: 10.1101/250696
Simons, Y. B., Bullaughey, K., Hudson, R. R., & Sella, G. (2018). A population genetic interpretation of GWAS findings for human quantitative traits. PLoS Biology , 16 (3). doi: 10.1371/journal.pbio.2002985
Sluis, A., & Hake, S. (2015). Organogenesis in plants: initiation and elaboration of leaves. Trends in Genetics : TIG , 31 (6), 300–306. doi: 10.1016/j.tig.2015.04.004
Souer, E., van Houwelingen, A., Kloos, D., Mol, J., & Koes, R. (1996). The no apical meristem gene of Petunia is required for pattern formation in embryos and flowers and is expressed at meristem and primordia boundaries. Cell , 85 (2), 159–170. doi: 10.1016/S0092-8674(00)81093-4
Stettler, R. F., & Bradshaw Jr, H. D. (1996). Evolution, Genetics, and Genetic Manipulation. In R. . Stettler, H. D. Bradshaw Jr, P. Heilman, & T. Hinckley (Eds.), Biology of Populus and its implications for management and conservation (pp. 1–6). Ottawa: NRC Research Press.
Street, N. R., Jansson, S., & Hvidsten, T. R. (2011). A systems biology model of the regulatory network in Populus leaves reveals interacting regulators and conserved regulation. BMC Plant Biology ,11 (1), 13. doi: 10.1186/1471-2229-11-13
Street, N. R., Sjödin, A., Bylesjö, M., Gustafsson, P., Trygg, J., & Jansson, S. (2008). A cross-species transcriptomics approach to identify genes involved in leaf development. BMC Genomics , 9 (1), 589. doi: 10.1186/1471-2164-9-589
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., … Mesirov, J. P. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences of the United States of America , 102 (43), 15545–15550. doi: 10.1073/pnas.0506580102
Sundell, D., Mannapperuma, C., Netotea, S., Delhomme, N., Lin, Y.-C. Y.-C., Sjödin, A., … Street, N. R. N. R. (2015). The Plant Genome Integrative Explorer Resource: PlantGenIE.org. New Phytologist ,208 (4), 1149–1156. doi: 10.1111/nph.13557
Sundell, D., Street, N. R. N. R., Kumar, M., Mellerowicz, E. J. E. J., Kucukoglu, M., Johnsson, C., … Hvidsten, T. R. (2017). AspWood: High-Spatial-Resolution Transcriptome Profiles Reveal Uncharacterized Modularity of Wood Formation in Populus tremula. The Plant Cell ,29 (7), 1585–1604. doi: 10.1105/tpc.17.00153
Traiser, C., Klotz, S., Uhl, D., & Mosbrugger, V. (2005). Environmental signals from leaves - a physiognomic analysis of European vegetation.New Phytologist , 166 (2), 465–484. doi: 10.1111/j.1469-8137.2005.01316.x
Tsuge, T., Tsukaya, H., & Uchimiya, H. (1996). Two independent and polarized processes of cell elongation regulate leaf blade expansion in Arabidopsis thaliana (L.) Heynh. Development , 122 (5), 1589–1600.
Tsukaya, H. (2005). Leaf shape: genetic controls and environmental factors. The International Journal of Developmental Biology ,49 (5–6), 547–555. doi: 10.1387/ijdb.041921ht
Vroemen, C. W., Mordhorst, A. P., Albrecht, C., Kwaaitaal, M. A. C. J., & de Vries, S. C. (2003). The CUP-SHAPED COTYLEDON3 gene is required for boundary and shoot meristem formation in Arabidopsis. The Plant Cell , 15 (7), 1563–1577. doi: 10.1105/TPC.012203
Wang, J., Ding, J., Tan, B., Robinson, K. M. K. M., Michelson, I. H. I. H., Johansson, A., … Ingvarsson, P. K. P. K. (2018). A major locus controls local adaptation and adaptive life history variation in a perennial plant. Genome Biology , 19 (1), 72. doi: 10.1186/s13059-018-1444-y
Wu, R., Bradshaw, H. D. H., Stettler, R. F. R., Jr., & Stettler, R. F. R. (1997). Molecular genetics of growth and development in Populus (Salicaceae). v. mapping quantitative trait loci affecting leaf variation. American Journal of Botany , 84 (2), 143–153. doi: 10.2307/2446076
Xiao, D., Wang, H., Basnet, R. K., Zhao, J., Lin, K., Hou, X., & Bonnema, G. (2014). Genetic Dissection of Leaf Development in Brassica rapa Using a Genetical Genomics Approach. PLANT PHYSIOLOGY ,164 (3), 1309–1325. doi: 10.1104/pp.113.227348
Yang, J., Weedon, M. N., Purcell, S., Lettre, G., Estrada, K., Willer, C. J., … Visscher, P. M. (2011). Genomic inflation factors under polygenic inheritance. European Journal of Human Genetics ,19 (7), 807–812. doi: 10.1038/ejhg.2011.39
Zhong, R., Allen, J. D., Xiao, G., & Xie, Y. (2014). Ensemble-Based Network Aggregation Improves the Accuracy of Gene Network Reconstruction. PLoS ONE , 9 (11), e106319. doi: 10.1371/journal.pone.0106319
Zhou, X., Carbonetto, P., & Stephens, M. (2013). Polygenic Modeling with Bayesian Sparse Linear Mixed Models. PLoS Genetics ,9 (2), e1003264. doi: 10.1371/journal.pgen.1003264
Zhou, X., & Stephens, M. (2014). Efficient multivariate linear mixed model algorithms for genome-wide association studies. Nature Methods , 11 (4), 407–409. doi: 10.1038/nmeth.2848

Data Accessibility Statement

The data are deposited in the European Nucleotide Archive (ENA) as accession PRJEB31491.

Author contributions

Planned and designed the research: NRS
Performed experiments: NM, BS, KMR, BKT
Contributed resources: NRS, SJ
Contributed analysis ideas: MV, MESB, TRH,
Analyzed data: NM, BS, KMR, BKT, CM
Wrote the paper: NRS with contributions from all authors
Table 1 . Broad sense heritability (H2), population differentiation (QST), Genotype by Environment interaction (GxE) and environmental correlations for leaf circularity, indent depth and area. Phenotypes were measured in two years (2008 and 2011) and two common gardens (Ekebo in southern Sweden and Sävar in northern Sweden). GxE interactions for each phenotype were tested using an ANOVA model with garden, year and garden x year as independent variables, and phenotype as the dependent variable. TheF ratio and P -value for the garden:year interaction are reported. Figure S1 details comparisons for each year/garden combination.