In this study, the accessions and experimental fields were the same as in the previous study (Sakurai et al., 2022). The diverse panel of 178 soybean accessions used were obtained from the gene bank of the National Institute of Agrobiological Sciences, Tsukuba, Japan ( Table S1 ). These 178 accessions mainly consisted of Japanese and global soybean minicore collections (Kaga et al., 2011; Kajiya-Kanegae et al., 2021). A total of 178 soybean accessions were grown in three years 2019, 2020, and 2021 in the same experimental field at the Arid Land Research Center, Tottori University, Japan (35°32’ N lat, 134°12’ E long, 14 m above sea level). The experimental field soil was sandy and retained high water permeability. We used four types of drought treatments: no watering (drought, Treatment D), watering for 5 d followed by no watering for 5 d (Treatment W5), watering for 10 d followed by no watering for 10 d (Treatment W10), and no watering for 10 d followed by watering for 10 d (Treatment D10). Treatments D, W5, and D10 were implemented in 2019. Treatments D, W10, and D10 were implemented in 2020 and 2021 ( Figure 1 ). Each treatment consisted of two ridges. Each ridge consisted of two rows and one irrigation tube, and microplots were placed parallel to each other on either side of the tubes. The ridge height and width were 30 and 136 cm, respectively. For each treatment, 178 accessions were randomly assigned to microplots per year. As each accession was assigned one microplot per treatment, there were no replicates per treatment. The reason for not taking replications but prioritizing the number of genotypes is that previous studies have shown that genomic predictive modeling and QTL analysis are better with a higher number of genotypes, even if the number of replications is set to one (Knapp and Bridges, 1990; Lorenz, 2013). Four plants of each accession were grown in each microplot. The distances between the rows, microplots, and plants were 50, 80, and 20 cm, respectively ( Figure 2 ). Fertilizers (13, 6.0, 20, 11, 7.0 g/m2 115 of N, P, K, Mg, and Ca, respectively) were applied to the field prior to sowing. Sowing was performed on 10 July 2019, 8 July 2020, and 6 July 2021. Two to three seeds were sown at each position, after which the germinated seedlings were thinned to one per position 2 weeks after sowing.

White mulch sheets (Tyvec; DuPont, Wilmington, USA) were laid over the ridges to prevent rainwater infiltration into the soil and control soil drought levels ( Figure 1 ). A watering tube was installed under the mulch sheets at the center of each row. The watering tube (JKC Agro, Kumamoto, Japan) irrigated at a flow rate of 1.1 L/h m. Watering was done for over 5 h daily (7:00–9:00, 12:00–14:00, and 16:00–17:00), starting the day after seedling thinning for Treatments W5, W10, and D10. The irrigation cycle for each treatment is shown in Figure S1 . Soil moisture was measured using a soil moisture meter (TDR-341F; Fujiwara Seisakusho, Tokyo, Japan) at 10 sites for each treatment over 30 d. Except for rainy days, no more than 3 d were allowed between soil moisture measurements. Soil moisture in each treatment was 127 calculated as the average of the treatments in the year ( Figure S1 ). Days to flowering (DTF), which was defined as the date that 50% of plants flowered in each microplot ( Figure S2 ), were determined. To record DTF, we visited each microplot every other day during the flowering period. In each microplot of the four plants, two centrally located plants were collected 62 d after sowing, and the fresh aboveground weight of each plant was measured. The phenotypic value of the fresh weight of each microplot was calculated as the average of the two plants measured for each microplot. The number of rainy days was 19, 12, and 22 during the field experiment in 2019, 2020, and 2021, respectively.

The MS image collection and image analysis referred to the method employed by Sakurai et al. (Sakurai et al., 2022). In each treatment of a year, MS images were collected using unmanned aerial vehicles (UAVs). In 2019 and 2020, MS images (1.0 cm/pixel) were collected using a four-eye MS camera (Xacti, Osaka, Japan) mounted on a quadcopter UAV (DJI Matrice 100; DJI, Shenzhen, China). The MS camera has four independent lenses and sensors attached to different filters (MidOpt, Palatine, USA), including a triple-bandpass filter (TB550/660/850) and a red-edge bandpass. The TB550/660/850 can collect spectral intensities at 550 nm (green), 660 nm (red), and 850 nm (near-infrared). Bi725 can collect the spectral intensity at 725 nm (red-edge). The MS camera was set for continuous data capture at two frames per second per lens for a total of eight frames per second. The overlap and sidelap rates were set to 90% and the flights were set to an altitude of 20 m. In 2021, MS images (0.74 cm/pixel) were collected using DJI Phantom 4 Multispectral (P4M; DJI, Shenzhen, China). P4M has one RGB sensor and five spectral-band sensors at 450 nm (blue), 560 nm (green), 650 nm (red), 730 nm (red-edge), and 840 nm (near-infrared). P4M continuously collected MS images every two seconds during each flight. The overlap and sidelap rates were set to 75% and the flights were set to an altitude of 15 m.

Each spectral reflectance was calculated as the ratio of each spectral intensity from a grey scale panel set in the experimental field ( Figure 1 ). All the flights were scheduled for 11:00-13:00 under clear sky conditions. Images were collected four, six, and seven times in 2019, 2020, and 2021, respectively ( Figure 3 ). UAVs measurements were scheduled before and after the irrigation treatment was switched on W10 and D10 to capture the response of plants to changes in irrigation. Orthomosaic images of each spectral reflectance were obtained using a Pix4Dmapper (Pix4D, Prilly, Switzerland). Sixteen ground control points (GCPs) were set up in the field annually. The positions of the GCPs were measured using Aeropoints (Propeller, Sydney, Australia). Using geolocation information from the GCPs, 792 microplots, each with a maximum of four plants, were segmented from the orthomosaic images. Plants were segmented from an image of each microplot using the NDVI-based segmentation method to extract and analyze the MS image data from only the plants. In this experimental field, it was reported that the NDVI differed significantly from the soil surface, white mulch sheets, and plants in 2019 (Sakurai et al., 2022); thus, plants were segmented with the NDVI threshold set to 0.15, as applied by Sakurai et al. (Sakurai et al., 2022). These analyses (extraction of spectral reflectance, calculation of NDVI, and NDVI-based segmentation) were performed using the OpenCV v3.3.1, library in Python v3.6.8.

Based on the spectral reflectance values of each plant pixel segmented in each microplot, we calculated two types of VIs: NDVI (Price and Bausch, 1995; Zarate-Valdez et al., 2015) and NDRE (Jorge et al., 2019). The equations for these VIs are listed in Supplementary Table S2 . As there were four plants in each microplot, MS data were collected as a community of four plants within each microplot without considering the overlap between plants. Although we attempted to segment only the plant pixels, background noise was not completely removed. In each microplot, the average VI value may have been heavily influenced by the background noise. To reduce this effect, the median of all the segmented plant pixels was used as the representative value of each microplot. The VIs were calculated using R v4.1.2.

The genome dataset was the same as that used by Sakurai et al. (Sakurai et al., 2022). Whole-genome sequence data were obtained for all accessions (Kajiya-Kanegae et al., 2021). All accessions were genotyped using an Illumina HiSeq X Ten or HiSeq 4000 (Illumina, San Diego, USA), and 4,776,813 single-nucleotide polymorphisms (SNPs) were identified. Among these SNPs, those that were heterozygous or those in which >95% of the individuals had missing data were excluded. Markers were also filtered for a minor allele frequency<0.025, and missing data were imputed based on the mean, after which they were filtered again for a minor allele frequency<0.05. Finally, linkage disequilibrium was computed only for SNP pairs for which the distance was<25,000 base pairs, and SNPs with linkage disequilibrium below 0.8 were selected, resulting in 173,583 SNP markers. Using these SNP markers, the additive numerator relationship matrix G was estimated using the ‘calcGRM’ function in the ‘RAINBOWR’ package in R v0.1.25 (Hamazaki and Iwata, 2020).

The MS data were collected at four, six, and seven time points in 2019, 2020, and 2021, respectively. The ratio of flowering at the timings of UAV measurements in each combination of treatments was calculated using DTF data ( Table S3 ). A time series of VI values was modeled using an RRM (Mrode, 2014) for each treatment of each year separately. The RRM takes the following form:

where y i t is the phenotypic value of each VI (NDVI or NDRE) at the time point t   ( t = 1 , … , 4 in 2019, t = 1 , … , 6 in 2020, t = 1 , … , 7 in 2021) for genotype i   ( i = 1 , … , 178 ) , b t is the fixed effect of each time point, n r is the order of Legendre polynomial for the genetic effects, ϕ k ( t ) is the k th ( k = 0 , … , n r ) Legendre polynomials for time point t , u i k is the genetic effect for the k th coefficients of Legendre polynomials, and e i t is the random residual effect. Vector u k = ( u 1 k , … , u 178 k ) T , u = ( u 0 T , … , u n r T ) T follows the multivariate normal (MVN) distribution, u ∼ MVN ( 0 , Q ⊗ I 178 ) where Q is ( n r + 1 ) × ( n r + 1 ) (co)variance matrix for the Legendre polynomials and I 178 is 178 × 178 identical matrix. Vector e k = ( e 1 k , … , e 178 k ) T , e = ( e 0 T , … , e n r T ) T follows the MVN distribution e ∼ MVN ( 0 , I 178 ( nr + 1 ) σ e 2 ) , where I 178 ( n r + 1 ) is the identical matrix and σ e 2 is the residual variance. To determine the order of n r , we built RRMs with three values of n r   ( n r = 0 , 1 , 2 ) using the data collected in 2021, because the number of time points was the largest in 2021 among other years. The goodness of the model fit was assessed by computing Akaike’s information criterion (AIC) (Akaike, 1974). The best value of n r is selected based on the lowest AIC value. This RRM model was built using the “ASREML” R package v4.1.0.154 (Gilmour et al., 2015).

We used CV (Francis and Kannenberg, 1978) of fresh weight as an indicator of plant’s stability under drought stresses. CV was calculated as follows:

where μ i is the mean fresh weight in all nine combinations of treatments by years for genotype i   ( i = 1 , … , 178 ) , σ i is the standard deviation of the fresh weight in the combinations. A low CV value indicates high environmental stability. Before calculating μ i and σ i , fresh weight was scaled by min-max normalization, ranging from 0 to 1 for each combination of treatments by years. This is because C V i is significantly affected by specific combinations of treatments by years, which have large mean values.

Simple genomic prediction models were built to calculate the genomic heritability (Litchfield et al., 2015) for three traits; fresh weight in each combination of treatments by years, the genetic random regression coefficient of RRMs, and CV. Also, a simple genomic prediction model was built to predict the genetic value of the CV. The simple genomic prediction model has the following form:

where y i is the phenotypic value for each trait of genotype i   ( i = 1 , … , 178 ) , μ is the overall mean, u i is the genetic random effect, and e i is the residual random effect. The vector u = ( u 1 , … , u 178 ) T follows the MVN distribution, u ∼ MVN ( 0 , G σ g 2 ) where G is the additive numerator relationship matrix, and σ g 2 is the additive genetic variance. Vector e = ( e 1 , … , e 178 ) T follows the MVN distribution e ∼ MVN ( 0 , I 178 σ e 2 ) , where I 178 is the identical matrix and σ e 2 is the residual variance. Based on the estimated parameters of genetic and residual variances, genomic heritability was calculated as h 2 = σ g 2 σ g 2 + σ e 2 . The model was implemented using the ‘EMM.cpp’ function in the ‘RAINBOWR’ package in R v0.1.25 (Hamazaki and Iwata, 2020).

For each combination of treatments by years, the MTM was built as a Bayesian multivariate Gaussian model (Montesinos-López et al., 2016; Lado et al., 2018) to predict the CV of a genotype. An MTM using M secondary traits takes the following form:

where y i 0 is the phenotypic value of the CV for genotype i   ( i = 1 , … , 178 ) ; y i m   ( m = 1 , … , M ) is the phenotypic value of the secondary traits; μ m is the overall mean for trait m   ( m = 0 , … , M ) ; u i m is the genetic random effect ( m = 0 , … , M ) ; and e i m is the residual random effect ( m = 0 , … , M ) . Vector u m = ( u 1 m , … , u 178 m ) T , u = ( u 0 T , … , u M T ) T follows the MVN distribution, u ∼ MVN ( 0 , Σ ⊗ G ) where Σ is ( M + 1 ) × ( M + 1 ) genetic (co)variance matrix across traits and G is the additive numerator relationship matrix. Vector e m = ( e 1 m , … , e 178 m ) T , e = ( e 0 T , … , e M T ) T follows the MVN distribution e ∼ MVN ( 0 , R ⊗ I 178 ) , where R is ( M + 1 ) × ( M + 1 ) residual (co)variance matrix across traits and I 178 is the identical matrix.

We built two different types of MTMs: (1) an MTM directly using the VI value at each time point as a secondary trait ( MT All   model ), and (2) an MTM using the genetic random regression coefficients of RRM (Equation 1) as a secondary trait ( MT RRM   model ).

We assumed three cases using MTM ( Figure 4 ) and compared the prediction accuracies of MT All   model and MT RRM   model . The prediction accuracy was evaluated using a 10-fold cross-validation with 10 replicates. Pearson’s correlations were calculated between the observed and predicted CV values in each replicate, and the average of these correlations was used as the prediction accuracy.

Fresh weight varied among the nine combinations of treatments by years ( Figure 5 ). As the fresh weight in Treatment D was smaller than that in the other treatments each year, Treatment D was the treatment with the most severe drought stress. For Treatment D10, the fresh weight in 2019 was higher than that in the other two years. Genomic heritability of fresh weight ranged from 0.18 to 0.64 and that of CV was 0.35 ( Table 1 ). Next, we calculated the phenotypic correlation between CV and fresh weight for each combination of treatments by years. No positive phenotypic correlation was observed between the CV and fresh weight. In addition, the phenotypic correlation between CV and fresh weight was negative in seven out of nine combinations of treatments by years. This result indicates a positive relationship between fresh weight and the stability under drought stress as assessed by CV and in the seven combinations of treatments and years.

We evaluated the goodness-of-fit of the RRMs using NDVI and NDRE values for each treatment in 2021. Under all treatments, the best model was based on NDVI values using linear Legendre polynomials, that is, the order of n r was 1 ( Table 2 ). The model with n r = 0 exhibited the best NDRE values under all treatments ( Table S4 ). This result indicates that only the intercept varied among genotypes in the RRM of NDRE values. Therefore, we employed the RRM of NDVI using n r = 1 in later analyses. We built each RRM (Equation 1) using the time-series NDVI data for 2019, 2020, and 2021. We then calculated the genetic effect for the intercept and 1st coefficients of the Legendre polynomials (L0 and L1) for each year.

We estimated the genetic correlations between CV and L0 and CV and L1, and estimated the genomic heritability of L0, L1, and CV. In seven out of nine combinations of treatments by years, negative genetic correlations between CV and L0 or CV and L1 were observed ( Table 3 ). L0 and L1 values for each genotype were associated with the intercepts and slopes of the time-series NDVI values. Therefore, these negative genetic correlations indicate that a low CV is associated with large intercepts and slopes of the time-series NDVI values. In 2019, there were no obvious patterns in genetic correlations or genomic heritability among the treatments ( Table 3 ). However, in 2020, treatments W10 and D10 showed higher genetic correlations and genomic heritability than Treatment D for all the traits. In 2021, the genetic correlation of L0 was -0.2 (p< 0.01) in Treatment D, whereas W10 and D10 showed higher genetic correlations of -0.4 (p< 0.001) and -0.51 (p< 0.001), respectively.

In Case1, we used CV, which was calculated from the fresh weight collected in nine combinations of treatments by years, as the same target trait and built MT RRM   model and MT All   model in each combination of treatments by years. The prediction accuracy of the simple genomic prediction model was 0.32. The prediction accuracy of MT RRM   model was higher than that of the simple genomic prediction model for all year treatments ( Figure 6A ). Three out of nine MT All   models evaluating plant’s stability under drought. We focused on the MT RRM   model . In 2019, there was no difference in the prediction accuracy among the treatments in MT RRM   model . However, the prediction accuracies of treatments W10 and D10 in 2020 were 25% and 30% higher, respectively than those of Treatment D in 2020. In addition, the prediction accuracies of treatments W10 and D10 in 2021 were higher by 14% and 11%, respectively than those of Treatment D in 2021. These results indicate that the time-series MS data collected during the treatment, which changed the irrigation pattern, were more useful than those in the continuous drought treatment for predicting the CV in 2020 and 2021. To evaluate the importance of using RRMs, we compared MT RRM   model and MT All   model . In the seven combinations of treatments by years, the prediction accuracies of MT RRM   model were higher than those of MT All   model ( Figure 6C ). In particular, the proportion of improvement was 14% for Treatment W10 in 2021.

To reduce the cost and labor required to create the prediction model, we changed the proportion of genotypes collected as secondary trait data in the training set. In 2020 and 2021, except for Treatment D in 2020, the prediction accuracy for all combinations of treatments by years increased as the proportion of genotypes with secondary trait data increased ( Figure 7 ). In Treatment D10 of 2020 and 2021, even when the proportions of data with secondary traits were 10%, the prediction accuracies of MT RRM   models were 0.40 and 0.42. These prediction accuracies were 23% and 30% higher than those of the simple genomic prediction model, respectively. Moreover, there was no obvious pattern in the proportion of improvement ( Figure S3 ).

In Case3, we evaluated the across-year predictions for each treatment. In all three treatments, MT RRM   models and MT All   models outperformed the simple genomic prediction model ( Figures 8A, B ). In cross-validations with the 2020 and 2021 datasets, the prediction accuracy of MT RRM   models for Treatments W10 and D10 were on average 42% and 42% higher than that of the simple genomic prediction model, respectively. The proportion of improvement in all MT RRM   models was greater than zero ( Figure 8C ). The proportion of improvement was, on average, 11, 6, and 6% for treatments W10, D10, and D, respectively. This result indicates that the time-series MS data should be modeled using RRM for across-year predictions.

In this study, we compared the prediction accuracy of a simple genomic prediction model and a genomic prediction model using time-series NDVI values as secondary traits to predict the CV ( Figures 6 – 8 ). In soybeans, a relationship between NDVI and drought tolerance has been reported (Zhou et al., 2020). Except for Case2 in 2019, all genomic prediction models using time-series NDVI values as a secondary trait showed higher prediction accuracy than a simple genomic prediction model ( Figures 6 – 8 ). These results suggest that time-series MS data are useful for predicting the CV.

In 2020 and 2021, the prediction accuracy of treatments W10 and D10, which changed irrigation treatments, was higher than that of Treatment D, which was the no watering treatment, in the same year ( Figures 6 , 7 ). A relationship between plant responses to changes in irrigation and drought tolerance has been reported. In soybeans, slow wilting resulted in drought tolerance (Fletcher et al., 2007; Pathan et al., 2014; Valliyodan et al., 2017; Ye et al., 2020) and a relationship between slow wilting and NDVI values under drought is observed (Zhou et al., 2020). Additionally, the speed of recovery from drought stress is associated with drought tolerance in several crop species (Hayano-Kanashiro et al., 2009; Kränzlein et al., 2021), including soybean (Hossain et al., 2014). However, no study has evaluated the plant’s stability under drought stress using time-series NDVI changes caused by changes in irrigation. This result indicates that plant responses to irrigation changes are useful for evaluating plant’s stability under drought stress.

In this study, the overall prediction accuracy of CV was lower in 2019 than that in the other years. There are two possible reasons for this result. One possible reason is that NDVI was measured in 2019 (four times) less frequently than in 2020 (six times) and 2021 (seven times). The small number of measurements may not capture the time-series changes in NDVI well. In addition, we collected NDVI values on the first day after the irrigation change in Treatment D10 in 2019 ( Figure 3 ). The second possible reason is the weather conditions. In 2020 and 2021, the soil moisture content increased during the irrigation period and decreased during the non-irrigation period ( Figure S1 ). In contrast, in 2019, treatments W5 and D10 showed the opposite trend of increasing soil moisture content during the last non-irrigation period. In 2019, it rained for 6 out of 10 days from 25 August to 3 September during the period of no irrigation in Treatment D10. The lack of drought stress during this period might have resulted in a failure to capture the plant response to irrigation change and low prediction accuracy in 2019.

In Case1, MTRRM model, which modeled time-series NDVI data and used L0 and L1 (genetic random regression coefficients) as secondary traits in the MTM, generally showed a higher prediction accuracy than MTAll model, which treated time-series NDVI data as independent traits and used each day’s NDVI value as a secondary trait in the MTM ( Figure 6 ). It was reported that the RRM is superior to the MTM in modeling time-series data (Kranis et al., 2007; Mota et al., 2013). In addition, there are problems with over-parameterization and high computational requirements when treating time-series data as independent traits and constructing MTMs (Misztal et al., 2000; Oh and See, 2008; Moreira et al., 2020). In terms of calculation efficiency, RRM is useful for modeling time-series data.

Campbell et al. (Campbell et al., 2018) built RRM using time-series projected shoot area (PSA) data for 179 rice lines obtained in their own study. Using the built-in RRM (Campbell et al., 2018), Campbell predicted time-series changes in PSA for new 178 rice lines collected in a different year (Campbell et al., 2017). This result suggests that RRM is useful for predictions across years. In our study, MTRRM model is superior to MTAll model in terms of prediction across years in Case3 ( Figure 8 ). MTRRM model can build a model using all-day data of NDVI value in each year, whereas MTAll model can only use dates that are common across all years. This difference in the amount of data may be the reason for the difference in prediction accuracy. When using MTAll model for predictions across years, it is necessary to match the data measurement dates. It was difficult to collect data on a desired date because of weather conditions or equipment problems. At this point, the RRM can be built without considering the measurement date and number of days of data measurement.

Environmental stability, also, the stability of phenotypes over environments, is an important target in breeding (Tollenaar and Lee, 2002; Sabaghnia et al., 2008; Mickelbart et al., 2015; Blum, 2018). The CV of phenotypic values over environments can be a good indicator for selecting stable genotypes. To evaluate genotypes’ stability across environments, it is necessary to conduct multi-environmental experiments which require high costs and labor (Fan et al., 2007; Karimizadeh et al., 2013; Akter et al., 2015; Lal et al., 2022). In this study, we proposed a method to measure the response to changes in irrigation by remote sensing and model it with RRM at a single location in multiple years. We found that genetic variation captured by the RRM was associated with plant’s stability under drought stress, even with only one location in three years. The results of Case1 validation suggest that genetic variation captured by measuring and modeling plant responses to irrigation in a single environment changes can predict the stability over nine environments ( Figure 6 ). The results of Case2 validation indicate that even when 10% of genotypes in a training set had secondary trait data, MTRRM model was more accurate than a simple genomic prediction model in treatments W10 and D10 in 2020 and 2021 ( Figure 7 ). Therefore, it is possible to reduce the number of genotypes measured for secondary traits in the training set. The results of Case3 validation also suggest that we can predict the plant’s stability across years ( Figure 8 ). Once a prediction model is constructed, the stability of unknown genotypes can be predicted based on time-series NDVI data. The results of Case2 and Case3 validations also demonstrate the flexibility of MTRRM model. When the stability of unknown genotypes under drought is evaluated in a single environment, it greatly reduces the time and cost, and thus streamlines breeding schemes. When calculating CVs, experimental designs with repetitions per treatment and year combination are usually used (Mohammadi and Amri, 2008; Das et al., 2010; Küchenmeister et al., 2012; Di Matteo et al., 2016). ¨ In this study, priority was given to increasing the number of genotypes over replications. This is because it is known that maximizing the number of genotypes (even with one replication) is the most effective when the number of plots is limited for QTL analysis and genomic prediction (Knapp and Bridges, 1990; Lorenz, 2013). Given the lack of replications and the small size of microplots, the fresh weight is probably influenced by non-genetic factors. However, the prediction accuracy of CV was improved by using information obtained from the methods for measuring and modeling changes in time-series NDVI data. Therefore, the effectiveness of the developed method was confirmed. In order to apply these measuring and modeling methods in breeding selection, it is necessary to associate the information obtained from the developed methods with drought tolerance-related traits and even yield-related traits in expanded field trials.

In this study, we developed a method to measure and model plant responses to irrigation changes and confirmed their association with plant’s stability under drought stress. Drone-based HTP allows us to capture time-series plant responses to environmental and irrigation changes (Chen et al., 2014; Dodig et al., 2021), and measuring and modeling plant responses to these environmental changes will provide important insights that have not yet been obtained previously (Arnold et al., 2019; Moreira et al., 2020). This developed method will contribute to the study of abiotic stress and genetic improvement of soybean.