Dynamics of plant phenome can be accurately predicted from genetic markers

doi:10.21203/rs.3.rs-4958737/v1

Download PDF

Article

Dynamics of plant phenome can be accurately predicted from genetic markers

https://doi.org/10.21203/rs.3.rs-4958737/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Molecular and physiological changes across crop developmental stages shape the plant phenome and render its prediction from genetic markers challenging. Here we present dynamicGP, an efficient computational approach that combines genomic prediction with dynamic mode decomposition to characterize temporal changes in the crop phenotype and to predict genotype-specific dynamics for multiple traits. Using genetic markers and data from high-throughput phenotyping of a maize multi-parent advanced generation inter-cross population, we show that dynamicGP outperforms a state-of-the-art genomic prediction approach for multiple traits. We demonstrate that the developmental dynamics of traits whose heritability varies less over time can be predicted with higher accuracy. The approach paves the way for interrogating and integrating the dynamical interactions between genotype and phenotype over crop development to improve the prediction accuracy of agronomically relevant traits.

Biological sciences/Plant sciences/Plant development

Biological sciences/Biological techniques/High-throughput screening

Biological sciences/Genetics/Quantitative trait

The phenome of a plant comprises the entirety of traits expressed at any given time, and is the integrated outcome of the effects of genetic factors, environmental conditions and their complex interactions. Understanding how the crop phenome changes over time can help predict individual traits at specific time points in crop development. However, this problem is challenging not only because of the intricate dependence between individual traits, but also due to differences in how the phenomes of specific genotypes change over the plant life cycle.

The classical approach of genomic prediction (GP) in crops trains machine learning models using data on traits measured in a population of genotypes at a specific time point based on genetic markers¹. These models can be used to predict the traits in other genotypes for which genetic markers are available, thus foregoing measurements; the accuracy of predictions depends on representativeness of the training set and on the genetic architecture of the predicted trait². As a result, GP accelerates genetic gain and reduces the need for labor-intensive field phenotyping in a test set of genotypes. While GP has been extended to allow and improve the simultaneous prediction of multiple correlated traits, comprising the crop phenotype^3–5, existing GP approaches have not yet addressed the problem of predicting the dynamics of plant phenotypes, i.e., the expression of multiple traits at different time points across the entire period of growth of the plant.

Advances in high-throughput phenotyping (HTP) technologies, including non-destructive approaches such as: hyperspectral, multispectral, fluorescence, and thermal imaging, have enabled the detailed capture of plant traits at various growth stages and have resulted in large-scale time-resolved heterogeneous data sets on crop phenotypes^6–10. These developments offer the possibility of predicting time-resolved crop phenotypes from genetic markers, with the potential to revolutionize GP. Recent efforts have cast the prediction of the dynamics of crop phenotypes as a spatio-temporal sequence prediction problem. The proposed solutions rely on deep learning approaches, such as spatio-temporal long short-term memory (ST-LSTM) or Memory-in-memory (MIM) networks¹¹, generative adversarial networks¹² (GAN), and U-net encoder-decoder¹³. However, rather than directly predicting time-resolved crop traits, these approaches aim to predict a time series of images that are then used to infer the traits of interest.

The time-resolved prediction of crop phenotypes can also be addressed by dynamic mode decomposition (DMD)¹⁴. DMD is an advanced data-driven method used in decomposing data from complex, time-dependent systems into spatio-temporal modes that provide low-dimensional description of the systems’ dynamics. These modes render DMD particularly useful for understanding and predicting the behavior of systems characterized by high-dimensional data. DMD has been applied across different fields, including: fluid mechanics, engineering, epidemiology¹⁴, and neuroscience¹⁵; however, no attempts have yet been made to apply DMD to data from HTP technologies or to combine this method with genetic markers.

Here we combined DMD with GP to predict the genotype-specific dynamics of crop phenotypes using HTP data measured on a maize population of recombinant inbred lines (RILs). Our approach, termed dynamicGP, showed consistent prediction accuracy for multiple traits and across multiple time points. When compared with a baseline GP approach, dynamicGP exhibited higher accuracies at the majority of time points and irrespective of traits, showing that it can effectively capture the modes, and therefore the developmental dynamics of untested genotypes based on genetic marker data.

DMD accurately forecasts time-resolved phenotypes

DMD is a data-driven diagnostic and time-series prediction method that has found diverse applications from engineering to finance¹⁴ and has only recently been applied in the biological sciences¹⁶. To illustrate the effectiveness of DMD, we applied it to predict the dynamics of growth-related traits derived from multi-modal HTP imaging measured at 25 time points, five days per week for five weeks, beginning at day 15 after sowing, in a maize multi-parent advanced generation inter-cross (MAGIC) population of 347 RILs of which genotyping data were available for 330 (see Methods). A one-day shift of the second of three HTP experiments resulted in a two-day gap and thus five consecutive measurements per week, as these were adjusted between the three experiments based on the days after sowing (DAS). By using network clustering, we selected 50 traits as representatives for the compendium of 498 top, side, and combined measurements of leaf-colour and -geometric traits (see Methods, Tables 1 and S1).

Table 1

Representative Traits
	Q & R
Minimum	plant average hue value in side images
Q1	Top view RGB blue pixels
Mean	Side view RGB convex hull
Q3	Top view RGB average saturation
Maximum	Side view fluor standard deviation of lightness

For a single maize line, the time-resolved phenotype was arranged into a $\:p\times\:T$ matrix $\:\mathbf{X}$, where $\:p$ is the number of traits and $\:T$ is the number of time points (Fig. 1a). Two submatrices, $\:{\mathbf{X}}_{1}$ and $\:{\mathbf{X}}_{2}$, which are offset by a single time point are derived from $\:\mathbf{X}$ (Fig. 1a), and then used to calculate a best fit linear $\:p\times\:p$ operator (matrix) $\:\mathbf{A}$, linking the phenotype at one time point with the phenotype in the following time point (Fig. 1a). We note that the usage of linear operator does not assume that the traits change linearly with time. The operator $\:\mathbf{A}$ is usually non-symmetric and can be computed directly using the classical DMD approach (see Algorithm 1, Methods).

We first tested how accurately the operator $\:\mathbf{A}$, calculated using the data from a specific genotype, was able to recreate the underlying time-resolved phenotypes. These tests were performed in two scenarios, iterative and recursive. In the iterative scenario, the traits at each time point were predicted based on the measured data on the traits from the preceding time point. In the recursive scenario, the traits at each time point were predicted based on the traits predicted from the preceding time point, with measurements available only for the initial time point (see Methods).

We found that the operator $\:\mathbf{A}$ derived from the classical DMD approach was able to near perfectly recreate the training data when used iteratively, with a mean prediction accuracy of ~ 1 across all traits over the time points within a five-day block of daily measurements, and a mean prediction accuracy of 0.70 (± 0.20) for the time points immediately following the two-day gap, for which no data were available (Fig. S1). In the recursive scenario, the classical DMD approach resulted in a perfect recreation of the training data for the first five-day block; however, following the first two-day gap, the prediction accuracy decreased rapidly as a result of error propagation (Fig. S1). These results indicated that the classical DMD approach yielded models which overfitted to the training data, and were not robust to slight deviations in the input data.

We additionally tested another DMD algorithm, named Schur-based DMD, which has improved numerical stability¹⁷. This algorithm uses the singular vectors associated with the $\:r$ largest singular values of $\:{\mathbf{X}}_{1}$ alongside the Schur decomposition to reconstruct another operator $\:{\mathbf{A}}_{\mathbf{r}}$ (Algorithm 2, see Methods), which is used in place of the operator $\:\mathbf{A}$ in the prediction procedure outlined above. We note that like the operator $\:\mathbf{A}$, $\:{\mathbf{A}}_{\mathbf{r}}$ is not necessarily symmetric. In comparison to $\:\mathbf{A}$ from the classical DMD, we found that use of $\:{\mathbf{A}}_{\mathbf{r}}$ in the iterative scenario resulted in decreased mean prediction accuracy across all traits and all time points of 0.84 (± 0.18). However, in the recursive scenario, which is particularly relevant for practical applications, we observed robust performance, supported by the mean prediction accuracy of 0.78 (± 0.16) and 0.79 (± 0.13) across all traits and all time points and for the last time point, respectively (Fig. S2). These findings demonstrate that the Schur-based DMD exhibited good predictive performance, and prompted our investigation of the extent to which the elements of $\:{\mathbf{A}}_{\mathbf{r}}$ and corresponding matrices used in its derivation can be predicted from genetic markers.

Building blocks of dynamicGP are heritable and can be predicted from genetic markers

For each genotype, the operators $\:\mathbf{A}$ and $\:{\mathbf{A}}_{\mathbf{r}}$ as well as the intermediate component matrices in the Schur-based DMD (Algorithm 2, see Methods) can be determined from phenotypic data. Here, we treat each individual entry of these matrices as a trait to be predicted by GP models for the operators $\:\mathbf{A}$ or $\:{\mathbf{A}}_{\mathbf{r}}$ (Fig. 1b). Our approach, dynamicGP, uses genetic markers to predict matrix entries, which can then be combined with selected phenotypic measurements to make longitudinal predictions of plant development for unseen genotypes.

Here, the use of the Schur-based DMD has the advantage that its intermediate components contain fewer entries than either $\:\mathbf{A}$ or $\:{\mathbf{A}}_{\mathbf{r}}$. This renders the components of Schur-based DMD a potential entry point for the computationally efficient prediction of $\:{\mathbf{A}}_{\mathbf{r}}$ from genetic markers. We examined the marker-based heritability of the entries of all the intermediate matrices of the Schur-based DMD using 70,846 single nucleotide polymorphisms (SNPs) from the MAGIC maize population (see Methods). These matrices include those of the singular value decomposition of $\:{\mathbf{X}}_{1}=\mathbf{U}\varvec{\Sigma\:}{\mathbf{V}}^{\mathbf{T}}$, the rank reduced representation of the operator $\:\mathbf{A}$, denoted by $\:\stackrel{\sim}{\mathbf{A}}$, the matrices $\:\mathbf{Q}$ and $\:\mathbf{R}$ resulting from its Schur decomposition as well as the projected DMD modes, $\:\varvec{\Phi\:}$ (Algorithm 2, Methods).

We first observed that non-zero heritability measures for the entries of these matrices were obtained by using the first two singular vectors ($\:r=2$, see Methods). Specifically, we observed that the mean heritability of $\:\stackrel{\sim}{\mathbf{A}}$ entries was 0.28, ranging from 0.13 to 0.43. Similarly, the entries of matrix $\:{\mathbf{U}}_{\mathbf{r}}$ exhibited mean heritability of 0.35, ranging from ~ 0 to 0.75. Moreover, the entries of $\:{\varvec{\Sigma\:}}_{\mathbf{r}}$ displayed mean heritability of 0.47, while those of $\:{\mathbf{V}}_{\mathbf{r}}$ showed mean heritability of 0.20, ranging from 0.03 to 0.39 (Fig. 2a). Regarding the Schur decomposition of $\:\stackrel{\sim}{\mathbf{A}}$, the mean heritability of $\:\mathbf{R}$ entries was 0.30, ranging from 0.20 to 0.39. In addition, the $\:\mathbf{Q}$ entries demonstrated a mean heritability of 0.04, ranging from ~ 0 to 0.09. Notably, the mean heritability of $\:\varvec{\Phi\:}$ entries was 0.10, ranging from ~ 0 to 0.42 (Fig. 2a). These findings indicated that substantial parts of the building blocks of dynamicGP exhibit moderate to high heritability, suggesting that they may be predictable from genetic markers, and we can therefore use them to predict operator $\:{\mathbf{A}}_{\mathbf{r}}$.

To test this hypothesis, we used 20 iterations of 5-fold cross validation to examine the genomic predictability of the components of the singular value decomposition of $\:{\mathbf{X}}_{1}$ and the entries of $\:\stackrel{\sim}{\mathbf{A}}$ using SNPs. The entries of $\:\stackrel{\sim}{\mathbf{A}}$ exhibited a mean prediction accuracy of 0.24 (± 0.15), while the entries of matrix $\:{\mathbf{U}}_{\mathbf{r}}$ demonstrated a slightly higher mean accuracy of 0.31 (± 0.16). Moreover, the entries of $\:{\varvec{\Sigma\:}}_{\mathbf{r}}$ showed a mean prediction accuracy of 0.39 (± 0.12), whereas those of $\:{\mathbf{V}}_{\mathbf{r}}$ displayed a slightly lower mean accuracy of 0.19 (± 0.14). Regarding the Schur decomposition of matrix $\:\stackrel{\sim}{\mathbf{A}}$, the entries of $\:\mathbf{R}$ exhibited a mean accuracy of 0.27 (± 0.13), while those of $\:\mathbf{Q}$ showed a lower mean accuracy of 0.08 (± 0.12); in addition, the mean accuracy of the entries of $\:\varvec{\Phi\:}$ was 0.15 (± 0.14) (Fig. 2b). These results suggested that the building blocks of dynamicGP can indeed be predicted from genomic data using standard genomic prediction models.

Differences in performance between the two versions of dynamicGP

We next used the predicted values for each entry in matrices $\:\varvec{\Phi\:}$ and $\:\mathbf{R}$ to recreate the predicted operator $\:{\mathbf{A}}_{\mathbf{r}}$ for each line (Fig. 1c). These predicted operators were used to make longitudinal predictions of all 50 traits over the entire time domain using the two versions of dynamicGP, iterative and recursive, defined similarly to the usage of DMD (see Methods).

We found that the iterative approach yielded more consistent accuracies across the investigated time points compared to the recursive approach (Fig. 3a). Across all traits and time points, the iterative approach resulted in a mean prediction accuracy of 0.34 (± 0.30) in 20 iterations of 5-fold cross validation. The average blue value of plant pixel colors, which quantifies the blue-yellow color space of visible light (see Table S1 for trait descriptions), emerged as the best-performing trait, with a mean prediction accuracy of 0.79 (± 0.07). In contrast, the plant average hue value in side images, quantified in the hue, saturation and value color space, showed the lowest prediction performance of -0.38 (± 0.12). Conversely, the recursive approach had a mean prediction accuracy of 0.17 (± 0.19) across all traits and time points. Here, the standard deviation of the y values in the fluorescent xyz color space from top images could be predicted with mean prediction accuracy of 0.52 (± 0.11), while the plant average hue value in side images, as quantified in the hue, saturation and value color space, again displayed the lowest mean prediction accuracy of -0.30 (± 0.13). The recursive predictions, as expected, tended towards zero at later time points as errors compounded.

We observed that traits exhibiting consistent heritability over time, indicated by smaller values for the coefficient of variation, tended to demonstrate higher prediction accuracies across all time points. This relationship was supported by a moderate negative Pearson correlation coefficient of -0.46 between the mean prediction accuracy for the traits from iterative dynamicGP across all time points and the coefficient of variation of the heritability estimates across the entire time points (Fig. 3d); the corresponding correlation coefficient for recursive dynamicGP was − 0.37 (Fig. S10). Therefore, in line with expectations, traits whose heritability does not vary during development were found to be better predicted by dynamicGP approach.

Performance of dynamicGP compared to baseline models

To assess the predictive performance of dynamicGP, we compared it with a GP baseline approach using RR-BLUP models, which is the standard in breeding programs. The RR-BLUP models were trained with data on each trait from the first time point in a training set and used to predict the trait in all subsequent time points for the testing set. Inspecting the difference of performance between dynamicGP and the baseline for traits with different predictability, we found that the iterative version of dynamicGP outperformed the baseline RR-BLUP models consistently across all investigated time points, while the recursive version performed competitively, outperforming the RR-BLUP models for all time points between day 16 and day 33 (Fig. S8). Specifically, the mean prediction accuracy across all traits for the first time point, on which the models were trained, was 0.24 (± 0.15). While this accuracy was maintained when the baseline RR-BLUP models were applied to predict the traits over the first four subsequent time points, already at the second time point, both versions of dynamicGP outperformed the baseline (Fig. S8). The recursive version of dynamicGP was superior to the baseline for 12 time points, after which the accuracies of the two models were closely aligned (Fig. S8). Moreover, the iterative version of dynamicGP outperformed both the baseline and the recursive version of dynamicGP for every time point.

We examined the differences in mean prediction accuracy across all time points between dynamicGP and the baseline models for each trait. The best performing trait using the recursive method, the standard deviation of the y values in the fluorescent xyz color space from top-down images, showed mean prediction accuracies above 0.5 of all time points, which was 60% higher than the best performing trait in the baseline models (Fig. 4). In the iterative method the best performing traits, namely the average blue value of pixels as well as the average brightness measurements of green pixels and pixels with 27% saturation (Table S1), showed mean accuracies above 0.7 of all time points, which was ~ 2.4-fold higher than the best trait from the baseline model predictions (t-test, p-values < 0.001 for all time points, Fig. 4). Furthermore, both versions of dynamicGP yielded mean squared errors much lower than the RR-BLUP baselines (Fig. S9).

Developing computational approaches that can predict the dynamics of crop growth and developmental plasticity for unseen genotypes in new environments holds the promise to revolutionize breeding practices. Such approaches can provide insights not only into how the genotype affects the time-resolved phenome, but also can highlight the temporal dependence between traits comprising the phenome. While there have been significant developments in GP of multiple-traits^3–5, they essentially operate with snapshot data and do not capture the dynamics of crop traits.

Here, we introduced dynamicGP, a computational approach that resolves the problem of predicting the growth dynamics across development in crops for which time-series phenome measurements for multiple genotypes are available from HTP platforms^6–10. DynamicGP essentially applies GP to the building blocks of DMD, thus facilitating the prediction of multiple traits over time of unseen genotypes. Using a time-resolved HTP data collected for five days per week for five weeks beginning at day 15 after sowing from RILs of a maize MAGIC population, we showed that both the iterative and the recursive versions of dynamicGP outperformed the baseline GP approach.

In addition, we showed that at the predictability threshold of 0.5 for two of the traits, recursive dynamicGP exhibited better performance compared to the baseline GP approach. These traits were both top-down photographs measured in the fluorescence color space. For the iterative dynamicGP, we found three traits above a threshold of 0.7 where dynamicGP outperformed the baseline, all measured in the visible colour space, namely the average blue-value of the plant pixel colors of images from above, average brightness of fraction of green-colored pixels in top-down images, and average brightness of fraction of pixels with mean 27% saturation in top-down images (Table S1). Further, we observed that the mean squared error of dynamicGP was lower for both the iterative and recursive versions than the baseline, with log-transformed MSEs of -3.29 and − 3.57 for iterative and recursive, respectively and − 0.08 for the baseline (Fig. S9). Therefore, dynamicGP yields predictions which are closer to the true values than the RR-BLUP baseline models.

The performance improvement of dynamicGP over the classical GP approach is due to the usage of time-resolved phenotypic data. While availability of such data is rapidly increasing due to deployment of HTP in controlled environments and on the field, dynamicGP is more resource intensive than the classical GP. However, dyanamicGP does not rely on model fitting (e.g., of growth curves) and subsequent usage of the model parameters in GP, thus avoiding propagation of error and often unfounded assumptions about models that best describe the dynamics of single traits. Instead, dynamicGP builds on the advantages of DMD that allows the simultaneous investigation of multiple traits over time.

Although dynamicGP is applicable with any time-resolved phenotypic data on multiple traits, we tested its performance with traits retained after correlation analysis. The selection of traits was conducted to avoid the consideration of highly correlated traits that may affect predictive performance. Like with any data-driven approach, we hypothesize that the performance of dynamicGP depends on the number of traits included and time points considered, which may vary between data sets. Both the recursive and iterative versions of dynamicGP offer the possibility to make predictions about plant phenotypes of unseen genotypes given measurements of specified time points in time, thus reducing the need for extensive temporal measurements. In this regard, while dynamicGP provides the means to overcome the challenge of predicting temporal phenomes for unseen genotypes, it inherits the limitations of GP, such as the transferability of the models to different populations and environments².

Future developments of dynamicGP can rely on extensions of DMD to consider effects of environmental factors¹⁸. These will facilitate further refinements of the proposed approach that are expected to have very substantial impact on breeding crop varieties adapted to particular regions as well as precision agriculture.

Data generation

In three independent high-throughput phenotyping experiments, 347 RILs and 9 founder lines of a MAGIC maize population19 were screened for natural variation in the HTP facility for large plants at the Leibniz Institute of Plant Genetics and Crop Plant Research in Gatersleben, Germany²⁰. The panel was divided into three equal-sized subgroups and analysed in three consecutive experiments: 2113MH (116 RILs and 9 founders), 2121MH (116 RILs and 9 founders), and 2137MH (115 RILs and 8 founders). To ensure compatibility of the data, each HTP experiment was supplemented by an additional randomly drawn overlap of 21, 21 or 27 RILs, respectively, half of which came from each of the other two experiments. The experimental design followed a complete randomised block design with three replicates per genotype (three carriers with two plants each) overlapping with three blocks corresponding to four of the twelve lanes of the HTP system. Image acquisition began at 15 DAS for 31 time points per experiment and six days per week. An image acquisition event was a carrier with two plants phenotyped by one image in top view and four images in side view at 22°, 45°, 112° and 135°. Altogether 347,740 images were taken (2113MH, 114710 images, 91768 side view, 22942 top view; 2121MH, 116840 images, 93472 side view, 23368 top view; 2137MH, 116190 images, 92952 side view, 23238 top view). Traits were derived from images by IAP version 2.3.0²¹ and were analysed by alignment by days after sowing (DAS). A shift of the sowing date by one day in 2121MH allowed alignment to DAS for only 5 days per week, so that 2 days were missing each week. Consequently, the time points were not consecutive, but were organized in five-day blocks, with one measurement taking place each day and no measurements for two days over five weeks.

The analysed data comprised measurements for 498 traits from side, top and combined categories derived from multi-colour-space image analysis (Fig. 3e, Tables S1²² Heuermann et al. 2024). A linear mixed model with random effect of genotype, experiment, lane, genotype-by-experiment interaction and replicate nested in experiment was built for each time point. Trait heritability was calculated by the ratio between genotype variance component and total variance component. Summation of best linear unbiased prediction (BLUP) value of genotype and fixed effect of intercept were used as the final trait values.

Genotyping data was available for 330 Maize RILs consisting of 79,557 high-quality SNPs derived from SPET genotyping²³. Missing genotyping data were imputed using Beagle 5.2²⁴. We then removed SNPs with a minor allele frequency < 0.05, resulting in a final set of 70,846 SNPs for GP models.

The phenotypic data sets used in this study are available at: https://doi.ipk-gatersleben.de/DOI/0c9c6237-41f2-411f-a51e-809eb23d1088/f844533e-d775-46dd-8523-d485591f6ea8/2/1847940088

Clustering and trait selection

Traits were first clustered using the Mantel correlations over all genotypes and time points, represented as a matrix. To this end, a network was created with nodes representing traits and edges denoting inter-trait Mantel correlations of value above 0.96. Modularity clustering, as implemented in the R package “igraph”²⁵, was used to create clusters of traits, from which the trait with the highest mean SNP-based heritability over the considered time points was selected as a cluster representative in the following analyses. This resulted in 50 representative traits from five categories (Fig. 3e). The mean Mantel correlation between all traits was 0.63. Using this clustering method, the mean within-cluster inter-trait Mantel correlation was 0.97, while the mean inter-cluster correlation was 0.51. Investigations of the hierarchical clustering based on the distance matrix resulting from the Mantel correlation resulted in clusters with mean within-cluster inter-trait correlations approximately equal to the inter-cluster correlations, and was therefore deemed unsatisfactory (see Table S2 for a comparison of modularity and hierarchical clustering). Data for each trait were then min-max normalized across all lines and time points and were used in modeling.

Dynamic mode decomposition

Dynamic mode decomposition determines a time-invariant best-fit linear operator $\:\mathbf{A}$ that transforms measurements of data at one time point into the measurements at the following time point, i.e.

$$\:{\mathbf{x}}_{\mathbf{t}+1}=\mathbf{A}{\mathbf{x}}_{\mathbf{t}}$$

where $\:\mathbf{x}$ is a column vector of traits and $\:\mathbf{A}$ is a $\:p\times\:p$ matrix. When the $\:\mathbf{x}$ vectors of all time points are concatenated into matrices, Eq. (1) becomes

$$\:{\mathbf{X}}_{2}=\mathbf{A}{\mathbf{X}}_{1}$$

where $\:{\mathbf{X}}_{1}$ and $\:{\mathbf{X}}_{2}$ are $\:p\times\:\left(T-1\right)$ matrices which are offset by a single time point (Fig. 1b). From $\:{\mathbf{X}}_{1}$ and $\:{\mathbf{X}}_{2}$, the matrix $\:\mathbf{A}$ can be directly calculated using

$$\:\mathbf{A}={\mathbf{X}}_{2}{{\mathbf{X}}_{1}}^{†}$$

where $\:{{\mathbf{X}}_{1}}^{†}$indicates the Moore-Penrose pseudoinverse of $\:{\mathbf{X}}_{1}$. The operator $\:\mathbf{A}$ allows us to predict the columns of $\:\mathbf{X}$ beginning with $\:{\mathbf{x}}_{1}$. In fluid mechanics $\:p$ is often very large (i.e. millions) and $\:T$ is in tens of thousands, and there are many algorithms which rely on lower-dimensional representations. We cannot use these representations since they rely on eigenvalue decomposition that renders the lower-dimensional approximation symmetric. In our case, with $\:p=50$ and $\:T=25$, we can compute $\:\mathbf{A}$ directly from the data. However, our tests demonstrated that this approach is sensitive to deviations in input data. To address these issues and to improve numerical stability, we used Algorithm 2, known as Schur-Based DMD¹⁷ to find the modes and to obtain an approximation of $\:\mathbf{A}$, denoted as $\:{\mathbf{A}}_{\mathbf{r}}$, detailed below:

Algorithm 1

(Classical DMD):

$$\:{\mathbf{x}}_{\mathbf{t}+1}=\mathbf{A}{\mathbf{x}}_{\mathbf{t}}$$

$$\:{\mathbf{X}}_{2}=\mathbf{A}{\mathbf{X}}_{1}$$

$$\:\mathbf{A}={\mathbf{X}}_{2}{{\mathbf{X}}_{1}}^{†}$$

;

Algorithm 2

(Schur-based DMD):

$\:{\mathbf{X}}_{1}=\mathbf{U}\varvec{\Sigma\:}{\mathbf{V}}^{\mathbf{T}}$ , (4) singular value decomposition of $\:{\mathbf{X}}_{1}$

$\:\stackrel{\sim}{\mathbf{A}}={\mathbf{U}}_{\mathbf{r}}^{\mathbf{T}}{\mathbf{X}}_{2}{\mathbf{V}}_{\mathbf{r}}{\varvec{\Sigma\:}}_{\mathbf{r}}^{-1}$ , (5) rank-reduced representation of $\:\mathbf{A}$ projected onto the $\:r$ POD modes of $\:\mathbf{U}$

$\:\stackrel{\sim}{\mathbf{A}}={\mathbf{Q}}^{\mathbf{T}}\mathbf{R}\mathbf{Q}$ , (6) Schur decomposition of $\:\stackrel{\sim}{\mathbf{A}}$

$\:\varvec{\Phi\:}={\mathbf{X}}_{2}{\mathbf{V}}_{\mathbf{r}}{\varvec{\Sigma\:}}_{\mathbf{r}}^{-1}\mathbf{Q}$ , (7) projected DMD modes

$\:{\mathbf{A}}_{\mathbf{r}}=\varvec{\Phi\:}\mathbf{R}{\varvec{\Phi\:}}^{†}$ . (8) reconstruction of truncated $\:{\mathbf{A}}_{\mathbf{r}}$

In Schur-based DMD, $\:{\mathbf{U}}_{\mathbf{r}}$, $\:{\mathbf{V}}_{\mathbf{r}}$ and $\:{\varvec{\Sigma\:}}_{\mathbf{r}}^{-1}$, denote the matrices from the singular value decomposition of $\:{\mathbf{X}}_{1}$ restricted to the first $\:r$ columns. To account for the two-day gaps in our data we modified our $\:{\mathbf{X}}_{1}$ and $\:{\mathbf{X}}_{2}$ such that $\:{\mathbf{X}}_{1}$ omitted the time point immediately prior to the gap, and $\:{\mathbf{X}}_{2}$ omitted the time point immediately following the gap so that the difference between the columns of the two matrices was always a single time step. An R implementation of Algorithms 1 and 2 is available at https://github.com/dobby978/dynamicGP

Heritability and genomic prediction analyses

To quantify the extent to which the selected traits and all the elements of all matrices defined in Algorithm 2, can be predicted by genetic markers, we calculated the SNP-based heritability using GCTA²⁶. The genomic relatedness matrix was generated using TASSEL 5²⁷. We then investigated the accuracy of genomic prediction using ridge regression best linear unbiased prediction (RR-BLUP) model as implemented in R package “rrBLUP”²⁸ to predict individual elements of these matrices. Prediction accuracy for the elements was quantified as the mean Pearson correlation coefficient between the true and predicted values in 20 iterations of 5-fold cross validations.

Selection of $\:r$ value in Schur-based DMD

The number of singular values $\:r$ was selected by analyzing the heritability and predictability of components in Schur-Based DMD using $\:r\in\:\left[2,\:3,\:4,\:5,\:6\right]$. The heritabilities of the singular values and entries of the first two singular vectors are relatively high (0.47 and 0.35 on average, respectively, Fig. S4). Most entries of the singular vectors associated with the smaller singular values exhibited near-zero heritability (Fig. S4). For the Schur decomposition of $\:\stackrel{\sim}{\mathbf{A}}$ with $\:r=2$ ($\:{\stackrel{\sim}{\mathbf{A}}}_{2}$), the entries of $\:\mathbf{Q}$ and $\:\mathbf{R}$ exhibited non-zero heritability (0.04 and 0.29 on average, respectively), while the entries of the projected DMD modes, $\:\varvec{\Phi\:}$, had mean heritability of 0.10 (Fig. S5). We observed similar findings regarding to predictability that the mean prediction accuracies for the first two singular values and the entries of the singular vectors as well as of $\:\stackrel{\sim}{\mathbf{A}}$ were at least 0.1, while for other values of $\:r$ the prediction accuracy for the entries of $\:{\stackrel{\sim}{\mathbf{A}}}_{\mathbf{r}}$ were near-zero (Fig. S6). Similarly, the entries of the matrices $\:\mathbf{Q}$ and $\:\mathbf{R}$ obtained from $\:{\stackrel{\sim}{\mathbf{A}}}_{2}$ showed non-zero prediction accuracy, while $\:\mathbf{Q}$ and $\:\mathbf{R}$ of $\:{\stackrel{\sim}{\mathbf{A}}}_{\mathbf{r}}$ with $\:r$ larger than 2, exhibited near-zero prediction accuracy. The matrix $\:\varvec{\Phi\:}$ had non-zero mean prediction accuracy for $\:r$ value of at most 3 (Fig. S7). Therefore, we used the first two singular vectors ($\:r=2$) in the subsequent analysis, as the matrix entries corresponding to larger values of $\:r$ were not heritable nor predictable with sufficient accuracy to contribute any useful information. This also had the benefit of acting as a noise-filter which reduced overfitting.

Implementation of dynamicGP

The predicted $\:\mathbf{Q}$ and $\:\mathbf{R}$ matrices, denoted $\:{\varvec{\Phi\:}}_{\mathbf{P}}$ and $\:{\mathbf{R}}_{\mathbf{P}}$, for each line were created by collecting the entries obtained from the RR-BLUP models and placing them in their corresponding location in the matrices; $\:{\varvec{\Phi\:}}_{\mathbf{P}}$ and $\:{\mathbf{R}}_{\mathbf{P}}$ were then used in the following

$$\:{\mathbf{A}}_{\mathbf{P}}={\varvec{\Phi\:}}_{\mathbf{P}}{\mathbf{R}}_{\mathbf{P}}{{\varvec{\Phi\:}}_{\mathbf{P}}}^{†}$$

to obtain a prediction of $\:{\mathbf{A}}_{\mathbf{r}}$, defined in Algorithm 2 (see Eq. (8)). For each iteration of the 5-fold cross validation, we predicted $\:{\mathbf{A}}_{\mathbf{P}}$ operators in the testing fold using the model from data of the training folds. The $\:{\mathbf{A}}_{\mathbf{P}}$ operators were then used to predict the trait values in testing fold. The accuracies of prediction were assessed for each trait using the Pearson correlation between the predicted trait values and the measured trait values. This resulted in 100 predictions of each trait at each time point, and the mean value was reported as the final predication accuracy.

Iterative and recursive versions of dynamicGP

Trait values in the following time points were predicted in two scenarios, namely iterative and recursive. In the iterative version of dynamicGP, we used the measured trait values at each time point to predict the values at the following time point, i.e.:

$$\:{\mathbf{x}}_{\mathbf{P},\mathbf{t}+1}={\mathbf{A}}_{\mathbf{P}}{\mathbf{x}}_{\mathbf{t}}$$

while in the recursive version of dynamicGP, we used the measured trait values at $\:t=1$ and then employed the predicted values to predict the next time point over the time interval of interest, namely

$$\:{\mathbf{x}}_{\mathbf{P},\mathbf{t}+1}={\mathbf{A}}_{\mathbf{P}}{\mathbf{x}}_{\mathbf{t}}\:\text{f}\text{o}\text{r}\:t=1,$$

(11)

$$\:{\mathbf{x}}_{\mathbf{P},\mathbf{t}+1}={\mathbf{A}}_{\mathbf{P}}{\mathbf{x}}_{\mathbf{P},\mathbf{t}}\:\text{f}\text{o}\text{r}\:t>1.$$

Meuwissen, T. H. E., Hayes, B. J. & Goddard, M. E. Prediction of Total Genetic Value Using Genome-Wide Dense Marker Maps. Genetics 157, 1819–1829 (2001).
Alemu, A. et al. Genomic selection in plant breeding: Key factors shaping two decades of progress. Mol Plant17, 552–578 (2024).
Gill, H. S. et al. Multi-Trait Multi-Environment Genomic Prediction of Agronomic Traits in Advanced Breeding Lines of Winter Wheat. Front Plant Sci12, (2021).
Karaman, E., Lund, M. S. & Su, G. Multi-trait single-step genomic prediction accounting for heterogeneous (co)variances over the genome. Heredity 124, 274–287 (2020).
Montesinos-López, O. A. et al. A singular value decomposition Bayesian multiple-trait and multiple-environment genomic model. Heredity (Edinb) 122, 381–401 (2019).
Li, D. et al. High-Throughput Plant Phenotyping Platform (HT3P) as a Novel Tool for Estimating Agronomic Traits From the Lab to the Field. Front Bioeng Biotechnol8, (2021).
Song, P., Wang, J., Guo, X., Yang, W. & Zhao, C. High-throughput phenotyping: Breaking through the bottleneck in future crop breeding. Crop J9, 633–645 (2021).
Sun, D. et al. Using hyperspectral analysis as a potential high throughput phenotyping tool in GWAS for protein content of rice quality. Plant Methods 15, 54 (2019).
Zhang, H., Wang, L., Jin, X., Bian, L. & Ge, Y. High-throughput phenotyping of plant leaf morphological, physiological, and biochemical traits on multiple scales using optical sensing. Crop J 11, 1303–1318 (2023).
Poorter, H. et al. Pitfalls and potential of high-throughput plant phenotyping platforms. Front Plant Sci14, (2023).
Wang, C. et al. Predicting Plant Growth and Development Using Time-Series Images. Agronomy 12, 2213 (2022).
Yasrab, R., Zhang, J., Smyth, P. & Pound, M. P. Predicting Plant Growth from Time-Series Data Using Deep Learning. Remote Sensing13, 331 (2021).
Chang, S., Lee, U., Hong, M. J., Jo, Y. D. & Kim, J.-B. Time-Series Growth Prediction Model Based on U-Net and Machine Learning in Arabidopsis. Front Plant Sci 12, (2021).
Kutz, J. N., Brunton, S. L., Brunton, B. W. & Proctor, J. L. Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems. (SIAM, 2016).
Brunton, B. W., Johnson, L. A., Ojemann, J. G. & Kutz, J. N. Extracting spatial–temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition. J Neurosci Meth258, 1–15 (2016).
Hasnain, A. et al. Learning perturbation-inducible cell states from observability analysis of transcriptome dynamics. Nat Commun 14, 3148 (2023).
Thitsa, M., Clouatre, M., Verriest, E., Coogan, S. & Martin, C. A Numerically Stable Dynamic Mode Decomposition Algorithm for Nearly Defective Systems. IEEE Control Syst Lett5, 67–72 (2021).
Proctor, J. L., Brunton, S. L. & Kutz, J. N. Dynamic Mode Decomposition with Control. SIAM J Appl Dyn Syst15, 142–161 (2016).
Dell’Acqua, M. et al. Genetic properties of the MAGIC maize population: a new platform for high definition QTL mapping in Zea mays. Genome Biol 16, 167 (2015).
Junker, A. et al. Optimizing experimental procedures for quantitative evaluation of crop plant performance in high throughput phenotyping systems. Front Plant Sci5, (2015).
Klukas, C., Chen, D. & Pape, J.-M. Integrated Analysis Platform: An Open-Source Information System for High-Throughput Plant Phenotyping. Plant Physiol165, 506–518 (2014).
Heuermann, M. C., Meyer, R. C., Knoch, D., Tschiersch, H. & Altmann, T. Strong prevalence of light regime-specific QTL in Arabidopsis detected using automated high-throughput phenotyping in fluctuating or constant light. Physiol Plantarum176, e14255 (2024).
Ferguson, J. N. et al. The genetic basis of dynamic non-photochemical quenching and photosystem II efficiency in fluctuating light reveals novel molecular targets for maize (Zea mays) improvement. Preprint at https://doi.org/10.1101/2023.11.01.565118 (2023).
Browning, B. L., Zhou, Y. & Browning, S. R. A One-Penny Imputed Genome from Next-Generation Reference Panels. Am J Hum Genet103, 338–348 (2018).
Csardi, G. & Nepusz, T. The igraph software package for complex network research. InterJournal Compl Syst, 1695 (2006).
Yang, J., Lee, S. H., Goddard, M. E. & Visscher, P. M. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet88, 76–82 (2011).
Bradbury, P. J. et al. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635 (2007).
Endelman, J. B. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome4, 250–255 (2011).

There is NO Competing Interest.

dynamicGPsupplementaryinformation.pdf

Download PDF

Version 1

posted

You are reading this latest preprint version

Dynamics of plant phenome can be accurately predicted from genetic markers

Status:

Version 1

Abstract

Figures

Introduction

Results

DMD accurately forecasts time-resolved phenotypes

Building blocks of dynamicGP are heritable and can be predicted from genetic markers

Differences in performance between the two versions of dynamicGP

Performance of dynamicGP compared to baseline models

Discussion

METHODS

Data generation

Clustering and trait selection

Dynamic mode decomposition

Heritability and genomic prediction analyses

Selection of \(\:r\) value in Schur-based DMD

Implementation of dynamicGP

Iterative and recursive versions of dynamicGP

(11)

References

Additional Declarations

Supplementary Files

Status:

Version 1