A graphical summary of the procedures explained hereon is presented in Figure 1.

According to the joint phenotypic analysis, the plot level repeatabilities were 0.23 for GY and 0.59 for PH, revealing the traits as lowly and moderately heritable, respectively. The LRT test (p < 0.05) showed the effects of the environment, block within environment, genotype, and genotype (check) by environment to significantly affect both GY and PH. The BLUPs of GY averaged 6.79 Mg ha⁻¹ ranging from 4.90 to 8.36 Mg ha⁻¹. For PH, the mean was 199.02 cm, with values varying between 170.69 and 217.74 cm (Supplementary Figure S1). More information on these genotypes can be found in Alves et al. (2019), Fritsche-Neto et al. (2019), and Galli et al. (2020a).

The regression metrics for the prediction of PH and GY using GBLUP, MLP, or CNN are presented in Table 1. Considering each prediction scenario (e.g., PH with MLP), the metrics (MAE or MSE) were generally consistent across the inner training, inner validation, and outer validation sets. Therefore, scenario comparisons were performed, having the outer validation set as a reference. The values of MAE varied from 0.0635 to 0.0770 for PH and from 0.0797 to 0.0850 for GY. The loss function (MSE) varied from 0.0083 to 0.0109 for PH and from 0.0114 to 0.0128 for GY. The correlations varied from 0.68 to 0.75 for PH and from 0.53 to 0.59 for GY.

Patterns arose from comparisons between the studied scenarios. Contrasting the GBLUP (standard model) with MLP and CNN, GBLUP yielded superior results across all metrics for PH. For GY, a similar pattern was observed when comparing GBLUP with MLP, except for MSE. However, CNN outperformed GBLUP for this trait concerning all metrics. Regarding MLP or CNN, the latter presented better results for both traits considering all estimated metrics, except for r in PH.

The classification metrics for the prediction of PH and GY using MLP or CNN are presented in Table 2 and Supplementary Table S1. Under each prediction scenario, the loss increased from inner training to inner validation to outer validation. Regarding the other metrics (with few exceptions), values were (in average across scenarios) greater in inner training, followed by inner validation and outer validation sets. In order to facilitate the evaluation, comparisons between scenarios were performed using the outer validation set as a reference.

The observed values of metrics varied depending on the prediction scenario. The TNR varied from 0.26 to 0.65 for PH and 0.66 to 0.97 for GY; the recall ranged from 0.54 to 0.96 for PH and from 0.17 to 0.67 for GY; the precision presented values from 0.65 to 0.92 for PH and from 0.35 to 0.70 for GY; the F1score showed results from 0.57 to 0.93 for PH and from 0.22 to 0.67 for GY; the accuracy varied largely presenting values from 0.52 to 0.88 for PH and 0.66 to 0.90 for GY; the variation of balanced accuracy varied from 0.55 to 0.61 for PH and from 0.57 to 0.68 for GY; at last, the AUC ranged from 0.61 to 0.67 for PH and from 0.71 to 0.74 for GY.

The effect of selection intensity (moderate and extreme) presented tendencies to the estimated metrics. TNR, precision, recall, and F1 score were higher at extreme selection intensity for PH. For GY, the opposite was observed. The accuracy was higher at extreme selection intensity for both traits. The balanced accuracies using GRM were equal for both selection intensities for PH and GY. However, when CNN was used, this metric was lower at extreme selection intensity for both traits.

At last, regarding the AUC, moderate-intensity presented better values for both traits. Regarding the effect of the prediction method, for predicting PH, using MLP showed better values of precision, accuracy, balanced accuracy, and AUC. However, for recall and F1 score, this was only observed at extreme intensity. For GY, using MLP generally presented better results at extreme selection intensity, while CNN was superior at moderate selection intensity. The exceptions were precision, where image-based models were better for both intensities, and recall, which presented the opposite behavior.

The neural network structures that minimized the loss function for each replication under each scenario are presented on Supplementary File S1. The classification scenarios were constitutionally composed of an input layer as the first, a dense layer as the second-last summarizing all the neurons of the previous layer, and an activation layer with the sigmoid function to generate the output probabilities. Similarly, the regression scenarios had an input layer as the first and a dense layer as the last to summarize all neurons to one output. Nevertheless, the network structures were generally different, with few exceptions. Among these coincidences, seven out of nine were of the same task type (regression or classification), four were of the same trait (PH or GY), and three were of the same selection intensity (extreme or moderate). However, the number of parameters varied greatly (from 24,533 to 23,589,764), typically higher when images were used.

Dealing with GRM or images requires networks with specific internal layers. The scenarios with MLPs presented a varying number of dense layers (1 to 4); normalization layers (0 to 4; present in about half of the networks); ReLU activation function (positioned after dense layers except the last one); and dropout (0 to 4; present in about 2/3 of networks). The CNNs were composed of 2-dimensional convolutions (1 to 4 in classifications and 2 to 6 in regressions; present in all networks); normalization layers (0 or 1; present in about 2/3 of the networks); 2-dimensional max/global max or average pooling (0 to 3; present in about 2/3 of networks); dropout (0 to 3; present in about 2/3 of networks); image processing filters (resize, random flip, contrast, rotation, translation, and concatenation; 0 to 4; present in about half of the networks; being more common for PH). Also, ResNet50 and Xception networks appeared within 1/3 of the classification networks (more common for PH).

Finally, the preferred optimizer was Adam; Adadelta and SGD also appeared in a limited number of cases. The most common learning rate was 0.001, followed by 0.01, 0.00001, 0.0001, and 0.1. The dropout regularization had values of 0.5 (most common) and 0.25.

Benchmark studies suggest the inconsistent performance of neural networks compared to standard GP methods, which depends on a series of factors. We contrasted our findings to reference studies and explored how these factors might have affected the results in this context. Concerning the regression analysis, the GBLUP method outperformed MLP for both traits. One of the factors reported determining the best methodology is how modeling is performed. GP was carried out as a two-stage analysis. Hence, the genotypic value of hybrids across environments was obtained before prediction. Accordingly, the environmental source of variation was absent and could not be captured by the ML methods. For instance, it has been extensively shown that linear models (e.g., GBLUP or BMTME) tend to be outperformed by MLP in a multi-environmental joint analysis if the genotype by environment factor is not modeled for the prediction of PH and GY in maize. This holds under both single (Montesinos-López et al., 2018b) and multi-trait (Montesinos-López et al., 2018a) modeling contexts. Accordingly, MLP was outperformed by GBLUP for both traits in our study, supporting the suggested effect of modeling to the comparative outcome for the studied GP methods.

The use of CNN presented better results than GBLUP and MLP for predicting GY. This contrasts with the findings of Azodi et al. (2019), who suggested that Ridge regression BLUP, a GBLUP-equivalent method, outperforms both MLP and CNN for predicting several traits on numerous crops, including PH and GY in maize inbred lines. In this case, CNN was the poorest performing method for both traits. The inconsistency between the results of these studies regarding the performance of the CNNs for GY could be attributed to the restrictive search space for hyperparameters given the computational requirements for the analysis of an astonishing amount of studied traits and species by Azodi et al. (2019); we tailored ML models to each scenario within each trait, which is known to improve NN performance (Montesinos-López et al., 2018a). Another factor that might have led to this discrepancy was the use of pre-processed genomic information (genomic images) in our CNNs. They opted for using the raw genomic matrix. At last, inbred lines were used in their work, while hybrids were used in ours; studies suggested that CNNs tend to have better performance (than linear methods) when strong nonlinear (e.g., dominance) effects are present (Bellot et al., 2018; Abdollahi-Arpanahi et al., 2020); which is the case of GY in population we studied (Alves et al., 2019).

Regarding the underperformance of NN methods at predicting PH, tangible reasons could be pointed out. It has been hypothesized that the occurrence of extreme allelic frequencies (e.g., only two genotypes are present for a given locus) favors linear models by enabling the capture dominance and epistatic variance (Azodi et al., 2019); however, this does not hold for this dataset (Galli et al., 2020a). Also, PH is predominantly governed by additive allelic interactions (Alves et al., 2019), which enables linear models to capture a considerable proportion of the genotypic variance; nevertheless, regardless of the nature of the effects governing the traits under study, ML should always be (at least) as good as linear models given their ability to model linear relationships (Azodi et al., 2019), which was not the case. At last, a cause could be the number of training samples, which might not have been enough for modeling linear and nonlinear interactions between markers by the NN. This is a problem of common occurrence in plant breeding given the usually low number of samples, a large number of markers, and heterogeneity of data (Abdollahi-Arpanahi et al., 2020; Pook et al., 2020).

Overall, CNN presented better results compared to MLP. This advantage might have been due to the processing of the genomic matrix into the additive GRM, in the case of MLP. In this case, only the linear relationship between genotypes was modeled, which might compromise the potential of MLP to identify nonlinear effects. The genomic matrix could be used for further studies at the expense of computational time to overcome this issue.

Given their complex genetic nature, most plant traits present continuous phenotypes. Moreover, traits that were previously discretized by means of measurement ease, such as resistance to biotic stresses, have had their continuous nature better explored by high-throughput phenotyping (Galli et al., 2020b). Therefore, regression tasks are an adequate fit for genetic analysis, including GP. Nevertheless, plant breeding is globally a classification problem in which genotypes are assigned classes (Ornella et al., 2014), usually selected and non-selected. Hence, we elaborate on this problem, unifying ranking and selection by using classifying predictors. These prediction machines were evaluated using metrics that assess the model’s ability to distinguish which genotypes should be selected.

A critical step on classification tasks is the discretization of the continuous variable; when applicable. Discretizing traits has its inherent degree of subjectivity, regarding, e.g., the number of classes and which threshold values are used to classify the data. Accordingly, these choices have been reported to influence the performance of prediction models (Ornella et al., 2014; González-Camacho et al., 2016). Furthermore, a greater level of subjectivity is introduced when genotype classification as true positives or negatives before the prediction is based on the empirical distribution rather than the absolute value of the trait. Predicting genotypes from a related population, classification would not be tied to the percentiles of the distribution on the training population but to the genotypic values and genetic variants under each class and the genetic similarities across populations. In this context, the algorithm might be targeting, e.g., plants with a height between 1.95 and 2.05 m, but not the 10% or 50% best yielding hybrids since the distribution of genotypic values of a new population is likely to differ from the training population.

The classification problem was approached considering two scenarios: one highly imbalanced, where the size of the classes differed substantially (extreme SI), and one nearly balanced, where each class contained about half of the individuals (moderate SI). Both scenarios are plausible and of common occurrence in plant breeding, depending on the program stage. Nevertheless, imbalanced datasets should be evaluated with further cautiousness (Fernández et al., 2011). TNR, precision, recall, F1 score, accuracy, and AUC are examples of metrics sensitive to class imbalance, meaning their results might not be directly interpretable for comparing predictions with differing selection intensities. This is also evidenced by the discrepancy between the accuracy and the balanced accuracy at extreme selection intensity. The selection intensities presented little influence over the balanced accuracy for the same trait and independent variable, except for GY when images were used.

The balanced accuracy is calculated by averaging the proportion of correct predictions in each class, meaning that the label (selected or non-selected) is not relevant. This metric varied from 0.55 to 0.61 for PH and 0.57 to 0.68 for GY. These results are inconsistent with the regression analysis, which showed higher predictability for PH according to all metrics, following the higher heritability of this trait. We postulate that this is associated with the region of the empirical density of genotypic values from which genotypes were regarded as “selected”. For PH, the distinction between the best and the worst individuals was non-directional, which might have difficulted the distinction between which hybrids should or not be selected by the models. For GY, as the selection is directional, this was not an issue. Overall, balanced accuracies were closer to 0.5 (random guess) than to 1 (all correct) for both traits, meaning that further improvements are required. Nevertheless, our results suggest the possibility of non-directional selection, as for PH, which is highly relevant for breeding programs.

Unlike the regression task, where the use of CNN usually presented the best results between machine learning methods, there was considerable inconsistency regarding the superiority of MLP or CNN in the classification task. The comparative performance of the neural network methodologies seemed highly conditioned to trait and selection intensity. Generally, MLP presented the best results for PH, while for GY, the best method heavily depended on the selection intensity. Therefore, it is reasonable to assume that the discretization process of PH and GY impacted the performance of CNN more than that of MLP; but further investigation is warranted. GP prediction regression tasks with machine learning models are already common, but studies comparing methods for predicting discretized variables are still limited.

The choice of the neural network hyperparameters has been a critical step for NN-based GP by extensive benchmarking (Azodi et al., 2019). Ergo, network search for a given task and dataset has been applied in recent ML-based GP studies (Montesinos-López et al., 2018b; Azodi et al., 2019; Abdollahi-Arpanahi et al., 2020). However, model tuning has been primarily performed using naïve approaches such as random (values sampled from distribution) or grid (discrete values) search, which may limit the number of hyperparameter combinations based on a set of user-defined a priori information (Jin et al., 2019). Due to recent advancements in computer science and technology, less restrictive, free, and easy-to-use hyperparameter search algorithms have been made available. Hence, we used Auto-Keras, an AutoML search algorithm with Bayesian optimization to identify (suitable) models. Overall, the algorithm yielded adequate performing neural networks despite the absence of the commonly required human intervention for adjustments.

It has been previously reported that different neural network hyperparameters can be obtained from network search algorithms for a given task (Bellot et al., 2018; Huang et al., 2020). The neural networks selected by the AutoML algorithm presented idiosyncrasies within replications of the same scenario (File S1). The lack of similarity between structures might arise from the ability of AutoML to adapt the network to the dataset (Jin et al., 2019), which changes due to sampling in repeated validation. Huang et al. (2020) suggest this event to be a consequence of insufficient data, but further confirmation is required. Additionally, this may also be associated with the sampling nature of the hyperparameter search system (Jin et al., 2019). Despite the inconsistency between structures, systematic regularities are suggested by the within scenario low standard deviation of the estimated metrics (Tables 1 and 2). Hence, the networks might be capturing similar features, yielding consistent predictions. This has relevant implications for the choice of (deep) neural networks, meaning that distinct but adequate network structures result in similar outcomes.

Further observations can be drawn from the chosen network structures: (i) Although limited, the cases where structures did match (within and between scenarios) suggest that: type of task (regression or classification) is determinant over structure since most matches were of the same type; matches across traits were common, suggesting that similar sources of information might have been captured, which is probably intrinsically associated to the genetic correlation between PH and GY in maize. (ii) Also, when images were used as the input for prediction, augmentation procedures (e.g., resize, flip, rotation) were allocated in the structure of about half of the chosen models despite the spatial structure in the genomic images created by the decomposition performed by DeepInsight; further inferences on this matter would require studying the implications of such procedures to the original images, which is not in the scope of this study. (iii) At last, the depth and number of parameters of the networks within scenarios were highly variable for both MLP and CNN, suggesting that simple architectures were as effective as the more complex ones. Simpler models also have the advantage of being generally quicker to train (van Dijk et al., 2021). (iv) Regarding overfitting, which is the tendency of a model to perform well on training but not on unseen data (van Dijk et al., 2021), some differences in performance could be observed between inner training, inner validation, and outer validation sets, but further investigation would be required to determine their extent and consequences. The dropout regularization, temporarily setting a percentage of random neurons to zero (Srivastava et al., 2014), was present on 2/3 of the chosen models, presumably acting on the overfitting issue (Montesinos-López et al., 2018b).

Overall, based on the empirical and experimental evidence, neural networks are especially competitive under the presence of strong nonlinear factors and interactions and hidden relationships between pieces of information. Accordingly, it is also dependent on the population type (e.g., lines or hybrids) and the consequent, non-mutually exclusive, genetic architecture of the trait (Bellot et al., 2018; Abdollahi-Arpanahi et al., 2020). The performance of NN is certainly conditioned on the choice of hyperparameters (Bellot et al., 2018; Zingaretti et al., 2020) and neural network type (MLP or CNN). It depends on how the input data is processed before prediction; consequently, special attention should be given to this step since valuable information could be lost. Also, it is presumably dependent on the number of samples and the sample to parameter ratio (Montesinos-López et al., 2018a; Azodi et al., 2019; Pérez-Enciso and Zingaretti, 2019; Abdollahi-Arpanahi et al., 2020). Therefore, it is the scientist’s discretion to test and identify the best performing method for their task. To this day, the only identified consistency regarding GP benchmarking is that no model performs best for all situations.

From experience, inferences on using images for GP could be drawn. In the original work by Sharma et al. (2019), DeepInsight was used for transforming RNA-seq, text, and artificial datasets into images. Our work is the first to apply such methodology in a GP context, and it is noteworthy that: (1) the algorithm can create images of different sizes. Image size, which is a hyperparameter, should be adapted to the available dataset and computational power. With the increasing size of the genomic matrix, there is a greater chance that a considerable amount of information would be lost as correlated markers would be tightly grouped, so larger images should be used (Sharma et al., 2019). Additionally, increasing the size of images consequently increases the number of parameters estimated in the neural network, requiring greater computational power. In this work, using 120 by 120 images seemed to be an adequate fit for ~30,000 genomic markers; (2) different dimensionality reduction techniques can be used: t-SNE and kPCA are implemented in the algorithm, but any other of interest can be implemented; further testing should elaborate on this matter; (3) images can have multiple layers: neural networks can model linear and nonlinear relationships between neurons, including other effects and layers of data, such as dominance, epistasis, g × e, transcriptome, and so on.; (4) the cost–benefit in terms of predictive gain and additional work, the use of images as input is arguable. Nevertheless, the methodology’s potential for GP is unprecedented; (5) simulations should provide new valuable and unbiased information.

At last, we discussed two prediction alternatives: regression and classification. Under the regression context, MLP and CNN presented competitive results. Under the classification context, we expected better performances. Nevertheless, we believe that the latter has great potential for plant breeding since it simplifies the pipeline. Neural networks are self-adaptable and aimed at prediction alone. This statement implies that understanding and exposing the events underlying the relationship between phenotypes and genotypes are not of particular interest but could be done if necessary (Azodi et al., 2020). This also implies that limited genetic knowledge of the trait is not a constraint for prediction. Coupled with a simpler processing, direct classification opens new possibilities regarding selecting traits where the ideotype points to intermediate phenotypes, e.g., plant height, ear height, and flowering time (under some circumstances) in maize. Hence, we believe this methodology deserves attention since it could further enhance the GP pipeline in breeding programs.