Figure 1. Populus tremula leaf physiology and morphology. A) A representative mature terminal leaf.B) A representative flushed short shoot bud with mature pre-formed leaves. C) Cross-section (15 µm thick) of a representative mature terminal leaf. D) Cross-section (10 µm thick) of a representative mature pre-formed leaf. Cross sections were embedded in 1,5 % plant agarose, sectioned using a vibratome and stained for one minute with Toluidine Blue; Magnification 200X.
Figure 2 . An overview of leaf circularity, indent depth and leaf area in the Swedish Aspen (SwAsp) collection. A) Example leaves from the SwAsp collection. Each leaf is from one genotype and the set of leaves was selected to reflect the extent of variation represented among all genotypes. B) Genetic correlations among measures of leaf circularity, median indent depth and leaf area for data recorded in two common gardens (Ekebo, southern Sweden and Sävar, northern Sweden) and two years (2008 and 2011). C) Heatmap representations of leaf circularity (left), median indent depth (centre) and leaf area (right) in the two common gardens and years. Each represented value is a median calculated from 10 leaves per clonal replicate (median n=3) of each genotype. In each heat map, genotypes are clustered by Euclidean distance with the population from which the genotype originates indicated by coloured bars to the left of the heat map. The values represented were used to perform the genome-wide association mapping results depicted in D. D) Intersection of the top 1,000 ranked SNPs from genome wide association mapping for the three traits measured in the two years and gardens.
Figure 3. The genomic context distribution of SNPs within the top ranked 1,000 genes identified using association mapping for leaf circularity and indent depth.
Figure 4. Selection of a set of high and low phenotype value genotypes. Genotypes with high and low phenotypic values of circularity (A ) and indent depth (B ) were selected on the basis of the top or bottom quartiles of the Best Linear Unbiased Prediction (BLUP) values for each trait. The colour intensity reflects the four quartiles of trait values with darker colours meaning higher values. The left density plots represent BLUPs for the two traits. The distribution plots in the centre depict the original trait measurements in the two gardens and the two years of sampling prior to the BLUP calculation. Gene expression assayed in young leaf buds (Mähler et al. 2017) was used to test for differential expression between genotypes within the low and high sets for each phenotypic trait. The Venn diagram shows the intersection between the significantly differentially expressed genes and the top 1,000 genes from the genome wide association study (GWAS) for the two traits.
Figure 5. Overview of gene expression in developmental series of terminal leaves. A) Principal component analysis of the 1,000 most variable genes. Data points are shaded by tissue type (Pf = Pre-formed; T = terminal) and leaf number in addition to different point styles being used. Ellipses around each group indicate a 95% confidence area for a particular sample group. B) Differential expression matrix of all sample points showing all significantly up-regulated genes (P adj < 0.05) in the lower triangular matrix and all significantly down-regulated genes (P adj < 0.05) in the upper triangular matrix. In the lower triangular, columns constitute the numerator of the differential test, while rows constitute the numerator in the upper triangular half. T = terminal.
Figure 6. Visualisation of differential expression within the terminal developmental series gene co-expression network. A)Transition between stage -4 and -3. B) Transition between stage -3 and -2. C) Transition between stage -2 and -1. D)Transition between stage -1 and 0. E) Transition between stage 0 and 1. F) Transition between stage 1 and 2. Nodes are sized inversely by adjusted P -value of a differential expression comparison at the relevant transition and shaded by the fold change (higher in red, lower in blue) with the earlier leaf number as the numerator. P- values ranged from 1.519197e-104 to 9.999962e-01 with the highest P -value nodes being 100x smaller than the lowest ones.
Figure 7. Developmental expression of genes differentially expressed between phenotypic extremes. The heatmap shows the expression of the differentially expressed genes depicted in Figure 4. Expression in the terminal leaf development data set is represented for genes differentially expressed between phenotypically extreme genotype sets for circularity (A ) and indent depth (B ). Genes that were found among the top 1,000 GWAS genes are indicated by coloured lines on the left side of the heatmap. Genes are clustered hierarchical clustering using average linkage. Samples are arranged by leaf age. Values depicted are VST (Variance-Stabilising Transformation) normalised expression. Clustering was performed using complete linkage and correlation-based distance.
Supplementary Figure S1 . Boxplots showing genotype x environment interactions for size and shape metrics in pre-formed leaves of the SwAsp collection in common gardens at Ekebo (56 °N) and Sävar (62 °N). Each box represents a common garden in a given year with 111 to 113 genotypes in clonal replicates. (A) Leaf circularity, (B) indent depth, and (C) leaf area. GxE analysis of variance F and P -values are provided in Table 1 Boxes annotated with the same letter are not significantly different at P < 0.05.
Supplementary Figure S2 . Manhattan and QQ-plots for leaf area, circularity, and indent depth.
Supplementary Figure S3 . Venn diagram of the top 1,000 genes from the genome wide association analysis for the three leaf traits.
Supplementary Figure S4 . Correlation between gene expression among genotypes of the Swedish Aspen collection and leaf shape and size traits. The left panel shows the distribution for all expressed genes for each phenotypic trait. The right panel separates genes that were differentially expressed between sets of genotypes with contrasting trait values from non-differentially expressed genes.
Supplementary Figure S5. Network clusters of genes within the community expression network inferred for the terminal leaf developmental series. Clusters were defined using InfoMap. The corresponding eigengene expression profile is shown for each partition as well as representative enriched PFAM and GO categories. * indicates an adjusted P -value for the terms between 0.05 and 0.01, ** indicates 0.01 - 0.001, *** indicates < 0.001. Only terms in the same range of significance were included for a cluster (if any). The node size is scaled according to the node’s PageRank.
Supplementary Figure S6 . Gene set enrichment analysis results for 1,000 top genes from the genome wide association (GWAS) and differentially expressed gene sets from the GWAS population (DEGs) within the terminal developmental gene co-expression network.
Supplementary Figure S7 . Heatmap of the scaled and centred gene expression values from the terminal leaf development series for all genes that were differentially expressed between genotype extreme groups for at least one of the traits, indicated by the green marks on the left.
Supplementary Figure S8 . Venn diagram of the genes differentially expressed between phenotypic extremes (high and low BLUP quartiles) of the three leaf traits.
Supplementary Figure S9. Scaled network centrality for genes that were differentially expressed between sets of genotypes with contrasting trait values and non-differentially expressed genes within the Swedish Aspen, terminal and pre-formed leaf development series co-expression networks. The Swedish Aspen co-expression network is from Mähler et al . (2017) and uses WGCNA-calculated network centrality scores. The leaf development network uses InfoMap scores. Both scores were scaled to uniform scale to enable comparison. Results were compared to randomised sets of genes to ensure that the observed signal was not a general pattern. Genes that were also among the top 1,000 genes identified by genome wide association mapping are indicated as coloured data points.
Supplementary Figure S10. Gene set enrichment analysis results for genome wide association and population differentially expressed gene sets within the Swedish Aspen gene co-expression network from Mähleret al. (2017).
Supplementary Figure S11. Tajima’s D values for the top 1,000 genes identified by genome wide association mapping and for genes that were differentially expressed between sets of genotypes with contrasting trait values. Tajima’s D values were obtained from Mähler et al.(2017). Differentially expressed genes that were also among the top 1,000 genes identified by genome wide association mapping are indicated as coloured data points.
Supplementary Table S1 . Top 1,000 genome wide associations for leaf area, circularity, and indent depth. The columns are snp: SNP name; chr: scaffold where SNP is located; pos: base-pair position of SNP on scaffold; pseudo_chr: pseudo-chromosome where SNP is located; pseudo_pos: base-pair position of SNP on pseudo-chromosome; minor_allele: minor allele of SNP; major_allele: major allele of SNP; maf: minor allele frequency of SNP; beta: effect size of the association; se: standard error of the beta; pvalue: likelihood ratio test P -value of the beta being different from zero; gene: gene that the SNP is located in/near (within 2 kbp), if any; feature: gene feature that the SNP is located in, if any; at_gene: best BLAST hit inArabidopsis thaliana ; potri_gene: best BLAST hit inPopulus trichocarpa .
Supplementary Table S2 . Top 1,000 genes from the genome wide associations for leaf area, circularity, and indent depth. The genes were selected by rank ordering the GWAS results by associationP -value and then picking genes from the top of the list until 1,000 unique genes had been selected.
Supplementary Table S3 . Differential expression results for the three leaf traits. In addition to the statistics from DESeq2, it has been indicated whether the gene is also an eGene in the SwAsp collection, whether the gene is part of the leading-edge subset in the gene set enrichment analysis for the trait, gene annotations, reciprocal best DIAMOND hit in A. thaliana , as well as the best DIAMOND hit in A. thaliana .
Supplementary Table S4 . Node centrality statistics of all genes in the terminal leaf network. The centrality statistics computed are PageRank, closeness, betweenness, strength, eigenvector, Katz, Kleinberg’s authority and Kleinberg’s hub score.
Supplementary Table S5 . Gene ontology and PFAM (protein family) enrichments for all clusters with more than five members in the terminal leaf network.
Supplementary Table S6 . Normalised and raw expression values and results of all pairwise differential expression comparisons containing mean expression, log2 fold changes, standard error of the log2 fold changes, the test statistic, the test P value and the adjusted P value for all genes in all comparisons.