We used four Japanese cultivars as parents of the Japan-MAGIC2 (JAM2) lines: Iwaidawara (IW), Toyomeki (TO), Akidawara (AK), and Tachiharuka (TH). We first crossed AK with TO and IW with TH to produce seeds of two types, called the AKTO and IWTH two-way hybrids. Then these hybrids were crossed to produce four-way recombinants. We finally produced 100 JAM2 lines by the single-seed descent (SSD) method. These JAM2 lines (F5 in 2022, F6 in 2023) were used in this study (Supplementary Figure 1A).

We cultivated the 100 JAM2 lines in a rice field in Tsukuba, Japan, in 2022 and 2023 using both regular (R) and delayed (D) transplanting protocols (Figure 1; Supplementary Table 1). The dates of sowing and transplanting to the paddy field in the delayed transplanting protocol were about a month later than those in the regular transplanting protocol (Figure 2). In this study, the data from 2022 were used for the main analysis, and the data of 2023 were used for evaluating the robustness of the trait prediction models.

Days to heading (DTH) was scored as the number of days from transplanting rice to the field to the appearance of the first panicle in more than half of the plants in each JAM2 line. CL was assessed as the length of the longest culm of each plant measured with a ruler more than 10 days after heading. For the measurement of panicle weight (PW) and stem and leaf weight (SLW), shoots of mature plants were air-dried from one to two months in a drying room and then cut 3 cm below the panicle base to separate the parts. The aboveground dried weight (ADW) was calculated as PW + SLW. The averages of five plants in the middle lane except for the plants at both edges were used as the CL, PW, SLW, and ADW values of each JAM2 line (Figure 1).

Similar to our previous studies (Ogawa et al., 2021b; Taniguchi et al., 2022), UAV-based aerial photography was conducted once a week to track the growth of the rice lines. For lines grown under both the regular transplanting and delayed transplanting protocols, aerial photography was conducted once a week after puddling and before the planting date until harvest, which occurred by 1 November. For the aerial photography, we used a Phantom 4 Pro (P4P) UAV equipped with an RGB camera to obtain data for calculating CH, and a DJI Phantom 4 Pro Multispectral (P4M) UAV, which is equipped with multispectral cameras, to obtain data for calculating CIg (Figure 1). The automatic flight settings for the P4M were slightly different from those for the P4P in order to shorten the flight time, thereby reducing the impact of short-term changes in incident light intensity on the spectral images (Ogawa et al., 2021a). The flight path and image shooting setting were programmed by using DJI GS Pro software. The detailed settings are in Supplementary Table 2. The band setting of P4M were as follows: blue is 450nm±16nm, green is 560nm±16nm, red is 650nm±16nm, red edge is 730nm±16nm and NIR is 840nm±16nm.

The P4P was manually raised to 10 m and the camera’s focus distance was adjusted on a region of the canopy before the aerial photography, then was started in automatic flight mode at an altitude of 10 m. The P4M was manually operated to a position 1–2 m above a Micasense Calibration Reflectance Panel (MicaSense Inc., Seattle, WA, USA) to capture spectral images of the gray plate with 50% reflectance before and after operation in automatic flight mode at an altitude of 20 m. We set painted black and white markers on paved surfaces at eight points surrounding the field as ground control points (GCPs) and then precisely measured the latitude, longitude, and altitude of each point with a TCRP1205 surveyor (Leica, Heerbrugg, Switzerland).

We obtained multispectral orthomosaic images and a crop surface model (CSM) from each set of aerial images with Agisoft MetaShape Professional v. 1.7.3 software (Agisoft, St. Petersburg, Russia). The CSM was generated from the high-resolution RGB images acquired by the P4P using the same steps as described previously (Ogawa et al., 2019): (1) Align photos (accuracy, high), (2) input GCPs, (3) build dense cloud (accuracy, high), (4) build mesh (surface type, height field; source data, dense cloud), and (5) build digital elevation model (DEM; source data, dense cloud). The DEM image of 16 May was defined as the height of the ground surface. Then, a CSM image representing the height of the rice plants was created by taking the difference between the ground surface height image obtained on another observation date and the DEM image. The multispectral orthomosaic images were generated from the P4M aerial images with an alternative procedure to reduce misalignment of the spectral images (Sakamoto et al., 2022). The following steps were repeated for each spectral image: (1) Align photos (accuracy, high), (2) input GCPs, and (3) calibrate reflectance using the sun sensor data and gray panel images. The sparse-point data and the GCP data for each spectral dataset were merged into a single chunk in the Meta Shape to perform camera calibration. Then, the following steps were conducted: (4) Build DEM (source data, sparse cloud) and (5) build orthomosaic images (surface DEM; blending mode, mosaic). The orthomosaic images and the DEM images were analyzed in ENVI v. 5.5 remote sensing software (Harris Geospatial, Boulder, CO, USA). The original spatial resolution was 3 mm/pixel for P4P and 11 mm/pixel for P4M. The map projection was converted to UTM zone 54 N (WGS-84), and then both DEM and orthomosaic images were resampled with a 2-mm/pixel resolution. The resampling resolution was set to a slightly higher resolution than P4P RGB camera images because of accounting for variations in ground resolution due to changes in actual UAV flight altitude. The converted image was rotated 66° clockwise to match the direction of the long side of the field with the lateral direction of the final output image. The image was then resized to a rectangle (28,000 pixels × 14,000 pixels).

JAM2 lines were cultivated in a small plot in a field: the distance between adjacent plants was 30 cm between columns and 18 cm between rows. The plots in which each line was planted were cut out by using ENVI software from the orthomosaic and DEM images. Each cut-out image corresponded to 90 cm × 126 cm on the ground and contained 21 plants (Supplementary Figure 1B). From the cut-out images, CH and CIg were calculated as follows. We defined CH as the difference between the canopy position, that is, the 95th percentile value in the cut-out DSM, and the ground level; the method used to obtain values were consistent with manually measured CH data (Taniguchi et al., 2022). CIg was calculated as the ratio of NIR to green reflectance values (Gitelson et al., 2003; Gitelson et al., 2005):

We defined the CIg of each line as the mean CIg in the plot corresponding to that line in the orthomosaic image.

We fitted time-series models to the CH and CIg time-series data to obtain the CH parameters and CIg parameters (Figure 1). In our previous study (Taniguchi et al., 2022), in which we analyzed time-series CH data, we applied a modified three-parameter logistic model. The three-parameter logistic model accounted for the CH decrease in the late growth period by using a quadratic curve,

where x is the days after transplanting. The parameters of the three-parameter logistic model were estimated by the same procedure of our previous study, the algorithm of which was implemented in the R package phenolocrop (Taniguchi et al., 2022). Parameter K is the maximum value of CH, r1 is the growth rate before the peak, d0 is the time point at which the growth rate is a maximum, d1 is the time point at which the maximum value of y is reached, and a is the rate at which CH decreases in the late growth period (Figure 3A). Different from CH, CIg time-series data described an S-shape with time; therefore, we adopted a double logistic model for the CIg data (Fisher and Mustard, 2007; Yang et al., 2012).

Here, parameter r2 is the growth rate before the peak, r3 is the rate of decrease after the peak, d2 is the time point at which the growth rate is a maximum, d3 is the time point at which the decreasing rate is a maximum, and ymax is the maximum value of CIg (Figure 3B). Here, we fitted both a three-parameter logistic model and a double logistic model to the time-series CH and CIg data. To estimate the parameters of the double logistic model, we first set the maximum value of the objective time-series variable (CH or CIg) to ymax and then estimated the other variables in the framework of the nonlinear least squares method implemented in the R function nls. The algorithm was set to nl2sol, and the initial values of r2, r3, d2, and d3 were 0.05, 0.05, 40, and 100, respectively.

To evaluate which model was best suited to CH or CIg, we calculated the coefficients of determination (R2) for each observation dataset. For CH, the R2 values were high for both the three-parameter logistic model (R2=0.982 on average) and the double logistic model (R2=0.985 on average). While the double logistic model was slightly better, we adopted the three-parameter logistic model for the time-series CH data because our objective is to evaluate the robustness of the methodology presented in our previous study (Taniguchi et al., 2022). In contrast, for CIg, we adopted the double logistic model because the R2 value (R2=0.941 on average) for that model was larger than that for the three-parameter logistic model (R2=0.796 on average; Supplementary Figure 2).

All of the manually measured traits and model parameters were characterized by summary statistics (mean, variance, max value, and min value), frequency distribution, and the calculated Pearson’s correlation coefficient (cor) between each trait and parameter. To determine how the CH and CIg parameters were related, we conducted principal component analysis (PCA).

For the prediction of manually measured traits from CH and/or CIg parameters, we conducted variable selection and constructed six linear regression models depending on which parameters were used as the predictors (Figure 1). The CH-selected model used K, d₀, and d₁ as predictors. These parameters were selected from the five CH parameters by a procedure called backward variable selection to prevent multicollinearity (Hastie et al., 2009). When the variance inflation factor (VIF) was calculated using the car package in R (Fox and Weisberg, 2019), the selected parameters, K, d₀, and d₁, had VIFs lower than 5 (Supplementary Table 3). The CIg-selected model used r2, d2, d3, and ymax as predictors without multicollinearity. These parameters were selected from the five CIg parameters in the same way (Supplementary Table 4). Then, to investigate whether the two models contained different information for trait prediction, we constructed a “Hybrid-selected model,” which used all of the predictors used by the CH-selected and CIg-selected models. We also constructed the corresponding full models (CH-full model, CIg-full model, and Hybrid-full model) to examine the effects of variable selection.

We evaluated the performance of the six prediction models using four validation schemes: model fitness, model accuracy, type-1 model robustness, and type-2 model robustness (Figure 4). To evaluate model fitness, each prediction model was fitted separately to the JAM2 data of R and D in 2022, and their coefficients of determination (R2), regression coefficients, and p-values were calculated to measure the goodness of fit and to identify those parameters that were important for trait prediction. To evaluate model accuracy, we conducted 10-fold cross validation (10-CV) for R and D separately in 2022. Because cross-validation results can fluctuate depending on the data-splitting process, we randomly repeated the 10-CV 100 times and calculated root mean squared errors (RMSEs) and cor values. We also compared prediction accuracies among the CH-selected, CIg-selected, and Hybrid-selected models by conducting a “Tukey-like” non-parametric multiple comparison among them with the nparcomp function implemented in the nparcomp package in R (Konietschke et al., 2015). To evaluate type-1 model robustness, we trained the prediction model with R data and tested the trained model against D data, and vice versa. To evaluate type-2 model robustness, we trained the prediction model with R data of 2022 and tested the trained model against R data of 2023. We also trained the prediction model with D data of 2022 and tested the trained model against D data of 2023. We quantified the two types of model robustness by calculating the RMSE and the cor values.

After predicting the phenotypic values of the test data in the evaluation procedure of model robustness, we calibrated the predicted phenotypic values by using the data of the parental cultivars of JAM2. The calibration model was constructed in the following five steps (Supplementary Figure 3).

As shown in Supplementary Figure 3A, for the type-1 model robustness, (1) we trained the prediction model with D data of JAM2. (2) We applied the trained prediction model to R data of JAM2 and obtained the predicted values of R data of JAM2. (3) We applied the trained prediction model to R data of parental cultivars of JAM2 and obtained the predicted values of the parental cultivars. (4) We trained a single regression model as the calibration model by comparing the predicted and observed values of the parental cultivars. (5) We applied the trained calibration model to the predicted values of the R data of JAM2 and obtained the calibrated values. Similarly, for the case of training the prediction model with R data of JAM2, we calibrated the predicted values of D data of JAM2. We quantified the calibration by calculating the RMSE and the cor values.

For the type-2 model robustness, as shown in Supplementary Figure 3B, (1) we trained the prediction model with R data of JAM2 in 2022. (2) We applied the trained prediction model to R data of JAM2 in 2023 and obtained the predicted values of R data of JAM2 in 2023. (3) We applied the trained prediction model to R data of the parental cultivars of JAM2 in 2023 and obtained the predicted values of the parental cultivars. (4) We trained a single regression model as the calibration model by comparing the predicted and observed values of the parental cultivars. (5) We applied the trained calibration model to the predicted values of R data of JAM2 in 2023 and obtained the calibrated values. Similarly, for the case of training the prediction model with D data of JAM2 in 2022, we calibrated the predicted values of D data of JAM2 in 2023. These calibration frameworks were applied to the Hybrid-selected model; the parameter and trait data are available in Supplementary Data 1.

We examined the distributions of five manually measured traits of the JAM2 lines grown under the regular and delayed transplanting protocols. The ranges of trait distribution under the regular transplanting protocol were as follows: days to heading (DTH: 59–100 days), culm length (CL: 51–110 cm), aboveground dried weight (ADW: 56–118 g), stem and leaf weight (SLW: 24–76 g), and panicle weight (PW: 21–58 g). Those under the delayed transplanting protocol were as follows: DTH (48–84 days), CL (58–121 cm), ADW (69–114 g), SLW (23–70 g), and PW (28–59 g). We found that the CL, ADW, SLW, and PW distributions, which were investigated at the maturation stage, were largely comparable between the regular and the delayed transplanting protocols. By contrast, the DTH distribution was clearly shorter under the delayed transplanting protocol (48–84 days) than under the regular transplanting protocol (59–100 days). Thus, the vegetative growth period was shortened under the delayed transplanting protocol (Figure 5A; Supplementary Table 5).

Positive correlations were detected between the regular and delayed transplanting protocols for the five manually measured traits, with correlation coefficients ranging from 0.98 (DTH) to 0.50 (PW: Supplementary Figure 4). In the case of DTH, which had the highest correlation coefficient, the transition from the vegetative phase to the reproductive phase was consistently earlier under the delayed transplanting protocol. The probable reason is that the JAM2 parental cultivars are adapted to the Japanese environment, so the transition occurs when the plants detect high temperatures together with a shorter day length at the end of July (Figure 2), conditions that are known to promote heading (Vicentini et al., 2023). The correlation coefficient of CL (0.96) was also high, but the difference between the transplanting protocols in this trait was small. The correlation for PW was small (cor = 0.50) but still significantly positive (Supplementary Table 4). The correlations for ADW and SLW were between these values at 0.65 and 0.81, respectively.

Among the five traits, DTH was consistently highly correlated to SLW (cor > 0.7; Supplementary Figure 5). The correlation of CL with SLW was smaller but still positive (about 0.6). Under the regular transplantation protocol, PW showed almost no correlation with DTH or SLW, whereas under the delayed transplanting protocol, PW was negatively correlated with DTH and SLW. Thus, PW exhibited distinctive characteristics under the delayed protocol.

Overall, these results indicate that these five traits have distinctive characteristics in terms of influence by different transplanting protocols and correlations among traits. It is therefore reasonable to focus on the growth pattern of the JAM2 lines during cultivation to predict the five manually measured traits.

We first focused on the vertical growth pattern through time-series CH monitoring. From temporal CH data of the JAM2 lines, we extracted the CH parameters from the time-series image data from the UAVs, which reflect the vertical growth of rice and consist of a, d0, d1, r1, and K (Figure 3A). The frequency distributions of the growth speed r1, the time point of maximum growth speed d0, and the time point of maximum CH d1 in the modified three-parameter logistic model of time-series CH data differed notably between the two protocols (Figure 5B).

Next, we focused on CIg, an index of the total chlorophyll content of the canopy that reflects vegetation amount and senescence. The pattern of the CIg curve was similar to that of the CH curve, but the maximum CIg was reached at an earlier date than the CH maximum (Figures 3A,B). We fitted a double logistic model to the time-series CIg data of the JAM2 lines to determine the five CIg parameters r2, r3, d2, d3, and ymax. The frequency distributions of all CIg parameters differed between the regular and delayed transplanting protocols, but the distributions of growth rate r2 and the time point of the maximum growth rate d2 differed more remarkably between the protocols than those of the other CIg parameters (Figure 5C). These results indicates that early growth might be sensitive to the cultivation protocol.

Focusing on PC1 in the PCA (Figure 6A), we found three parameter clusters common to both the regular and delayed transplanting protocols: cluster I (consisting of d0, d1, and d3), cluster II (K and ymax), and cluster III (r1, r2, and r3). The parameters in each cluster were positively correlated (Figure 6B). Notably, the parameters of clusters I and III, which reflect key time points and rates of growth and development, respectively, were negatively correlated. Parameter d2 was also negatively correlated with r2. These results reflect the tendency for faster growth to be associated with a shorter period of growth and development and vice versa. The reproductive phase parameters a and r3 were weakly correlated. In addition, K and ymax, both in cluster II, were weakly correlated but contained different information, as can be seen by examining PC3 and PC4 (Figure 6A). Considered together, these results show that the CH and CIg parameters had both similar and distinct characteristics.

We evaluated the model fitness based on the goodness of fit of the prediction models to the training data by considering R2 (Figures 4, 7). Furthermore, we examined which CH parameters had large effects on model fitness. With regard to the prediction of DTH, R2 was between 0.7 and 0.8 for both protocols (Supplementary Table 6). The coefficient of d1 on DTH had the largest absolute value and was significantly positive (Figure 7A; Supplementary Tables 7, 8). For the prediction of CL, R2 was more than 0.9 under both protocols. The coefficient of K on CL had the largest absolute value and was significantly positive. For the prediction of SLW and ADW, R2 was between 0.5 and 0.8 under both protocols. The coefficient of d1 on SLW had the largest absolute value and was significantly positive. For the prediction of ADW, the absolute value of the coefficient of d1 was largest under the regular transplanting protocol, whereas the effect of K was largest under the delayed transplanting protocol. For the prediction of PW, R2 was close to zero for all CH parameters (Figure 7B; Supplementary Table 6). These results indicate that the model using CH parameters was fit to predict DTH, CL, ADW, and SLW but not PW. More specifically, d1 and K contained rich information for the prediction model to fit the data of DTH, CL, ADW, and SLW in the JAM2 lines, but there were no CH parameters to fit the data of PW.

We also asked if CIg would be useful for the prediction models to fit to the training data of manually measured traits. Since CIg had a higher correlation to PW than NDVI, especially during the early growth period (Supplementary Figure 6), we adopted CIg as the focal VI. With regard to the prediction of DTH, R2 was consistently greater than 0.7 and the absolute value of the coefficient of d3 was largest and significantly positive (Figure 7A; Supplementary Tables 6, 9, 10). For the prediction of CL, R2 was 0.17 for the regular transplanting protocol and 0.39 for the delayed transplanting protocol; therefore, CIg parameters were not appropriate for predicting CL. For the prediction of ADW and SLW, R2 was between 0.6 and 0.8 for both protocols, and the absolute value of the coefficient of ymax was largest for the regular transplanting protocol, whereas that of r2 was largest for the delayed transplanting protocol. For the prediction of SLW, the coefficient of d3 consistently had the largest absolute value and was significantly positive. For the prediction of PW, R2 was 0.44 for the regular transplanting protocol and 0.35 for the delayed transplanting protocol. The coefficient of ymax had the largest absolute value and was significantly positive. The coefficient of d3 on PW had the second largest absolute value and was significantly negative. These results indicate that the model using CIg parameters was fit to predict DTH, ADW, SLW, and PW but not CL. Specifically, r2, ymax, and d3 in the CIg parameters contained information for the prediction model to fit the data of DTH, ADW, SLW, and PW in the JAM2 lines, but there were no CIg parameters to fit the data of CL.

To improve predictions of the five manually measured traits, we attempted to use both CH and CIg parameters for prediction (Hybrid-selected model). For all five traits obtained under regular and delayed transplanting protocols, the R2 of the Hybrid-selected model was higher than that of either the CH-selected or the CIg-selected model (Figure 7B; Supplementary Table 6). This result indicates that the use of CH and CIg parameters together increases model fitness. We also compared CH-selected, CIg-selected, and Hybrid-selected models with CH-full, CIg-full, and Hybrid-full models. The effects of variable selection on model fitness are described in the Discussion.

Model accuracy was evaluated based on the ability to predict test data obtained under the same conditions (transplanting protocol and year) as the training data (Figure 4). The model accuracy of each trait by the CH-selected, CIg-selected, and Hybrid-selected models was consistent with the model fitness both in terms of cor and RMSE between predicted and observed values (Figure 8; Supplementary Figure 7 and Supplementary Table 11). For the prediction of PW, the CIg-selected model showed significantly better performance than CH-selected model (Supplementary Figure 7).

Model robustness was evaluated based on the ability to predict test data obtained under different transplanting protocols (type-1) or different years (type-2) from the training data (Figure 4). For DTH, CL, and SLW, the correlations between predicted and observed values in type-1 and -2 model robustness were consistent, but RMSEs tended to fluctuate and become larger (Figure 9). For example, when DTH was used to evaluate type-1 model robustness, CH-selected, CIg-selected, and Hybrid-selected models had similarly high cor values, whereas the RMSEs were different; CIg-selected and Hybrid-selected models had larger RMSEs than the CH-selected model (Supplementary Table 12), possibly because the values predicted by the former two models are biased to larger values (Supplementary Figure 8).

In the case of the prediction of ADW using CH-selected, CIg-selected, and Hybrid-selected models, cor values were high in terms of type-1 model robustness but not in terms of type-2 model robustness (Supplementary Tables 12, 13). For the prediction of PW by CIg-selected and Hybrid-selected models, which had high model fitness and accuracy, cor values were not always high (Supplementary Table 13). A bias in the predicted values of ADW and PW was also observed (Supplementary Figure 8).

In summary, except for the prediction of PW and ADW in some cases, the prediction models with high model accuracy also had high model robustness in terms of cor, but they also tended to have larger RMSEs because of prediction biases derived from different transplanting protocols and years. The effect of variable selection is considered in the Discussion.

To mitigate the problem of large RMSEs when the training and test data were from different transplanting protocols or years, we evaluated the improvement by a calibration procedure using training and test data (Figure 10). In terms of type-1 model robustness, the calibration reduced the RMSE in predicting DTH by the Hybrid-selected model from 11.3 to 6.81 (trained using R data) and from 16.7 to 5.90 (trained using D data; Supplementary Table 14). The calibration also reduced the RMSE values for the prediction of ADW, SLW, and PW (Supplementary Table 14). There was one case where calibration resulted in a larger RMSE: the prediction of CL under the delayed transplanting protocol. In terms of type-2 model robustness, calibration reduced the RMSE values for PW (Supplementary Table 15), but the calibration did not always work well for CL, ADW, DTH and SLW. In these cases, the phenotypic data of the four parental cultivars did not cover the full range of phenotypic variance of the JAM2 lines.