A non-destructive gel-based optical tomography imaging platform captured 3D RSA at day 9 after germination ( Figure 1A ). We used the RSA-GiA3D software to measure univariate root traits such as surface area, total root length and solidity (Galkovskyi et al., 2012; Topp et al., 2013), Dynamic Roots to measure local individual traits such as first order lateral root number and lateral root soil angle (Symonova et al., 2015; Jiang et al., 2019), and a TDA approach with PH to quantify the topology of the root structure (Li et al., 2017; Li et al., 2019; Figure 1B ). Because some traits are mathematically and/or phenotypically correlated, we employed the statistical approach variance inflation factor (VIF, Stine, 1995) to remove the redundant traits. We chose VIF because it has been effectively used in high-throughput plant phenotyping studies for this purpose for at least a decade (Chen et al., 2014; Falk et al., 2020), although advances in the statistics of high-dimensional data analysis have since proposed more sophisticated methods that can increase understanding of the most important response variables (Cheng et al., 2022) (see Conclusion). The remaining traits after VIF are highlighted in Supplementary Table 1 from the full list of all phenotypic traits. Then we performed principal component analysis (PCA) on the remaining univariate traits and PH traits and treated PCs as the multivariate traits ( Figure 1C ). Like our rationale for using VIF, we chose to use PCA because of its common use in plant phenotyping studies (Duc et al., 2023), but subsequent variations such as sparse PCA may increase interpretability of the multivariate traits (Zou et al., 2006, see Conclusion). We made GWA using univariate traits, PH traits, and multivariate traits by a Multi-Locus Mixed Model (MLMM), evaluating both the optimal and maximum models ( Figure 1D ; Ziegler et al., 2018). We then compared the results by co-aligning across the maize genome and evaluating correspondences in different window sizes.

Root systems are commonly measured by some intuitive topological descriptors such as number of tips and geometric descriptors such as root lengths and root-soil angles, which are useful but do not capture the entirety of 3D topological structure. Recently, TDA has demonstrated its wide application and success to extract complementary and comprehensive plant phenotypic traits such as 2D leaf shape, 2D root architecture and 3D inflorescence architecture (Li et al., 2017; Li et al., 2018a; Li et al., 2019). We used the TDA method PH to quantify the 3D RSA more comprehensively. We first extracted the surface voxels from the 3D model and manually cleaned the topological noise (e.g. branching touches) as much as we could. Then we assigned each surface voxel a value showing the curved distance from the voxel along the root to its top at the surface of the gel, which is called geodesic distance. Geodesic distance provides us the near true length of each root which is biologically important and known to be useful to distinguish different root architectures. The root examples in the Figures 2A, C, E show the colormap of the geodesic distance where red indicates high value, blue means low value. Starting from the voxel with the largest geodesic distance value, we plotted a bar to record the connected component information. Every isolated piece/blob is treated as one connected component. Then we kept adding new voxels as we continuously decrease the geodesic distance level. If the added voxels are connected to one of the components, the bar corresponding to that component will elongate; if the added voxels are the beginning of a new component (tip of a root), we start to plot to a new bar with the birth at the corresponding geodesic distance level; and if the added voxels merge two components, the shorter bar will die and the longer bar will elongate. At the end, the distance level is decreased to 0 and all the bars die except one since the root system is a single connected component. The bar graph forms a persistence barcode to record the topological information of the root structure ( Figures 2A, C, E ).

The persistence barcode associated with geodesic distance can record the topological structure and also some geometric information. For example, the number of bars is equal to the number of roots. The birth and death of each bar shows the depth of the root tip and emergence. The length of each bar indicates the length of the branch. However, it does not measure the angle. Thus, we combined it with a second persistence barcode that included a depth function to capture some of the angle information ( Figures 2B, D, F ). Because the depth of the root tip is highly related to the angle of the root relative to the soil surface, horizontal roots will have much shorter bars in the depth barcode ( Figure 2F ). For convenience, we simply name the two different persistence barcodes as geo-barcode and depth-barcode.

Barcodes need to be quantitatively compared to measure the similarity, quantify the variation and other statistics. Bottleneck distance is a robust metric to measure the distance between any two persistence barcodes (Cohen-Steiner et al., 2007). Intuitively, it measures how much minimum energy that it needs to take for reassembling a barcode to be the same as the other barcode. For each root, we computed the geo-barcode and depth-barcode. We measured the similarity between two roots by the square root of summation of bottleneck distance squared between two geo-barcodes and bottleneck distance squared between two depth-barcodes:

This measurement itself is a distance due to Cauchy-Schwarz inequality. Then we calculated the similarity matrix for the entire population and performed a multidimensional scaling method to get MDS scores using MATLAB function cmdscale(). For a non-Euclidean distance matrix, MDS will map the data into Euclidean space, while preserving the pairwise distance as well as possible. Then we performed PCA on the MDS scores which returned the same scores, but also provided some additional information such as the percentage variance per PCs. Thus we named it as PHGeodesicDepth_MDSPCs (GH_PCs). These GH_PCs are treated as PH traits. Similarly, G_PCs are the PH traits when only using geo-barcode to compute the similarity matrix. Please see step 5 of persistent homology in the Method section.

We introduced persistence barcode because it is more intuitive for biologists to understand. But the following feature is more straightforward to be calculated from another descriptor, persistence diagram. Persistence diagram is an equivalent descriptor to the persistence barcode ( Supplementary Figure 1A, B ). It is a 2D scatter plot: along the x-axis is the birth value for each bar and along the y-axis is the death value, unioning with the diagonal line. We turned this scatter plot into a Gaussian density estimator ( Supplementary Figure 1C , Adams et al., 2017) and performed PCA on the vectorized Gaussian density estimator. The PC scores are treated as another set of PH traits, named as PHDiagramKDE_PCs (PDD_PCs). Table 1 lists all the PH traits and their descriptions. In total, we have three groups of PH traits: G_PCs, GH_PCs, and PDD_PCs. G_PCs are the MDS-PC scores derived from the bottleneck distance of geo-barcode, GH_PCs are the MDS-PC scores derived from the combined distance of both geo-barcode and depth-barcode, PDD_PCs are the PC scores of the vectorized density estimator from geodesic persistence diagram.

Some traits could be strongly correlated. To reduce collinearity among explanatory variables ( Supplementary Figure 2A ), we applied a VIF. Simply, VIF_j = 1/(1-R_j ²) where the VIF for variable j is the reciprocal of the inverse of R² from the regression. To pick a proper threshold (see Materials and Methods), we calculated a sequence of thresholds and recorded the number of remaining traits and treated the thresholds as x values and remaining trait numbers as y values. We observed this data can fit well with a logarithm function (y=a*ln(x-b)+c, where a=15.71, b=0.8853, c=57.88, Supplementary Figure 1B ). We found the threshold at slope = 1 which means the VIF threshold increases at the same rate with the increasing of the traits number. We treated it as a high threshold. We picked a median threshold which is half of this threshold (value = 8; Supplementary Figures 2B, C ). Note that the first 29 G_PCs were removed because they are highly correlated to GH_PCs. The VIF chose to keep GH_PCs rather than G_PCs which may imply the depth function contributed to capture features (e.g. angle difference) that makes GH_PCs to be less correlated to other features. G_PC30 to G_PC36 were kept which may imply the depth function and geodesic distance captured some detailed differences which have much smaller variance compared to major topological differences. Then we performed PCA on the remaining traits including commonly used geometric univariate traits and PH traits. The output PC scores are recorded as mPCs, multivariate traits. The PH traits as well as multivariate traits are comprehensive but less intuitive compared to commonly used univariate traits. One way to interpret what a trait measures is to look at the genotypes with extreme values for that trait, but where other traits that are close to their mean values ( Figure 3 ). We also plotted out the distribution for each of the VIF remaining traits and first five mPCs to show the normed spread of the data ( Supplementary Figure 3 ).

Previous studies reported that the Illinois High Protein inbred (IHP1; derived from the IHP cycle 90 population) exhibits elevated N uptake and assimilation relative to the Illinois Low Protein inbred (ILP1; derived from the ILP cycle 90 population) (Moose et al., 2004; Lucas, 2013). By observing the RSA of both lines in gel, we found that IHP1 seedlings have nearly twice as many lateral roots than ILP1 with a high statistical significance (p=0.0002; Supplementary Figure 4 ). This finding is congruent with root adaptive responses to nitrogen availability (Drew et al., 1973; Robinson et al., 1999) and suggests that RSA may indeed have been inadvertently changed during recurrent selection for seed protein.

To further explore this possibility and map the genetic basis of root phenotypes using our expanded phenome, we performed a Genome Wide Association Study (GWAS). GBS (genotyping by sequencing; Elshire et al., 2011) was used to generate, 60,418 SNPs (single nucleotide polymorphisms). Using MLMM (Multi-Locus Mixed Model; Segura et al., 2012), one stringent model was chosen as the final/optimal with two different ways for multiple comparison, Bonferroni correction and E-BIC (Segura et al., 2012). The optimal model was picked after evaluating a few models with both forward and backward stepwise regressions. However, a few studies have shown that the Bonferroni correction increases the probability of producing false negatives (Segura et al., 2012; Ziegler et al., 2018). Therefore, we also employed a modified maximum model with only the forward stepwise regression and present both results ( Figures 4 , 5 ; Supplementary Figure 5 , Supplementary Tables 2 , 3 ).

For the optimal model, using a 1Mb window size (Yang et al., 2015b; Hu et al., 2017) as the QTL co-localization boundary, 102 total loci were detected for PH and all univariate traits in this study. Among them, only one 1Mb chromosome region is shared by both, representing <1% of QTL ( Figure 4A ). This locus (SNP29940) is located at ~170Mb on maize chromosome 4 (B73_Ref_v2) corresponding to PH traits, GH_PC5 and WidthDepthRatio3D ( Figure 5 ). The SNP hits the second exon of the ethylene receptor ETR2-like gene (GRMZM2G075368), also known as the ZmETR3, which was previously shown to be involved in root growth by regulation of ABA and/or auxin accumulation in root tips (Yang et al., 2015a; Li et al., 2018b). In addition, the same region was also found involved in regulation of another univariate median lateral root soil angle and the multivariate mPC3. although not the same SNP ( Figure 5 , SN29980, SNP29987). As described above, PCA was performed to generate multivariate traits ( Supplementary Figure 2 ). 62 QTL were identified for those multivariate traits with optimal model in MLMM ( Figure 4A ). 10 of these co-localized with the 102 loci identified in with the PH and univariate traits in the 1Mb window, representing 6.5% of total QTL. The maximum model identified more than twice as many QTL as did the optimal model, but only slightly increased the percentage of QTL shared by univariate and PH (0.98% to 3.3%) or multivariate and univariate + PH (6.5% to 7.3%; Figure 4B ). The list of trait-associated SNPs (TAS) for each phenotype class, model, and 1Mb overlaps thereof are provided in Supplementary Tables 2 – 4 . We also performed a similar analysis with no window size and got the majority of non-overlapped QTL ( Supplementary Figure 5 ). These results provide evidence of an expanded G to P map using both a TDA and a multivariate statistical approach, which are complementary and not exclusive.

An allelic effect size is how much of the total variation for a given phenotype the allele explains in the statistical model, and its direction (positive is larger and negative is smaller numerical value). Prior work showed that including multivariate and TDA of phenotypes could not only identify new loci ( Figure 4 ), but also loci of large effect size (Topp et al., 2013; Li et al., 2018a). To compute the effect size and direction in for each trait, we averaged the major allele trait values first, then we divided the estimate values of each QTL by the average values [effect size =estimate/mean(major allele)]. In this way, the positive value indicates increases on the major allele while negative values indicate increases on the minor allele phenotype. We found that most univariate traits have a narrow range of allele effect size (-0.25 to +0.25) while the effect size of both PH and multivariate traits ranged widely (-0.91 to +0.90) ( Figure 6 ; Supplementary Figure 6 ). These results reinforce the enhanced ability of data-driven phenotypes to identify large-effect QTL.

The maize IPSRI (Illinois Protein Strain Recombinant Inbreds) mapping population was initially described in Lucas et al. (2013). Initially, five plants from cycle 70 of Illinois Low Protein were crossed to five plants from cycle 70 of Illinois High Protein. After seven generations of random mating 200 individuals to reduce linkage disequilibrium (Dudley et al., 2007), 500 randomly selected individuals were self-pollinated for six generations. The 500 resulting recombinant inbred lines were genotyped with a set of 500 SNPs as described in Dudley et al. (2007). Clustering analysis using Euclidean distance of variation in marker genotypes produced 138 distinct groups where one IPSRI line was selected from each group to form a representative core set for phenotyping. These 138 core IPSRI inbreds were used in the root phenotyping analysis presented here. On average each genotype was replicated 3 times ( Supplementary Figure 7 ).

Seed preparation and growth conditions followed Jiang et al. (2019): “the growth medium was made with a modified 1/2x Hoagland solution pH 6.0 and solidified with gellan gum. The seeds were sterilized with 35% hydrogen peroxide for 20 minutes and rinsed four times with RO (reverse osmosis) water. After imbibing in RO water for 8 hours at 29°C in the dark, the seeds were sterilized again with 35% hydrogen peroxide for 10 minutes and rinsed four times with sterile water. The seeds were germinated at 29°C in the dark until the radicle reached 1-2 cm in length, approximately 48 hours. One seedling was planted into a glass growth cylinder sealed with Saran wrap - this constitutes a biological replicate. The cylinders were placed on a dark shelf at ambient conditions overnight for acclimation before moving them into a growth chamber starting at 4 DAG. The plants were lit with 315W Philips Ceramic Metal Halide bulbs, with a light intensity at the top of each jar of 700 µmol/m²/s. Humidity in the chamber was maintained at 50%, although the jars were sealed with Saran wrap. Temperatures were set to 28°C during the day and 24°C at night, with a 16/8h day/night cycle”.

As described in Griffiths et al. (2023), but with a different lens: the imaging setup consisted of an Allied Vision Manta G-609 machine vision camera (Allied Vision Technologies GmbH, Stadtroda, Germany) with a Kowa LM50SC 50mm 1” f/2.0 lens (Kowa, Japan) and an electronic turntable. The turntable operated in a water-filled tank to correct for light diffraction when imaging the glass cylinders. The glass cylinders were partly submerged to when placed in the center of the turntable. An LED flat panel light was used as a backlight to produce near binary images of the roots with a black silhouette of roots in the foreground against a white background. Root imaging took approximately 2 min per plant with 72 images collected over a 360-degree rotation. We studied 520 3D models and used the RSA-GiA3D software to measure global root traits (Galkovskyi et al., 2012; Topp et al., 2013) and the Dynamic Roots software to measure local individual traits (Symonova et al., 2015; Jiang et al., 2019). For the TDA analysis, we first extracted the surface voxel from the 3D model and saved them as.ply files. Then we manually cleaned some of the topological noises such as branching touches and loops as much as possible in Meshlab (Cignoni et al., 2008).

The pipeline for converting 3D root imaging data into persistent homology traits could involve five main steps: Extracting network (e.g. vertices and edges connecting vertices) from images, Assigning function values onto the network, Computing persistent barcodes with a filtration, Computing pairwise bottleneck distance for entire population, Performing statistical approach such as multidimensional scaling. The pipeline was conducted in MATLAB R2017a (The MathWorks Inc, 2017). A similar approach and code can also be found in Li et al. (2019).

The R function rcorr() in Hmisc package (Harrell, 2017) was used to compute the significance levels for Pearson correlations analysis. Both the correlation coefficients and the p-value of the correlation for all possible pairs of columns in the data table are returned. To create a graphical display of a correlation matrix and to highlight the most correlated variables (p< 0.05), the function corrplot() in the package of the same name (Wei and Simko, 2017) has been used in the study. Positive correlations are displayed in blue and negative correlations in red color. Color intensity and the size of the circle are proportional to the correlation coefficients ( Supplementary Figure 2A ).

To identify collinearity among explanatory variables, variance inflation factor (VIF) was used. Simply, VIF_j = 1/(1-R_j ²) where the VIF for variable j is the reciprocal of the inverse of R² from the regression. For each variable, one VIF value is calculated and variables with high values are removed. Instead of picking an arbitrary values in the range of 5-10 that are commonly used, we determined the definition of ‘high’ VIF is the x value at the slope is 1 for fitting a log function (y=a*ln(x-b)+c, where a=15.71, b=0.8853, c=57.88). Thus, the ‘high’ VIF=16.5953. We picked half of the ‘high’ VIF as the ‘median’ VIF (=8) in our study ( Supplementary Figure 2B ). After VIF, 33 univariant traits and 54 PH traits were left for further analysis ( Supplementary Figure 2C ).

The base R function prcomp() was used to perform the PCA with ‘scale =TRUE’ in the study. The number of PCs that captured > 90% variances totally and their rotated data (the scaled data multiplied by the rotation matrix) returned with ‘retx = TRUE’ were retained for subsequent analysis.

Broad sense heritability ( Supplementary Figure 8 ) was estimated as follows, H² = Vg/(Vg+Vr/nrep), which is the proportion of genetic variance out of total phenotypic variance. Vg indicates genetic variance, Vr represents the residual variance, while nrep is the mean number of repetitions for each genotype in the experiment. To generate the variance components, the R function lmer() from package lme4 (Bates et al., 2015) was used.

The multi-locus mixed-model (MLMM) approach was used to perform the genome wide association study. We chose MLMM because root system architecture traits are highly polygenic and MLMM has been demonstrated to identify small effect loci in structured populations that may also have large effect loci, although there is some risk of false-positives (Segura et al., 2012; Baseggio et al., 2021). Both forward and backward stepwise linear mixed-model regressions were used in the model, where the genetic variance and residual variance are estimated before each step. They then are used as follows to obtain generalized least-square (GLS) effect size estimates and F-test P values for each SNP: the SNP with the most significant association is then added to the model as a cofactor for the next step, and the P values for all cofactors are re-estimated together with the variance components. We used both the extended BIC (Bayesian Information Criterion) and the multiple Bonferroni criterion (mBonf) (alpha = 1.0 x e-5) for model selection. The mBonf was used to pick the optimal model. In addition, the maximum model which only includes the forward stepwise regression was also performed with 20 steps to ensure all potential trait-associated SNPs (TAS) were captured. To help determine the colocalization of different TASs, 1Mb window size has been applied which has been commonly used in some other previous studies (Yang et al., 2015b; Hu et al., 2017).

Effect estimates of significant TASs from the optimal model in MLMM were used to calculate the allele effect size. The fractions of the effect estimates and their corresponding average values of the major allele trait were used here to calculate the allele effect size. In this way, the positive value represents the QTL that have increases on the major allele, while negative values indicate QTL that have increases on minor allele.