3D reconstruction identifies loci linked to variation in angle of individual sorghum leaves

Image acquisition and annotation

A set of 366 individual plants from 236 genotypes of the sorghum association panel (Casa et al., 2008) were grown and imaged using the University of Nebraska-Lincoln’s automated phenotyping facility (Ge et al., 2016). Growth conditions and details of the image acquisition were previously described in (Gaillard et al., 2020b). RGB images were taken on April 11th (47 days after planting), 13th (49 days after planting) and 16th (52 days after planting) of 2018, using a Basler pia2400-17gc camera with a c6z1218m3-5 Pentax TV zoom lens. These images were taken at five different side view angles: 0°, 72°, 144°, 216°, and 288°side views and with a single additional image taken from directly overhead. Each image had a resolution of 2,454 × 2,056 pixels. The total number of side view images was dictated by the overall speed of the imaging system as the slowest imaging step dictated the time required to image all plants in the greenhouse. At the time, five side views was the maximum number which could be acquired by the RGB camera in the time required to capture one hyperspectral image which was the rate limiting step on data acquisition (Ge et al., 2016). Given this constraint, the angles 0°, 72°, 144°, 216°, and 288°were selected to capture five equidistant perspectives around a full 360°range of rotation. Images were calibrated as previously described (Gaillard et al., 2020b). Briefly, the calibration we applied corrects for the fact that pots are not perfectly centered on the axis of rotation of the turntable used to capture photos of the same plant from multiple angles and the axis of rotation for the turntable does not correspond to the optical center of the camera.

Two sorghum lines representing extremely different leaf architectures were selected and grown for paired manual and automated phenotyping: Btx623 (erect leaf architecture) and AS 4601 Pawaga (non-erect leaf architecure). Manual measurements of leaf angles were taken with a protractor on January 13th of 2021. Plants were imaged at three time points January 12th (70 days after planting), 13th and 14th 2021. These images were employed for 3D reconstruction and automated leaf angle measurements. Between 2018 and 2021 the physical RGB camera and lens employed by the Nebraska automated phenotyping greenhouse was replaced increasing the resolution and changing the focal length.

Plant voxelization and 3D skeletonization

A 3D volumetric grid was reconstructed for each plant on each day using the five side view images from different angles and single top view image described above and an extended voxel carving algorithm (Gaillard et al., 2020b). The algorithm employed segmented images from each view to carve a voxel grid with a resolution of 512³ voxels. The voxels remaining after carving were used to generate a 3D skeleton of the plant (Gaillard et al., 2020a).

Briefly, the generic voxel thinning algorithm implemented in the DGTal library was first run on the voxel set $\mathcal{V}$. This algorithm iteratively removes voxels from the reconstructed plant $\mathcal{V}$ until only its curve skeleton ${\mathcal{S}}_{raw}$ remains. This curve skeleton is a collection of curves whose shape faithfully represents the reconstructed plant. However, the curve skeleton ${\mathcal{S}}_{raw}$ includes some small spurious branches that do not correspond to separate organs of the plant being reconstructed. A machine learning classifier was run on the curve skeleton ${\mathcal{S}}_{raw}$ to distinguish between spurious branches and branches that correspond to real leaves with the spurious branches being discarded (Gaillard et al., 2020a). It is possible that the voxel grid contains more than one connected component i.e., the reconstructed plant is disconnected. In this case, the skeletonization algorithm will automatically connect branches together to form the plant skeleton. We denoted the final plant skeleton $\mathcal{S}$. Finally, the skeleton $\mathcal{S}$ is post-processed to identify the stem $\mathcal{T}$ and each individual leaf $\mathcal{L}=\left\{L_1,L_2,\dots ,L_n\right\}$. To differentiate voxels that are part of the stem or individual leaves, we consider that 1) any voxel that is shared by at least two branches is part of the stem, and 2) any voxel that is part of only one branch belong to the corresponding leaf (Gaillard et al., 2020a). Leaves were numbered based on the height of each leaf’s leaf-stem junction with the first leave being the leaf which joins the stem at the lowest point and the last leaf being the detected leaf which joins the stem at the highest point.

Leaf angle measurement from 3D skeletons

The angle of individual leaves is measured in a three step process. In the first step, the principal direction $\overrightarrow{v_{1,\mathcal{T}}}$ and the second direction $\overrightarrow{v_{2,\mathcal{T}}}$ of the stem are identified and measured by analyzing the 3D coordinates of voxels assigned to the stem of the plant skeleton $\mathcal{T}$ using PCA (Principal Component Analysis). In the second step, the principal direction $\overrightarrow{v_{L_n}}$ of the first 20 voxels of each leaf (the 20 voxels after the leaf-stem junction) L_n are computed, also through the use of the principal component analysis. In the final step, the principal directions of the stem and leaf are compared to compute two angles: (1) The polar angle θ ∈ [0°, 180°], which is the angle in the vertical plane, and (2) the azimutal angle ϕ ∈ [0°, 360°], which is the angle in the plane formed by $\overrightarrow{v_{1,\mathcal{T}}}$ and $\overrightarrow{v_{2,\mathcal{T}}}$. The angle θ is computed as the angle between the two 3-D vectors $\overrightarrow{v_{L_n}}$ and $\overrightarrow{v_{1,\mathcal{T}}}$. The angle ϕ is computed as the angle between $\overrightarrow{v_{2,\mathcal{T}}}$ and the projection of the leaf direction in the horizontal plane $\overrightarrow{v_{L_n}}-\left(\overrightarrow{v_{L_n}}\cdot \overrightarrow{v_{1,\mathcal{T}}}\right)\overrightarrow{v_{1,\mathcal{T}}}$.

Leaf angle measurements from 2D images

For each plant at each time point, all captured viewing angles were segmented into plant and not-plant using the published PlantCV protocol (Fahlgren et al., 2015; Gehan et al., 2017). Briefly, images were converted to grayscale, blurred, cropped to remove pots and backgrounds. These gray scale images were then converted to binary segmented images using the PlantCV thresholding function. Plant width was calculated as the distance between the left-most and right-most plant pixels in each image, and for each plant on each day the viewing angle with the greatest plant width was selected as the one most likely to represent a plane of phyllotaxy roughly perpendicular to the viewing angle of the camera. This widest image was skeletonized, the skeleton divided into leaf segments and the insertion angles of each leaf segment were calculated following the “morphology” protocol workflow from PlantCV’s documentation (Gehan et al., 2017).

Heritability and genome-wide association

Broad sense heritability was calculated using 48 replicated genotypes included in the study using the equation $H^2={\sigma }_G^2/\left({\sigma }_G^2+{\sigma }_e^2/2\right)$, where ${\sigma }_G^2$ is the total amount of variance that is due to genetics and ${\sigma }_e^2$ is the total residual variance. The variance components were extracted from the linear model: y_i = μ + t_i + e_i, where y_i is the mean yield of the genotype, µis the overall mean, t_i is the effect of genotype i and e_i is the residual error of genotype i. This was implemented in the software package lme4 (v1.1-23)(Bates et al., 2015).

Best linear unbiased predictors (BLUPs)(Robinson et al., 1991) were calculated for each set of leaf angle measurements used for GWAS with leaf angle treated as a random effect using the same linear model: y_i = μ + t_i + e_i. Any genotype with a BLUP more than 5 standard deviations away from the mean was treated as an outlier and removed from downstream analysis. Each outlier value was checked against image data to confirm the extreme values resulted from the plant itself and not reconstruction errors. In total one genotype was removed as an outlier for leaf 1, 1 was removed as an outlier for leaf 2, and one was removed from the aggregated leaf angle across leaves 1–4.

The genetic marker dataset employed in this study consisted of 569,306 SNP markers published as part of Miao et al. (2020a). The dataset was filtered to remove markers with a minor allele frequency of <0.05 and a heterozygousity of >0.05 among the 236 genotypes phenotyped in this study. This produced a set of 232,113 SNP markers employed for GWAS. GWAS was conducted using both mixed linear model (Yu, Pressoir & Briggs, 2006) and FarmCPU GWAS algorithms. MLM based GWAS was conducted using the GEMMA software package (v0.98.1) (Zhou & Stephens, 2012). FarmCPU based GWAS was conducted using the algorithm as implemented in the rMVP software package (v1.0.2) (Yin et al., 2020). In each model three principal components calculated from genetic marker data were incorporated as covariates. In addition, for the analysis conducted using GEMMA, a kinship matrix was fit as a random effect. For GEMMA, a statistical significance threshold of 6.39 × 10⁻⁷ equivalent to 0.05/78251 SNPs was employed. The number of SNPs used in calculating this statistical significance threshold was the effective number of independent tests represented in this genetic marker dataset calculated using GEC v.02 (Li et al., 2012). For FarmCPU, a subsampling strategy was employed to identify stable genetic associations. 100 random subsets of 87% of the data equivalent to 207/236 genotypes, were generated, and each was analyzed separately using the FarmCPU algorithm as described above. Individual SNPs which were assigned p-values more significant than 2.15 x 10⁻⁷ equivalent to 0.05/232,113 SNPs were treated as positive identifications. SNPs which were positively identified as linked to leaf angle in at least 10 of the 100 resampling analyses were classified as significant and stable associations.

Measurements from 3D reconstructions of sorghum plants

Initially two lines representing phenotypic extremes for leaf angle in the sorghum association population were selected to evaluate leaf angle measurement: BTx623, the reference genotype for sorghum which produces erect leaves (Fig. S1A), and AS 4601 Pawaga which produces non-erect leaves (Fig. S1B). Six plants were grown to late vegetative stage, the angles of the first four non-senesced leaves were manually measured, and the plants were imaged using a previously described automated greenhouse imaging system (Ge et al., 2016). Sets of 2D images collected from each plant from multiple viewing angles were processed to generate 3D voxel reconstructions using the method described in Gaillard et al. (2020a) (Fig. 1A). Three dimensional skeletonization and pruning was employed to identify individual leaves (Gaillard et al., 2020a) (Fig. 1B). And the angle of individual leaves relative to the principle direction of the stem was quantified from these segmented skeletonized leaves (Fig. 1C). The manual and Automated 3D measurements of median leaf angles across the first four leaves were highly positively correlated r = 0.98 with manual measurements taken with a digital caliper (Figs. S3A; S2). Manual and automated measurements for individual leaves exhibited a modestly lower but still high correlation r = 0.86 (Fig. S3B).

Having evaluated the performance of this 3D approach to leaf angle measurement on a small set of plants with ground truth manual measurement, the algorithm was next applied to a much larger experiment from which ground truth measurements had not been collected: An image dataset collected at multiple time points from a total of 366 sorghum plants representing 236 unique genotypes were used to create 1092 3D reconstructions, each representing a single plant on a single day. A total of 119 3D reconstructions failed. In 39 cases corresponding to 13 pots on each of the three time points no plant was present in any image, presumably the result of failed germination or the inclusion of calibration pots. In 35 cases calibration was unsuccessful. In the remaining 45 cases either the concordance between the original segmented 2D images and simulated 2D segmented images generated from 3D reconstructions was below 0.8 (Gaillard et al., 2020b), or the 3D skeletonization algorithm was not able to confidently identify the soil-stem junction (Gaillard et al., 2020a). Of the remaining 973 3D reconstructions, 2 were removed based on unexpectedly extreme leaf angle values. Manual evaluation of these two outlier reconstructions suggested that one resulted from a small plant with a large tiller (Fig. S4A). The reason for the poor reconstruction of the second outlier plant was not immediately apparent (Fig. S4B).

Sorghum plant reconstructions for the remaining 971 unique plant/timepoint combinations contained between 3 and 19 predicted leaves, with a median of 10 (Fig. S5). Median leaf angle declined, e.g., became more erect, starting from the bottom most leaf and proceeding upwards to leaf six (Fig. S6). Leaf angles measured from 3D skeletons for leaves above leaf six exhibited a much wider distribution of values, which may reflect the difficulty of accurately quantifying leaf angle for leaves still emerging from the leaf whorl where the ligule/auricle junction is not yet visible. It was not possible to automate the separation of mature and emerging leaves, so it was unclear how much accuracy was lost for uppermost leaves.

As this analysis was conducted using a previously collected sorghum image dataset for which ground truth leaf angles were not manually scored, it was not possible to directly evaluate how well this larger set of measured leaf angles represented the angles which would have been measured by hand from the exact same plants. However, because the genotypes used in this study were largely homozygous inbred lines, it was possible to grow plants with identical or near identical genomes to the same stage of development and collect ground truth measurements from these plants. Independent grow outs of the same set of genotypes do not produce perfectly equivalent phenotypes. Measurements of plant leaf angle collected in Iowa and Nebraska from the same population were compared in order to estimate how much correlation in leaf angle would be expected between two independently grown sets of plants at the same stage of development—reproductive stage in this case—collected in different environments. The correlation between these two manual datasets data was significant (r = 0.60; p = 0.0002; pearson). Leaf angle measurement were manually scored from a set of 79 lines from the sorghum association panel grown to a similar late vegetative stage of development under greenhouse conditions were significantly correlated with automated 3D measurements (r = 0.53; p = 0.0016; pearson) (Fig. S7). However, analysis of the same image data using conventional 2D approaches to quantifying leaf angle from single images produced values which were not significantly correlated with either manual measures of the same genotypes at the same stage, or leaf angles measured from 3D reconstructions using the same image data (Fig. S7). Manual measurements of plant leaf angle between Iowa and Nebraska datasets were anticipated to be higher than the automatic 3D versus manual measurements of the greenhouse plants where different plants of the same genotypes were compared. The leaf angles for each genotype in both the Iowa datasets and Nebraska dataset of mature plants were the means of up to 18 and 2 replicates of genotypes respectively. This would reduce environmental, genotype by environment and many other residual effects. The dataset composing of 3D reconstructed plants and that of the 79 lines under greenhouse conditions generally represented measurements of multiple leaves from only a single plant.

Assuming random error in quantification, as opposed to genetically controlled error (Liang et al., 2017), if the angles of upper leaves are measured with lower accuracy, the heritability—i.e., the proportion of total variance attributable to differences between genotypes—of leaf angle measurements for these leaves would likely be lower. A total of 48 sorghum genotypes were replicated between two and eight times each among the 366 unique sorghum plants imaged. Only the first seven leaves of each plant were considered in our analysis. Higher leaves were generally not mature as evident by the auricle not being extended from the sheath. This would limit the accuracy of leaf angle measurements of those after the seventh leaf. The heritability of 3D leaf angle measurements tended to be higher for the first five leaves than for leaves six or seven on each of the three time points measured (Fig. S8A). The median measurement for the angle of a particular leaf across three time points tended to be more heritable than measurements of the same leaf at individual time points, consistent with expectations for repeated measurements with independent error (Fig. S8A). Aggregating measurements across sequential leaves also tended to increase heritability relative to single leaf measurements (Fig. S8B). Of the aggregated leaf ranges within the aggregated time points,the median of leaves 1 to 4 had the highest heritability value of 0.72 (Fig. S8B). The heritability of measurements of leaf angle obtained from conventional 2D analysis of the same set of sorghum images were generally much lower than 3D measurements with no genetic effect detected in a number of cases (Fig. S9).

Genetic loci associated with leaf angle variation in the SAP

A published set of SNP markers for the sorghum association panel (Miao et al., 2020a) was employed to conduct a genome wide association study for both the 3D and 2D aggregated leaf angle data across the first four leaves and all three time points. Analysis with GEMMA (Zhou & Stephens, 2012), a mixed linear model based algorithm, did not identify any significant markers for the automated 2D angles (Fig. S10) but for 3D derived angles, identified a single large peak for leaf angle located on chromosome 7 between 59.8 and 59.9 MB. This position corresponds to the genomic location of dwarf3 (Sobic.007G163800), an ATP-binding cassette (ABC) transporter involved in auxin export (Multani et al., 2003) that has been previously associated with variation in leaf angle in sorghum RIL populations (Truong et al., 2015) and an association study (Mantilla-Perez et al., 2020). As dwarf3 is a large effect locus with functionally variable alleles segregating at high frequencies in this population, we speculated that the effects of additional loci on angle might be masked in the MLM-based GWAS results (Fig. 2A). The FarmCPU algorithm, a multi-locus GWAS algorithm that can control for large effect loci to increase the likelihood of identifying additional true positive signals, was employed to analyze the same trait and genotype datasets. The significance of FarmCPU results were controlled for using a resampling-based strategy (Valdar et al., 2009). A total of eight markers clustered at six locations in the genome showed significant associations (resample model inclusion probability (RMIP) > 0.1) with median angle of leaves 1–4 across all three time points in > 10% of resampled FarmCPU GWAS results, including the dwarf3 locus (Fig. 2B). None of the five non-dwarf3 clusters were associated with a set of 19 a priori candidate sorghum genes identified based on syntenic orthology to genes known to influence leaf angle in maize or rice. However, three of the six clusters associated with variation in leaf angle were consistent with QTL detected in a nested association mapping study for leaf erectness using a population of 2,200 sorghum RILs (Olatoye, Hu & Morris, 2020).

The reconstruction of full architectures of individual plants provided the opportunity to assess how much overlap was present in the loci showing significant associations with variation in leaf angle for individual leaves. Separate analyses were conducted for the angle of the first, second, third, and fourth identified leaves of each sorghum plant. It must be noted that these likely do not correspond to the leaves which would be called leaf 1, leaf 2, leaf 3, and 4 as early sorghum seedling leaves senesce and die as the plant continues to develop. However, the first identified leaf should always be an earlier developing leaf than the second and third and and so on. Eight total markers showed significant associations in >10% of sub-sampled GWAS results for the angle of the first identified leaf although three of these were grouped together in the pericentromeric region of chromosome 3 and may represent a single underlying signal, five clusters of ≥ 1 SNPs showed significant associations in >10% of sub-sampled GWAS results for the second identified leaf, five showed significant associations in >10% of sub-sampled GWAS results for the third identified leaf (Fig. 3), and no significant signals were identified for the fourth identified leaf. The signal associated with the dwarf3 locus was detectable for the second and third leaves but not the first leaf. One of the six significant clusters identified in GWAS for the angle of the first leaf was located on sorghum chromosome 2, 420 kilobases from the sorghum ortholog (Sobic.002G272400) of ZmTac1 (Ku et al., 2011). One of the five significant associations for the second leaf was located at 54.99MB, 503 kilobases from the sorghum ortholog (Sobic.006G190400) of the maize gene UPA2/ZmRAVL1 (Tian et al., 2019). In addition to the dwarf3 locus, a second shared signal was identified between GWAS for the angle of the second identified leaf and GWAS for the angle of the third identified leaf at position 52.75MB on sorghum chromosome 10. This hit did not correspond to either published loci associated with leaf erectness in the sorghum NAM population or sorghum orthologs of candidate leaf angle regulatory genes from maize or rice. Mixed linear model based, leaf by leaf GWAS results, showed 3 significant associations for leaf 2 and one for leaf 3 (Fig. S11). Two of the associations for leaf 2, one on chromosome 7 and the other on 10 was in close proximity to signifcant SNPs via the leaf by leaf FarmCPU resampling (5.031 and 29.997 kilobase differences respectively). The association for leaf 3, on chromosome 6, was also in close proximity to a significant SNP found in the FarmCPU resampling method (34.572 kilobase difference).