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.