ORIGINAL RESEARCH article
Morphometrics Reveals Complex and Heritable Apple Leaf Shapes
- 1Department of Plant, Food and Environmental Sciences, Faculty of Agriculture, Dalhousie University, Truro, NS, Canada
- 2Donald Danforth Plant Science Center, St. Louis, MO, United States
- 3Department of Horticulture, Michigan State University, East Lansing, MI, United States
- 4Computational Mathematics, Science and Engineering, Michigan State University, East Lansing, MI, United States
Apple (Malus spp.) is a widely grown and valuable fruit crop. Leaf shape is important for flowering in apple and may also be an early indicator for other agriculturally valuable traits. We examined 9,000 leaves from 869 unique apple accessions using linear measurements and comprehensive morphometric techniques. We identified allometric variation as the result of differing length-to-width aspect ratios between accessions and species of apple. The allometric variation was due to variation in the width of the leaf blade, not the length. Aspect ratio was highly correlated with the first principal component (PC1) of morphometric variation quantified using elliptical Fourier descriptors (EFDs) and persistent homology (PH). While the primary source of variation was aspect ratio, subsequent PCs corresponded to complex shape variation not captured by linear measurements. After linking the morphometric information with over 122,000 genome-wide single nucleotide polymorphisms (SNPs), we found high SNP heritability values even at later PCs, indicating that comprehensive morphometrics can capture complex, heritable phenotypes. Thus, techniques such as EFDs and PH are capturing heritable biological variation that would be missed using linear measurements alone.
Apples (Malus spp.) are one of the world's most widely grown fruit crops, with the third highest global production quantity of over 84 million tons in 2014 (Food Agriculture Organization of the United Nations, 2017). The shape and size of apple leaves can impact access to light, ultimately influencing flower number and fruit quality (Wünsche and Lakso, 2000; Dennis, 2003). Apple leaves are generally simple, with an elliptical-to-ovate shape. Previous studies in apple used linear measurements, such as length and width, to quantify leaf shape (Liebhard et al., 2003; Bassett et al., 2011). The length-to-width aspect ratio is a major source of variation in leaf shape. Differing aspect ratios lead to a disproportionate increase or decrease in length relative to width, resulting in allometric variation, in leaves (Gurevitch, 1992; Chitwood et al., 2013). While linear measurements such as leaf length and width are useful, they fail to capture the full extent of leaf shape diversity. Failing to measure leaf shape comprehensively also limits our ability to discern the total underlying genetic contributions.
Eigenshape analysis is a traditional morphometric technique for comprehensively quantifying the outline of a shape, using a number of semilandmarks placed at equal distances around the shape (Lohmann, 1983). In comparison, Elliptical Fourier descriptors (EFD) analysis first converts the outline to a chaincode, a lossless data compression method that encodes shape by a chain of numbers, in which each number indicates step-by-step movements to reconstruct the pixels comprising the shape. A Fourier decomposition is subsequently applied to the chaincode, quantifying the shape as a harmonic series. EFDs have been used extensively to quantify leaf shape in diverse species, such as grape (Chitwood et al., 2014), tomato (Chitwood et al., 2012), and Passiflora (Chitwood and Otoni, 2017). Previous work used EFDs to assess apple fruit shape (Currie et al., 2000), but this technique has not yet been applied to apple leaves. A newly developed morphometric technique, persistent homology (PH), provides another method for estimating leaf shape. PH, like EFDs, is normalized to differences in size, but it is also orientation invariant. PH treats the pixels of a contour as a 2D point cloud before applying a neighbor density estimator to each pixel. A series of annulus kernels of increasing radii is used to isolate and smooth the contour densities. The number of connected components is recorded as a function of density for each annulus, resulting in a curve (a reduced version of persistent barcode) that quantifies shape as topology. The topology-based PH approach can also be applied to serrations and root architecture, allowing a similar framework to be used across disparate plant shapes and structures (Li et al., 2017a,b,c).
Comprehensively measuring leaf shape, using approaches such as EFDs and PH, may reveal shape features associated with fruit shape or quality traits. For example, recent work examining persimmon (Diospyros kaki) found weak but significant correlations between EFD measurements for leaf and fruit shape (Maeda et al., 2018). Similarly, in previous work on apple, several leaf traits such as area and perimeter were correlated with fruit size (Khan et al., 2014). There are also several cases of unique leaf characteristics providing an early marker for other genetic differences in apple. For example, the gene underlying red fruit flesh color may lead to anthocyanin accumulation in the leaves, causing red foliage (Chagné et al., 2007; Espley et al., 2009) while columnar tree architecture may be accompanied by an increase in leaf number, area, weight per unit area and length-to-width ratio (Talwara et al., 2013). Leaf pH has also been proposed as an early indicator of low acid fruit (Visser and Verhaegh, 1978).
In addition to serving as early markers for other traits, leaf shape and size may influence the amount of light a tree receives, and light exposure is crucial for flowering in apple. Light penetration results in higher levels of flowering, while leaf injury or defoliation can reduce flowering (Dennis, 2003). Thinning apple trees to a particular leaf-to-fruit ratio is a common practice to attain optimal fruit color and size. Contrastingly, trees with fewer fruit may increase vegetative growth and thus leaf area (Wünsche et al., 2000). Clearly, there is an important relationship between the leaves and the fruit of an apple tree, and comprehensively quantifying the variation in leaf shape is a crucial component to understanding this relationship.
Leaves are the main photosynthetic organs of apple, but the genetic basis underlying their shape and size remains unknown. Although there are examples of a single locus controlling major variation in leaf shape (Kimura et al., 2008; Cartolano et al., 2015; Andres et al., 2017), in most instances leaf shape appears to be controlled by numerous small-effect loci (Langlade et al., 2005; Tian et al., 2011; Chitwood et al., 2013). There are limited examples of genomic analyses of leaf shape in apple. However, a previous bi-parental linkage mapping study found two suggestive quantitative trait loci for leaf size (Liebhard et al., 2003). Work by Khan et al. (2014) measured several leaf traits such as area, perimeter and circularity, in 158 apple accessions. The study linked these measurements with 901 single nucleotide polymorphisms (SNPs), but found no significant genotype-phenotype relationships. Thus far, efforts have not been made to estimate the genetic heritability of comprehensive morphometric leaf phenotypes, such as those described using EFDs and PH. It therefore remains unclear to what extent these methods are capturing biologically meaningful, heritable variation.
To fully understand the genetic basis of leaf shape, it is essential to include both linear and morphometric estimates of shape. Decreasing sequencing costs and access to a large and diverse germplasm collection allowed us to analyze approximately 9,000 leaves from over 800 unique accessions which we linked to over 122,000 genome-wide SNPs. We present the first comprehensive analysis of leaf shape in apple, revealing that both accessions and species show allometric variation due to differences in the width of the leaf blade. While the primary axis of variation in apple using EFDs and PH is due to this allometric variation, we find high SNP heritability values even at later principal components, indicating that comprehensive estimates of shape capture heritable variation which would be missed by linear estimates alone.
Materials and Methods
Apple trees in Kentville, Nova Scotia, Canada were budded onto M.9 rootstocks in spring 2012. In the fall, the trees were uprooted and kept in cold storage until spring 2013, when trees were planted in an incomplete block design (see “REstricted Maximum Likelihood (REML)” below). Leaves from over 900 trees were collected from August 24th to September 16th 2015. Ten mature leaves were collected from each tree and efforts were made to collect leaves from portions of the branch representing third year growth. No leaves were collected from trees which had been heavily pruned due to the presence of disease. Leaves were flattened and placed to avoid touching, then scanned using Canon CanoScan (LiDE 220) Colour Image Scanners. Leaves were then dried for 48 h at 65°C and weighed to estimate the total dry weight (g) for each tree.
Leaf scans were converted into a separate binary image for each leaf using custom ImageJ macros, which included the “make binary” function (Abràmoff et al., 2004). Images were converted to RGB .bmp files and a chain code analysis was performed using SHAPE (Iwata and Ukai, 2002). The chain code was used to calculate normalized elliptical Fourier descriptors (EFDs) in SHAPE. The normalized EFDs were read into Momocs v1.1.5 (Bonhomme et al., 2014) in R (R Core Team, 2016) where harmonics B and C were removed to eliminate asymmetrical variation in leaf shape.
The binary leaf images were also analyzed using PH (Li et al., 2017c). To numerically estimate the shape of the leaves using PH, we extracted the leaf contour using a 2D point cloud (Figure 1A). After centering and normalizing the contour to its centroid size, we used a Gaussian density estimator (Figure 1B), which assigns high values (red) to pixels with many neighboring pixels, and low values (blue) to pixels with fewer neighboring pixels. We multiplied the density estimator by an annulus kernel, or ring (Figure 1C), which emphasizes the shape in an annulus at the centroid and is thus invariant to orientation (Figure 1D). The resulting function can also be visualized from the side view (Figures 1E,F). As we moved a plane from the top to the bottom, we recorded the number of connected components above the plane, forming a curve. With each new component this value increased, and each time components were merged, it decreased (Figure 1G). For each leaf we computed 16 curves corresponding to 16 expanding rings. For computational purposes, each curve was divided into 500 numbers, ultimately resulting in the shape of each leaf being represented by 8,000 (16*500) values.
Figure 1. Visualization of persistent homology technique for annulus kernel 7. Binary images were converted into a 2D point cloud (A) which was then normalized using a Gaussian density estimator (B). For each leaf, 16 annulus kernels were used. Annulus kernel 7, indicated in purple (C) is used as an example for this visualization. The density estimator is multiplied by ring 7 (D). The function can also be visualized from the side view (E,F). As a plane moves from top to bottom, the number of connected components is recorded along the curve (G). Below (G) are five visualizations of curves that are represented as red vertical dotted lines in (G).
Only leaves for which both EFDs and PH shape estimations were successfully calculated were included in subsequent analyses. Additionally, only trees with 8–10 leaves were included, as leaves were sometimes removed due to tears, folding, or the absence of a petiole, which did not allow for accurate quantification of shape. The final dataset consisted of 915 trees with 8–10 leaves, which included 869 unique accessions and 8,995 leaves.
EFDs and PH values were averaged across leaves from an individual tree. The contribution of EFD harmonics 1–15 to the mean leaf shape across all trees was visualized using the “hcontrib” function in the Momocs R package (Figure 2). To allow for discrimination between accessions based on leaf shape, principal component analysis (PCA) was performed using the Momocs “PCA” function (Bonhomme et al., 2014) for EFDs, and the “prcomp” function in R for PH values, which center but do not scale the data. The resulting principal component (PC) values were adjusted using REstricted Maximum Likelihood (see below). Subsequently, we identified the accession with the minimum and maximum value along each of the first 5 PCs.
Figure 2. Contribution of elliptical Fourier descriptor harmonics to leaf shape. The leaf shapes depicted are the mean leaf shapes based on all 915 trees. Harmonics 1–15 are represented on the x-axis and each harmonic is multiplied by the amplification factor on the y-axis to visualize their contribution to mean leaf shape. An amplification factor of 0 indicates the removal of the harmonic; a factor of 1 results in the normal shape; and values above 1 exaggerate effects to better visualize the harmonic's contribution to the final shape.
In addition to estimating the contour of the leaf using EFDs and PH, we used several more metrics to describe the leaves. Using ImageJ, we automated the measurement of leaf surface area (cm2), length (cm) of the leaf and width (cm) of the leaf as well as major (blade length) and minor (blade width) axes of the best fitting ellipse—which excluded the petiole—through batch processes (Abràmoff et al., 2004). Throughout the manuscript, we use “major” when referring to the length of the leaf blade, and “minor” when referencing the width of the leaf blade. We also calculated the aspect ratio of the leaf, by dividing the major axis by the minor axis. Additionally, leaf mass per area was calculated for 780 trees where we possessed surface area data for all 10 leaves, by calculating the ratio of dry weight to surface area (g/cm2).
While linear phenotypes were calculated as an average value for a particular tree, we also estimated variance within a tree for aspect ratio, length, width, major and minor axis, and surface area. Variance was calculated as the coefficient of variation using the “cv” function in the raster package (Hijmans, 2016) in R to estimate within-tree variability in leaf size, which is indicated as “var” throughout this manuscript.
REstricted Maximum Likelihood (REML) Adjustment of Phenotype Data
The orchard sampled in this study is an incomplete block design with 1 of 3 standards per grid. The standards, or “control trees”—‘Honeycrisp,’ ‘SweeTango,’ and ‘Ambrosia’—are replicated across the grid. Leaves from these trees were sampled multiple times across the orchard which allowed us to correct for positional effects. Each phenotype was adjusted using a REstricted Maximum Likelihood (REML) model which resulted in one adjusted value per accession, even when multiple trees were measured. The impact of row grid (rGrid), column grid (cGrid), and rGrid × cGrid effects were adjusted for using the following REML model:
We fit a linear mixed-effects model via REML using the “lmer” function in the lme4 package in R (Bates et al., 2015) and then calculated the least squares means using the “lsmeans” function in the lsmeans R package (Lenth, 2016).
Thus, while the initial phenotype data were collected for 915 trees, following REML adjustment, one value remained per unique accession, resulting in 869 accessions. REML-adjustment was applied directly to all size, weight and variance estimates. For PH and EFDs, we applied the REML following PCA and thus the percent contribution for each PC was calculated using unadjusted values. The adjusted data for all 24 phenotypes are included in Table S1.
The correlation between leaf phenotypes was calculated using Pearson's correlation and p-values were Bonferroni-corrected for multiple comparisons. The resulting heatmap was visualized using the “geom_tile” function in ggplot2 in R (Wickham, 2009). Next, we examined the leaves for allometry using the “SMA” function in the smartr R package (Warton et al., 2012) to estimate if the slope between the log-transformed minor and major axis differed from 1.
Accessions were labeled as either Malus × domestica Borkh. or Malus sieversii Lebed. based on information provided by the United States Department of Agriculture (USDA) Germplasm Resources Information Network website (https://npgsweb.ars-grin.gov) (Table S2). We used a Mann-Whitney U-test to test if any phenotypes differed between species and Bonferroni-corrected all p-values for multiple comparisons.
DNA was extracted using commercial extraction kits. Genotyping-by-sequencing (GBS) libraries were prepared using ApeKI and PstI-EcoT221I restriction enzymes according to Elshire et al. (2011). GBS libraries were sequenced using Illumina Hi-Seq 2000 technology. Reads which failed Illumina's “chastity filter” were removed from raw fastq files. Remaining reads were aligned to the Malus × domestica v1.0 pseudo haplotype reference sequence (Velasco et al., 2010) using the Burrows-Wheeler aligner tool v0.7.12 (Li and Durbin, 2009) and the Tassel version 5 pipeline (Glaubitz et al., 2014). Tassel parameters included a minKmerL of 30, mnQS of 20, mxKmerNum of 50000000 and batchSize of 20. The kmerlength was set to 82 for ApeKI and 89 for PstI-EcoT22I based on the max barcode size. The minMAF for the DiscoverySNPCallerPluginV2 was set to 0.01. All other default parameters were used. Non-biallelic sites and indels were removed using VCFtools v.0.1.14 (Danecek et al., 2011). Variant Call Format (VCF) files for both enzymes were then merged using a custom perl script, preferentially keeping SNPs called by PstI-EcoT22I at overlapping sites, since those sites tended to have higher coverage.
Missing genotypes were imputed using LinkImputeR v0.9 (Money et al., 2017) with global thresholds of 0.01 for minor allele frequency (MAF) and 0.70 for missingness. We examined depths of 3–8 and selected a case for imputation with a max position/sample missingness of 0.70, a minimum depth of 5, and an imputation accuracy of 94.9%. The VCF was converted to a genotype table using PLINK v1.07 (Purcell et al., 2007; Purcell, 2009).
Of the 869 accessions assessed in this study, 816 had genomic data following imputation and filtering and were included in downstream analyses. The resulting genotype table consisted of 816 accessions and 197,565 SNPs. Subsequently, a 0.05 MAF filter was applied using PLINK, after which 128,132 SNPs remained. SNPs with more than 90% heterozygous genotypes were removed. The final genotype table consisted of 816 samples and 122,596 SNPs.
Prior to performing PCA, SNPs were pruned for linkage disequilibrium (LD) using PLINK. We considered a window of 10 SNPs, removing one SNP from a pair if R2 > 0.5, then shifting the window by 3 SNPs and repeating (PLINK command: indep-pairwise 10 3 0.5). This resulted in a set of 75,973 SNPs for 816 accessions. PCA was performed on the LD-pruned genome-wide SNPs using “prcomp” in R. The first 2 genomic PCs were visualized using ggplot2 in R (Wickham, 2009).
We performed a genome-wide association study (GWAS) using the mixed linear model in Tassel (version 5) for each phenotype, adjusting for relatedness among individuals using a kinship matrix as well as the first 3 PCs for population structure (Bradbury et al., 2007; Zhang et al., 2010). The threshold for significance was calculated using simpleM (Gao et al., 2008, 2010) which estimates the number of PCs needed to explain 0.995 of the variance, or the number of independent SNPs. The inferred Meff used to calculate the significance threshold was 91,667 SNPs.
We searched the regions surrounding any significant GWAS SNPs using the Genome Database for Rosaceae GBrowse tool for Malus × domestica v1.0 p genome (Jung et al., 2014). We used a window of ±5,000 bp (10 kb) surrounding the significant SNP to check for genes, and when identified, we used the basic local alignment search tool (BLAST) from NCBI to search for the mRNA sequence and reported the result with the max score (Altschul et al., 1990).
Genomic prediction was performed using the “x.val” function in the R package PopVar (Mohammadi et al., 2015). The rrBLUP model was selected and 5-fold (nFold = 5) cross-validation was repeated 3 times (nFold.reps = 3) with no further filtering (min.maf = 0) from the set of 122,596 SNPs used for GWAS. All other default parameters were used. In addition to performing genomic prediction on the main 24 phenotypes examined in this study, we performed genomic prediction on all 40 PCs for EFDs and on the first 40 PCs for PH values. We also used the “rnorm” function in R to generate 1,000 random phenotypes with a mean of 0 and a standard deviation of 1, and performed genomic prediction using these random phenotypes to obtain the range of genomic prediction accuracies one can expect at random. Lastly, we used genome-wide complex trait analysis (GCTA) v.1.26.0 which estimates the genetic relationships between individuals based on genome-wide SNPs and uses this information to calculate the variance explained by these SNPs. The ratio of additive genetic variation to phenotypic variance is used to calculate SNP heritability (h2) of a trait (Yang et al., 2011). We used GCTA to estimate heritability for each phenotype, including the first 40 PCs for EFD and PH. We also estimated the correlation between genomic prediction accuracy (r) and SNP heritability (h2) using a Pearson's correlation.
Variation in Apple Leaf Shape
We examined 24 phenotypes related to apple leaf shape and size including length, width, surface area, dry weight, leaf mass per area, within-tree variance, and overall shape estimated using PCs derived from EFD and PH data (see Materials and Methods and Figures 1–2). The distribution and raw values for each phenotype are provided (Figure S1, Table S1).
To visualize the primary axes of morphometric variation, we chose a representative leaf from accessions with the minimum and maximum values along the first 5 PCs for EFDs and PH (Figure 3A). The accessions with extreme values along PC1 for both methods are similar. In fact, ‘Binet Rouge’ has the lowest value along PC1 for EFD and PH, with the axis clearly representing a decrease in the length-to-width (aspect) ratio. The annulus kernels most strongly contributing to PH PC1 (Figure S2) provide further evidence that this PC captures variation in aspect ratio. Variation in leaf shape captured by higher-order PCs is more complex and cryptic, and is thus not captured using linear measurements alone. In addition, while the primary axis of variation (PC1) using EFDs and PH may explain similar aspects of leaf morphology, the morphospaces resulting from the two techniques differ (Figure 3B).
Figure 3. Examples of leaf shape across PCs derived from EFDs and PH. Binary images of leaves from accessions with minimum and maximum values along PCs 1–5 for EFD and PH estimates. PCs were calculated using averages estimated across 8–10 leaves but only a single representative leaf is displayed. PCs were REML-adjusted based on tree position in the orchard. The accession name is also listed (A). Visualization of PC1 vs. PC2 for EFD and PH data. Accessions with minimum and maximum values along PC1 and PC2 are indicated (B).
Next, we examined the correlation between all measured traits (Table S3). By assessing the correlation of PCs resulting from a classical morphometric technique such as EFDs with a novel, topology-based morphometric approach like PH, we reveal how complementary the methods are (Figure 4, Figure S3). While there is a highly significant correlation between PC1 for both methods (R2 = 0.949, p < 1 × 10−15), later PCs are often not significantly correlated, with the most notable exception being EFD PC2 and PH PC3 (R2 = 0.432, p < 1 × 10−15), although several other PCs also show weak correlations. Thus, while the primary axis of variation (PC1) is consistent and highly correlated between methods, each method captures distinct aspects of leaf morphology in subsequent PCs.
Figure 4. Correlations among leaf phenotypes. Values above the diagonal are colored according to the Pearson's correlation coefficient, and those below the diagonal indicate Bonferroni-corrected p-values. The box enclosed by the dotted lines indicates comparisons between phenotypes captured by comprehensive morphometric analyses.
Many of the leaf phenotypes show a strong correlation with each other (Figure 4). In particular, aspect ratio is highly correlated with PH PC1 (r = −0.878, p < 1 × 10−15), EFD PC1 (r = −0.855, p < 1 × 10−15) and minor axis (leaf blade width) (r = −0.734, p < 1 × 10−15). The correlation between the minor axis of a leaf and surface area (r = 0.939, p < 1 × 10−15) is higher than the correlation between the major axis (blade length) and surface area (r = 0.810, p < 1 × 10−15). As expected, leaf surface area is also highly correlated with average leaf dry weight (r = 0.934, p < 1 × 10−15), indicating that larger leaves are heavier.
Allometry in Apple Leaves
The high correlation between aspect ratio and PC1 for both EFD and PH methods indicates that length-to-width ratio is the primary source of variation in apple leaf shape. If there is an allometric relationship between the minor and major axis, and thus, the length and width of a leaf do not increase at equal rates, a slope significantly differing from 1 is expected. We find that the slope between the two measurements is significantly >1 (95% CI = 1.506–1.678, R2 = 0.343, p < 1 × 10−15), indicating that the minor axis increases at a greater rate than the major axis. While there is no significant correlation between the major axis (blade length) and EFD PC1 (R2 = 0.001, p = 1) or PH PC1 (R2 = 0.002, p = 1), there is a significant correlation for the minor axis (blade width) and EFD PC1 (R2 = 0.541, p < 1 × 10−15) and PH PC1 (R2 = 0.573, p < 1 × 10−15) (Figure 5). As PC1 explains 80.23% of the variation in the leaf shape for EFDs, and 62.20% for PH, it is apparent that the width of the leaf blade, and not length, is the major source of leaf shape variation in apple. In fact, the aspect ratio, calculated as the ratio of the major axis to minor axis, is even more strongly correlated with EFD and PH PC1, with an R2 of 0.732 for EFD PC1 (p < 1 × 10−15) and an R2 of 0.771 for PH PC1 (p < 1 × 10−15). Given the significant correlation between EFD PC1 and PH PC1 (Table S3, Figure S3), it is not surprising that aspect ratio is highly correlated with both.
Figure 5. Correlation between the primary axis of variation (PC1) captured using EFD and PH values and leaf shape measures. EFD PC1 is plotted against the major axis (length of leaf blade) (A), minor axis (width of leaf blade) (B) and aspect ratio (ratio of length-to-width of blade) (C). PH PC1 is plotted against the same measures in (D–F). The percent variances explained by PC1, prior to REML-adjustment, is shown in parentheses. All p-values are Bonferonni-corrected based on the number of comparisons in Figure 4. A regression line from a linear model with a shaded 95% confidence interval is also shown.
In addition to variation between accessions, we investigated differences in leaf shape and size between species by comparing M. domestica, the domesticated apple, with its primary progenitor species, M. sieversii (Table S4). PCA of the genome-wide SNP data reveals a primary axis of genetic variation that separates M. domestica and M. sieversii, although separation is incomplete (Figure 6A). The major axis (p = 0.975) of the leaves does not differ between species (Figure 6B). However, the minor axis (p = 4 × 10−4) of M. domestica leaves are significantly larger than M. sieversii (Figure 6C) and the aspect ratio (p = 0.023) is significantly less (Figure 6D). Thus, there is allometric variation both within (Figure 5) and between (Figure 6) Malus species.
Figure 6. Genetic and phenotypic comparison of the domesticated apple and its wild ancestor. PCs 1 and 2 were derived from 75,973 genome-wide SNPs and samples are labeled as M. domestica (purple), M. sieversii (green) or unknown (gray) (A). M. domestica leaves do not differ from M. sieversii leaves along the major axis (B), but they have a larger minor axis (C) and aspect ratio (D). P-values reported are Bonferroni-corrected based on multiple comparisons (Table S4). Species labels are based on USDA classification.
The Genetic Basis of Leaf Shape in Apple
GWAS of the 24 leaf phenotypes examined in this study yielded few significant results (Figure S4). We identified 70 significant SNPs representing 5 phenotypes which are reported in Table S5. Manhattan plots for 4 phenotypes which include 69 of the 70 significant SNPs are shown in Figure 7. We examined the regions surrounding significant SNPs for candidate genes using the GBrowse tool (Table S6; Jung et al., 2014). We searched within a ±5,000 bp window, which should capture any linked causal variation given the rapid LD decay observed in a diverse collection of apples that is largely replicated in the germplasm studied here (Migicovsky et al., 2016). However, no strong candidate genes were identified.
Figure 7. Manhattan plots of GWAS results for traits of interest. Results are shown for EFD PC3, PH PC4, leaf mass specific area, and width variance. P-values are log-transformed. The horizontal dotted line represents the simpleM-corrected significance threshold. Chromosome R indicates SNPs found on contigs unanchored to the reference genome. All remaining Manhattan plots are included in Figure S4.
While GWAS examines the genome for single, large-effect loci, genomic prediction estimates our ability to predict a phenotype using genome-wide marker data. We complimented our GWAS with genomic prediction and observed prediction accuracies (r) ranging from −0.10 to 0.52 (Table S7, Figure S5A). Aspect ratio is the primary source of variation in leaf shape (Figure 5C) and it was also the leaf measurement that had the highest genomic prediction accuracy (0.52). Other phenotypes highly correlated with aspect ratio, such as leaf width (0.51), minor axis (0.49), EFD PC1 (0.48), and PH PC1 (0.47), all had relatively high prediction accuracies. PH PC3 (0.51) was also among the most well-predicted using genetic data.
Similarly, estimates of SNP heritability (h2) calculated using GCTA (Yang et al., 2011) ranged from 0 to 0.75, with the highest heritability observed for aspect ratio (0.75) followed by leaf width (0.71), EFD PC1 (0.71), minor axis (0.69), and PH PC1 (0.65) (Figure 8, Table S8). Heritability estimates were highly correlated with genomic prediction accuracies (Figure S5B, R2 = 0.936, p < 1 × 10−15), which is not surprising given that both techniques involve predicting a phenotype from genome-wide SNP data. None of the phenotypes measuring variance within the 8–10 leaves sampled had heritability estimates significantly different from 0.
Figure 8. SNP heritability (h2) of leaf phenotypes. Values represent the additive genetic variance (Vg) divided by the phenotypic variance (Vp) with a standard error as calculated using GCTA (Yang et al., 2011). The dotted lines are found at h2 = 0, indicating that none of the phenotypic variance is explained by the genetic data. The proportion of the total phenotypic variance explained by each PC is indicated in parentheses.
While the primary axis of variation in leaf shape detected by EFDs and PH is aspect ratio, we were also interested in determining if higher-order PCs, which capture variation not readily visible to the eye, are extracting information that is biologically meaningful. Using genomic prediction and heritability estimates, we found evidence of a genetic basis for these complex “hidden phenotypes,” which are unmeasurable using linear techniques. For example, the heritability of phenotypes such as PH PC6 (0.48), PH PC9 (0.35), PH PC10 (0.33), and EFD PC9 (0.33) is similar to traditionally measured phenotypes such as leaf length (0.44) and leaf mass per area (0.40). While higher PCs may have relatively high heritability values, after a certain point the values (± standard error) overlap with 0, indicating that they are not heritable. The cutoff for morphometric PCs with a heritable genetic basis is approximately PC17. These results suggest that by making use of morphometric techniques that measure shape comprehensively, we are describing biologically meaningful, heritable phenotypes which would be missed by simple measurements such as leaf length, width and surface area.
Leaves play a crucial role in the growth and development of apple trees. To elucidate the genetic basis of variation in apple leaves, we quantified leaf shape using traditional linear measurements and comprehensive morphometric techniques. Our work offers the first comparison between the novel topology-based technique, PH, and EFDs, which we find are complementary but distinct methods. For both methods, PC1 was highly correlated with the aspect ratio, providing evidence that the primary axis of variation in apple leaf shape can be captured using linear measurements. The minor axis, or width of the leaf blade, was also highly correlated with PC1, while the major axis was not. Thus, variation in the aspect ratio is due to variation in the leaf blade width, not length. Leaf surface area was also more highly correlated with the minor axis than the major axis. Variation in leaf width is therefore essential to both the size and shape of apple leaves, similar to previous work in tomato (Schwarz and Kläring, 2001).
The width of the leaf blade is not only the source of variation between apple accessions, but also between M. domestica and M. sieversii. The presence of the same allometric relationship within and between species suggests that the genetic loci controlling intra-specific leaf shape variation within M. domestica may be the same as those controlling the divergence in leaf shape observed between the domesticated apple and its wild ancestor. For example, in birds, while PC1 and PC2 of bill shape explain the majority of variation across 2,000 species, they are also consistently associated with the variation between higher taxa (possessing >20 species) (Cooney et al., 2017). Our results suggest that the increase in leaf size since domestication has not been an overall increase in leaf size but specifically an increase in blade width leading to larger leaves with a reduced length-to-width ratio.
Our work provides evidence that allometry is the primary source of morphometric variation in apple leaves. These findings are consistent with work reported in other species such as tomato, where the length-to-width ratio was the major source of shape variation (>40%) (Chitwood et al., 2013). Similarly, work in Passiflora and Vitis species performed using two independent morphometric techniques identified allometric variation as the primary source of variation in PC1, which explained at least 40% of the variation in leaf shape (Chitwood and Otoni, 2017; Klein et al., 2017). Thus, linear measurements—in particular aspect ratio—are an important source of information when describing leaf shape. However, linear measurements are not sufficient for capturing the full spectrum of diversity. In our study, PC1 accounts for 62.20% or 80.23% of the variation, depending on the technique used. By simply quantifying apple leaves using linear measurements, we would miss nearly 40% of the variation in some cases. While PC1 is highly correlated with aspect ratio, later PCs represent orthogonal variation that can likely only be captured through morphometric techniques such as EFDs and PH. Therefore, to fully quantify variation in leaf shape, comprehensive morphometric techniques are essential.
To discern the genetic contributions to leaf shape, we paired both linear and comprehensive morphometric estimates of shape with genome-wide SNP data. There are examples of a simple genetic basis of leaf shape, such as in Arabidopsis thaliana, where the ANGUSTIFOLIA and ROTUNDIFOLIA3 independently control leaf width and length (Tsuge et al., 1996). In barley, transcript levels of BFL1 limit leaf width, with overexpression resulting in narrower leaves and loss of BFL1 function resulting in a reduced length-to-width ratio (Jöst et al., 2016). Using GWAS, we found no robust associations with shape phenotypes, observed a low ratio of significant SNPs to the number of phenotypes examined, and found that significant SNPs were sparsely distributed across multiple chromosomes. In addition, the small number of significant SNPs are likely spurious associations due to poor correction for cryptic relatedness, as evidenced by the QQ plots (Figure S4). These observations suggest that leaf shape is likely polygenic and controlled by a large number of small effect loci, such as in tomato and maize (Tian et al., 2011; Chitwood et al., 2013). In comparison, GWAS on apple fruit phenotypes, such as color and firmness, have revealed strong associations resulting from a small number of large effect loci (Migicovsky et al., 2016). However, it is possible that large effect loci were missed in the present study, either because of poor reference genome assembly or inadequate marker density due to rapid LD decay. Improvements in genome assembly and increases in marker number will aid to further reveal the genetic architecture of apple leaf shape variation.
Lastly, we investigated the degree to which leaf shape is heritable and can be predicted using genome-wide SNP data. We find that the genomic prediction accuracies of the primary axes of leaf shape variation are similar to previously reported estimates for fruit width (0.48) and length (0.47), indicating that leaf shape is as heritable as fruit shape (Migicovsky et al., 2016). In combination with few significant GWAS results, high SNP heritability estimates support a polygenic basis for leaf shape. Aspect ratio was identified as the primary source of variation in leaf shape in apple and had the highest genomic prediction and heritability estimates, indicating that there is a genetic, heritable basis for allometric variation in apple. Further, although the first 5 PCs for both EFDs and PH explain the majority of the variation in apple leaf shape, most PCs from 1 to 14 have heritability estimates above 0.20 and may still represent crucial differences in leaf shape from an ecological, evolutionary, or agricultural perspective. Thus, while our ability to detect the primary axes of variation in leaf shape using genome-wide data is expected, our observation that higher level PCs are also heritable confirms that these comprehensive morphometric methods capture biologically meaningful variation that would be missed by linear measurements alone.
It is clear from our work that variation in apple leaf shape and size are under genetic control. Further, high genomic prediction and heritability estimates for higher morphometric PCs indicate that techniques such as EFDs and PH are capturing heritable biological variation that will be missed if researchers restrict leaf shape estimates to linear measurements. Additionally, a better understanding of the variation in leaf shape and size in apple could have important implications for canopy management, where light exposure is crucial to flowering (Dennis, 2003). Ultimately, through the first in-depth study of leaf shape in apple, we uncover allometry between accessions and species, as well as evidence that complex and heritable phenotypes can be captured using comprehensive morphometric techniques.
ZM, DC, and SM designed the research. ZM and ML performed the research. ZM analyzed the data. DC and SM supervised the experiments. ZM, ML, DC, and SM wrote the article.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to acknowledge Gavin Douglas and Sherry Fillmore for their help setting up the statistical design of the orchard and SNP-calling pipeline. We also thank all past and present members of the Myles Lab for their work in maintaining the apple orchard. This research was supported in part, thanks to funding from the Canada Research Chairs program (SM), the National Sciences and Engineering Research Council of Canada (SM, ZM), Genome Canada (SM) and a Killam Predoctoral Scholarship from Dalhousie University (ZM). This research was published originally as part of the doctoral thesis of ZM (Migicovsky, 2017).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2017.02185/full#supplementary-material
Figure S1. Distribution of leaf phenotypes following REML-adjustment. N is equal to the total number of unique samples.
Figure S2. Visualization of contributions of each annulus kernel to PH PC1. Annulus kernels 6, 7 and 16 contribute the most to leaf shape according to PH PC1. The placement of each annulus kernel is visualized on a leaf representing the minimum and maximum value along PC1 (A). The contribution to PC1 of each of the 16 annulus kernels is also shown (B).
Figure S3. Comparison of morphometric EFD and PH PCs 1 to 5. Correlation between first 5 PCs, estimated using Pearson's correlation, including R2 and Bonferroni corrected p-values based on Figure 4, Table S3.
Figure S4. GWAS results for all 24 leaf phenotypes examined. Manhattan and QQ plots are included for each phenotype. The QQ-plot shows both the results of a naive GWAS (Pearson correlation) and the results from applying the mixed model. P-values are log-transformed and the threshold for significance is simpleM-corrected and indicated by a horizontal dotted line. Chromosome R indicates SNPs found on contigs unanchored to the reference genome.
Figure S5. Genomic prediction accuracy (r) (A) and correlation between genomic prediction results and SNP heritability estimates (h2) for all leaf phenotypes (B). Genomic prediction accuracies represent the average correlation (± standard deviation) between observed and predicted phenotype scores, based on 5-fold cross-validation with 3 iterations. Dotted red lines indicate the minimum and maximum prediction average accuracy (r) achieved using 1,000 randomly generated phenotypes. The percent variance explained by each PC was calculated prior to REML-adjustment and is indicated in parentheses.
Table S1. All leaf phenotypes assessed in apple, following REML-adjustment. Accessions are identified by their unique “apple id.” Species information for these accessions is available in Table S2.
Table S2. Species (M. domestica, M. sieversii, or unknown) for all accessions assessed in this study.
Table S3. Correlation between leaf phenotypes as well as Bonferroni-adjusted p-values. Pearson's product moment correlation coefficients are reported. These results are visualized in Figure 2.
Table S4. Comparison of leaf phenotypes between accessions based on species. Bonferroni-adjusted p-values resulting from a Mann-Whitney U-test estimating the difference between accessions based on species (M. domestica/M. sieversii) for the leaf phenotypes examined.
Table S5. Positional information for significant GWAS results. Additional information about significant SNPs are included such as p-value, marker R2, minor and major allele, minor and major effect and MAF.
Table S6. Genes found within ±5 kb of SNPs with significant GWAS associations. Results are listed according to the Genome Database for Rosaceae GBrowse (accessed January 27, 2017). Overlapping mRNA, length, contig, GO category, GO term accession, GO term name, InterPro Term, InterPro Description, and NCBI sequence with Max Score when BLASTed using NCBI are reported.
Table S7. Genomic prediction accuracies (r) for leaf phenotypes. r_avg represents the average correlation between observed and predicted phenotype scores, based on 5-fold cross-validation with 3 iterations. The standard deviation (r_sd) is also reported. These results are visualized in Figure S5A.
Table S8. SNP heritability (h2) for leaf phenotypes. h2 represents the genetic variance (Vg) divided by the phenotypic variance (Vp). The standard error (SE) is also reported. These results are visualized in Figure 8.
Andres, R. J., Coneva, V., Frank, M. H., Tuttle, J. R., Samayoa, L. F., Han, S. W., et al. (2017). Modifications to a LATE MERISTEM IDENTITY1 gene are responsible for the major leaf shapes of Upland cotton (Gossypium hirsutum L.). Proc. Natl. Acad. Sci. U.S.A. 114, E57–E66. doi: 10.1073/pnas.1613593114
Bassett, C. L., Glenn, D. M., Forsline, P. L., Wisniewski, M. E., and Farrell, R. E. (2011). Characterizing water use efficiency and water deficit responses in apple (Malus × domestica Borkh. and Malus sieversii Ledeb.) M. Roem. HortScience 46, 1079–1084.
Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308
Cartolano, M., Pieper, B., Lempe, J., Tattersall, A., Huijser, P., Tresch, A., et al. (2015). Heterochrony underpins natural variation in Cardamine hirsuta leaf form. Proc. Natl. Acad. Sci. U.S.A. 112, 10539–10544. doi: 10.1073/pnas.1419791112
Chagné, D., Carlisle, C. M., Blond, C., Volz, R. K., Whitworth, C. J., Oraguzie, N. C., et al. (2007). Mapping a candidate gene (MdMYB10) for red flesh and foliage colour in apple. BMC Genomics 8:212. doi: 10.1186/1471-2164-8-212
Chitwood, D. H., Headland, L. R., Kumar, R., Peng, J., Maloof, J. N., and Sinha, N. R. (2012). The developmental trajectory of leaflet morphology in wild tomato species. Plant Physiol. 158, 1230–1240. doi: 10.1104/pp.111.192518
Chitwood, D. H., Kumar, R., Headland, L. R., Ranjan, A., Covington, M. F., Ichihashi, Y., et al. (2013). A quantitative genetic basis for leaf morphology in a set of precisely defined tomato introgression lines. Plant Cell 25, 2465–2481. doi: 10.1105/tpc.113.112391
Chitwood, D. H., and Otoni, W. C. (2017). Morphometric analysis of Passiflora leaves: the relationship between landmarks of the vasculature and elliptical Fourier descriptors of the blade. Gigascience 6, 1–13. doi: 10.1093/gigascience/gix070
Chitwood, D. H., Ranjan, A., Martinez, C. C., Headland, L. R., Thiem, T., Kumar, R., et al. (2014). A modern ampelography: a genetic basis for leaf shape and venation patterning in grape. Plant Physiol. 164, 259–272. doi: 10.1104/pp.113.229708
Cooney, C. R., Bright, J. A., Capp, E. J. R., Chira, A. M., Hughes, E. C., Moody, C. J. A., et al. (2017). Mega-evolutionary dynamics of the adaptive radiation of birds. Nature 542, 344–347. doi: 10.1038/nature21074
Currie, A. J., Ganeshanandam, S., Noiton, D. A., Garrick, D., Shelbourne, C. J. A., and Oraguzie, N. (2000). Quantitative evaluation of apple (Malus × domestica Borkh.) fruit shape by principal component analysis of Fourier descriptors. Euphytica 111, 219–227. doi: 10.1023/A:1003862525814
Dennis, F. Jr. (2003). “Flowering, pollination and fruit set and development,” in Apples: Botany, Production and Uses, eds D. C. Ferree and I. J. Warrington (Wallingford, UK: CABI Publishing), 153–166.
Elshire, R. J., Glaubitz, J. C., Sun, Q., Poland, J. A., Kawamoto, K., Buckler, E. S., et al. (2011). A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE 6:e19379. doi: 10.1371/journal.pone.0019379
Espley, R. V., Brendolise, C., Chagné, D., Kutty-Amma, S., Green, S., Volz, R., et al. (2009). Multiple repeats of a promoter segment causes transcription factor autoregulation in red apples. Plant Cell 21, 168–183. doi: 10.1105/tpc.108.059329
Food and Agriculture Organization of the United Nations (2017). FAOSTAT. Available online at: http://www.fao.org/faostat/en/#data (Accessed Mar 9, 2017).
Gao, X., Becker, L. C., Becker, D. M., Starmer, J. D., and Province, M. A. (2010). Avoiding the high Bonferroni penalty in genome-wide association studies. Genet. Epidemiol. 34, 100–105. doi: 10.1002/gepi.20430
Gao, X., Starmer, J., and Martin, E. R. (2008). A multiple testing correction method for genetic association studies using correlated single nucleotide polymorphisms. Genet. Epidemiol. 32, 361–369. doi: 10.1002/gepi.20310
Glaubitz, J. C., Casstevens, T. M., Lu, F., Harriman, J., Elshire, R. J., Sun, Q., et al. (2014). TASSEL-GBS: a high capacity genotyping by sequencing analysis pipeline. PLoS ONE 9:e90346. doi: 10.1371/journal.pone.0090346
Iwata, H., and Ukai, Y. (2002). SHAPE: a computer program package for quantitative evaluation of biological shapes based on elliptic Fourier descriptors. J. Hered. 93, 384–385. doi: 10.1093/jhered/93.5.384
Jöst, M., Hensel, G., Kappel, C., Druka, A., Sicard, A., Hohmann, U., et al. (2016). The INDETERMINATE DOMAIN protein BROAD LEAF1 limits barley leaf width by restricting lateral proliferation. Curr. Biol. 26, 903–909. doi: 10.1016/j.cub.2016.01.047
Jung, S., Ficklin, S. P., Lee, T., Cheng, C. H., Blenda, A., Zheng, P., et al. (2014). The Genome Database for Rosaceae (GDR): year 10 update. Nucleic Acids Res. 42, D1237–D1244. doi: 10.1093/nar/gkt1012
Khan, M. A., Olsen, K. M., Sovero, V., Kushad, M. M., and Korban, S. S. (2014). Fruit quality traits have played critical roles in domestication of the apple. Plant Genome 7, 1–18. doi: 10.3835/plantgenome2014.04.0018
Kimura, S., Koenig, D., Kang, J., Yoong, F. Y., and Sinha, N. (2008). Natural variation in leaf morphology results from mutation of a novel KNOX gene. Curr. Biol. 18, 672–677. doi: 10.1016/j.cub.2008.04.008
Klein, L. L., Caito, M., Chapnick, C., Kitchen, C., O'hanlon, R., Chitwood, D. H., et al. (2017). Digital morphometrics of two North American grapevines (Vitis: Vitaceae) quantifies leaf variation between species, within species, and among individuals. Front. Plant Sci. 8:373. doi: 10.3389/fpls.2017.00373
Langlade, N. B., Feng, X., Dransfield, T., Copsey, L., Hanna, A. I., Thébaud, C., et al. (2005). Evolution through genetically controlled allometry space. Proc. Natl. Acad. Sci. U.S.A. 102, 10221–10226. doi: 10.1073/pnas.0504210102
Li, M., Frank, M. H., Coneva, V., Mio, W., Topp, C. N., and Chitwood, D. H. (2017c). Persistent homology: a tool to universally measure plant morphologies across organs and scales. bioRxiv. 1–25. doi: 10.1101/104141
Liebhard, R., Kellerhals, M., Pfammatter, W., Jertmini, M., and Gessler, C. (2003). Mapping quantitative physiological traits in apple (Malus × domestica Borkh.). Plant Mol. Biol. 53, 511–526. doi: 10.1023/A:1024886500979
Maeda, H., Akagi, T., and Tao, R. (2018). Quantitative characterization of fruit shape and its differentiation pattern in diverse persimmon (Diospyros kaki) cultivars. Sci. Hortic. 228, 41–48. doi: 10.1016/j.scienta.2017.10.006
Migicovsky, Z., Gardner, K. M., Money, D., Sawler, J., Bloom, J. S., Moffett, P., et al. (2016). Genome to phenome mapping in apple using historical data. Plant Genome 9, 1–15. doi: 10.3835/plantgenome2015.11.0113
Mohammadi, M., Tiede, T., and Smith, K. P. (2015). PopVar: a genome-wide procedure for predicting genetic variance and correlated response in biparental breeding populations. Crop Sci. 55, 2068–2077. doi: 10.2135/cropsci2015.01.0030
Purcell, S. (2009). PLINK v.1.07. Available online at: http://zzz.bwh.harvard.edu/plink/
Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A., Bender, D., et al. (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795
Talwara, S., Grout, B. W. W., and Toldam-Andersen, T. B. (2013). Modification of leaf morphology and anatomy as a consequence of columnar architecture in domestic apple (Malus × domestica Borkh.) trees. Sci. Hortic. 164, 310–315. doi: 10.1016/j.scienta.2013.08.025
Tian, F., Bradbury, P. J., Brown, P. J., Hung, H., Sun, Q., Flint-Garcia, S., et al. (2011). Genome-wide association study of leaf architecture in the maize nested association mapping population. Nat. Genet. 43, 159–162. doi: 10.1038/ng.746
Tsuge, T., Tsukaya, H., and Uchimiya, H. (1996). Two independent and polarized processes of cell elongation regulate leaf blade expansion in Arabidopsis thaliana (L.) Heynh. Development 122, 1589–1600.
Velasco, R., Zharkikh, A., Affourtit, J., Dhingra, A., Cestaro, A., Kalyanaraman, A., et al. (2010). The genome of the domesticated apple (Malus × domestica Borkh.). Nat. Genet. 42, 833–839. doi: 10.1038/ng.654
Visser, T., and Verhaegh, J. (1978). Inheritance and selection of some fruit characters of apple. II. The relation between leaf and fruit pH as a basis for preselection. Euphytica 27, 761–765. doi: 10.1007/BF00023712
Warton, D. I., Duursma, R. A., Falster, D. S., and Taskinen, S. (2012). smatr 3- an R package for estimation and inference about allometric lines. Methods Ecol. Evol. 3, 257–259. doi: 10.1111/j.2041-210X.2011.00153.x
Wünsche, J. N., Palmer, J. W., and Greer, D. H. (2000). Effects of crop load on fruiting and gas-exchange characteristics of “Braeburn”/M. 26 apple trees at full canopy. J. Am. Soc. Hortic. Sci. 125, 93–99.
Keywords: apple, leaf shape, morphometrics, elliptical Fourier descriptors, persistent homology, Malus domestica
Citation: Migicovsky Z, Li M, Chitwood DH and Myles S (2018) Morphometrics Reveals Complex and Heritable Apple Leaf Shapes. Front. Plant Sci. 8:2185. doi: 10.3389/fpls.2017.02185
Received: 06 October 2017; Accepted: 12 December 2017;
Published: 04 January 2018.
Edited by:Catherine Anne Kidner, University of Edinburgh, United Kingdom
Reviewed by:Fabrizio Costa, Fondazione Edmund Mach, Italy
Marie Monniaux, UMR5667 Reproduction et Developpement des Plantes (RDP), France
Copyright © 2018 Migicovsky, Li, Chitwood and Myles. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zoë Migicovsky, email@example.com