# Enhancing adaptation of tropical maize to temperate environments using genomic selection

Nicole Choquette, Teclemariam Weldekidan, Jason Brewer, Scott Davis, Randall J. Wisser, James Holland

► To cite this version:

Nicole Choquette, Teclemariam Weldekidan, Jason Brewer, Scott Davis, Randall J. Wisser, et al. Enhancing adaptation of tropical maize to temperate environments using genomic selection. G3, 2023, 13 (9), 10.1093/g3journal/jkad141. hal-04264611

HAL Id: hal-04264611

https://hal.inrae.fr/hal-04264611

Submitted on 30 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Enhancing adaptation of tropical maize to temperate environments using genomic selection

Nicole Choquette, Teclemariam Weldekidan, Jason Brewer, Scott Davis,

Randall J. Wisser, James Holland

► To cite this version:

Nicole Choquette, Teclemariam Weldekidan, Jason Brewer, Scott Davis, Randall J. Wisser, et al. Enhancing adaptation of tropical maize to temperate environments using genomic selection. G3, 2023, 13 (9), 10.1093/g3journal/jkad141. hal-04264611

HAL Id: hal-04264611

https://hal.inrae.fr/hal-04264611v1

Submitted on 30 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
<figure><img src="image_1.png" /></figure>

# Enhancing adaptation of tropical maize to temperate environments using genomic selection

Nicole E. Choquette,<fnref n="1" /> Teclemariam Weldekidan,<fnref n="2" /> Jason Brewer,<fnref n="3" /> Scott B. Davis,<fnref n="2" /> Randall J. Wisser,<fnref n="2" />,<fnref n="4" />

James B. Holland <fnref n="1" />,<fnref n="3" />,*

## Abstract

Tropical maize can be used to diversify the genetic base of temperate germplasm and help create climate-adapted cultivars. However, tropical maize is unadapted to temperate environments, in which sensitivities to long photoperiods and cooler temperatures result in severely delayed flowering times, developmental defects, and little to no yield. Overcoming this maladaptive syndrome can require a decade of phenotypic selection in a targeted, temperate environment. To accelerate the incorporation of tropical diversity in temperate breeding pools, we tested if an additional generation of genomic selection can be used in an off-season nursery where phenotypic selection is not very effective. Prediction models were trained using flowering time recorded on random individuals in separate lineages of a heterogenous population grown at two northern U.S. latitudes. Direct phenotypic selection and genomic prediction model training was performed within each target environment and lineage, followed by genomic prediction of random intermated progenies in the off-season nursery. Performance of genomic prediction models was evaluated on self-fertilized progenies of prediction candidates grown in both target locations in the following summer season. Prediction abilities ranged from 0.30 to 0.40 among populations and evaluation environments. Prediction models with varying marker effect distributions or spatial field effects had similar accuracies. Our results suggest that genomic selection in a single off-season generation could increase genetic gains for flowering time by more than 50% compared to direct selection in summer seasons only, reducing the time required to change the population mean to an acceptably adapted flowering time by about one-third to one-half.

**Keywords:** exotic germplasm, genomic prediction, GenPred, quantitative genetics, flowering time

## Introduction

Maize yields have increased consistently since the 1940s due to improved genetics, agronomic management intensification, and a synergistic interaction between modern hybrids and agronomy (Duvick 2005). Globally, maize production primarily relies on germplasm derived from temperate U.S. germplasm sources (Smith et al. 2022). Hybrid maize cultivars in the U.S. contain relatively little of the total genetic diversity found in maize globally (Liu et al. 2003; Pollak 2003; Mikel and Dudley 2006; Romay et al. 2013). This stems from the initial founding of North American breeding pools from only a few pre-adapted temperate populations followed by decades of intense selection in breeding programs (Hufford et al. 2012). In commercial programs, high-performing germplasm is continuously reused with little incorporation of outside or exotic germplasm (Mikel and Dudley 2006; White et al. 2020). Using diverse genetic resources to introduce novel alleles into modern temperate breeding pools can contribute to better disease resistance and agronomic performance, and potentially to better climate change resilience (Holland et al. 1996; Pollak 2003; Nelson and Goodman 2008; Cortés and

López-Hernández 2021). Facing the pressure of a growing global population and a changing climate with emerging pests and diseases, new strategies to capitalize on exotic diversity may help to maintain genetic gain and minimize genetic vulnerability (Rogers et al. 2022).

Among exotic germplasm sources, tropical maize is considered the most useful for integrating into temperate breeding programs (Goodman 2004). Tropical maize harbors greater allelic richness than temperate maize (Liu et al. 2003; Romay et al. 2013). Both private and public maize breeding programs recognize the potential of tropical maize. For example, exotic alleles for resistance to stalk rot derived from tropical germplasm have been deployed in elite commercial cultivars (Frey et al. 2011). Additionally, the United States Department of Agriculture collaborates with private companies and public sector researchers to integrate tropical germplasm into elite commercial genetic backgrounds through the Germplasm Enhancement of Maize (GEM) program (Rogers et al. 2022).

Tropical maize is adapted to shorter photoperiods and warmer nighttime temperatures than those that occur in the growing season of temperate U.S. production environments. Flowering time in

<footnote n="1">¹Department of Crop and Soil Sciences, North Carolina State University, Raleigh, NC 27695, USA</footnote>

<footnote n="2">²Department of Plant and Soil Sciences, University of Delaware, Newark, DE 19716, USA</footnote>

<footnote n="3">³USDA-ARS Plant Science Research Unit, Raleigh, NC 27695, USA</footnote>

<footnote n="4">⁴Laboratoire d'Écophysiologie des Plantes sous Stress Environmentaux, INRAE, University of Montpellier, L'Institut Agro, Montpellier, FR 34000, USA</footnote>

<footnote>*Corresponding author: USDA-ARS Plant Science Research Unit, Raleigh, NC 27695, USA. Email: jim.holland@usda.gov</footnote>
<figure><img src="image_1.png" /></figure>

# Enhancing adaptation of tropical maize to temperate environments using genomic selection

Nicole E. Choquette,<fnref n="1" /> Teclemariam Weldekidan,<fnref n="2" /> Jason Brewer,<fnref n="3" /> Scott B. Davis,<fnref n="2" /> Randall J. Wisser,<fnref n="2" />,<fnref n="4" />

James B. Holland <fnref n="1" />,<fnref n="3" />,*

## Abstract

Tropical maize can be used to diversify the genetic base of temperate germplasm and help create climate-adapted cultivars. However, tropical maize is unadapted to temperate environments, in which sensitivities to long photoperiods and cooler temperatures result in severely delayed flowering times, developmental defects, and little to no yield. Overcoming this maladaptive syndrome can require a decade of phenotypic selection in a targeted, temperate environment. To accelerate the incorporation of tropical diversity in temperate breeding pools, we tested if an additional generation of genomic selection can be used in an off-season nursery where phenotypic selection is not very effective. Prediction models were trained using flowering time recorded on random individuals in separate lineages of a heterogenous population grown at two northern U.S. latitudes. Direct phenotypic selection and genomic prediction model training was performed within each target environment and lineage, followed by genomic prediction of random intermated progenies in the off-season nursery. Performance of genomic prediction models was evaluated on self-fertilized progenies of prediction candidates grown in both target locations in the following summer season. Prediction abilities ranged from 0.30 to 0.40 among populations and evaluation environments. Prediction models with varying marker effect distributions or spatial field effects had similar accuracies. Our results suggest that genomic selection in a single off-season generation could increase genetic gains for flowering time by more than 50% compared to direct selection in summer seasons only, reducing the time required to change the population mean to an acceptably adapted flowering time by about one-third to one-half.

**Keywords:** exotic germplasm, genomic prediction, GenPred, quantitative genetics, flowering time

## Introduction

Maize yields have increased consistently since the 1940s due to improved genetics, agronomic management intensification, and a synergistic interaction between modern hybrids and agronomy (Duvick 2005). Globally, maize production primarily relies on germplasm derived from temperate U.S. germplasm sources (Smith et al. 2022). Hybrid maize cultivars in the U.S. contain relatively little of the total genetic diversity found in maize globally (Liu et al. 2003; Pollak 2003; Mikel and Dudley 2006; Romay et al. 2013). This stems from the initial founding of North American breeding pools from only a few pre-adapted temperate populations followed by decades of intense selection in breeding programs (Hufford et al. 2012). In commercial programs, high-performing germplasm is continuously reused with little incorporation of outside or exotic germplasm (Mikel and Dudley 2006; White et al. 2020). Using diverse genetic resources to introduce novel alleles into modern temperate breeding pools can contribute to better disease resistance and agronomic performance, and potentially to better climate change resilience (Holland et al. 1996; Pollak 2003; Nelson and Goodman 2008; Cortés and

López-Hernández 2021). Facing the pressure of a growing global population and a changing climate with emerging pests and diseases, new strategies to capitalize on exotic diversity may help to maintain genetic gain and minimize genetic vulnerability (Rogers et al. 2022).

Among exotic germplasm sources, tropical maize is considered the most useful for integrating into temperate breeding programs (Goodman 2004). Tropical maize harbors greater allelic richness than temperate maize (Liu et al. 2003; Romay et al. 2013). Both private and public maize breeding programs recognize the potential of tropical maize. For example, exotic alleles for resistance to stalk rot derived from tropical germplasm have been deployed in elite commercial cultivars (Frey et al. 2011). Additionally, the United States Department of Agriculture collaborates with private companies and public sector researchers to integrate tropical germplasm into elite commercial genetic backgrounds through the Germplasm Enhancement of Maize (GEM) program (Rogers et al. 2022).

Tropical maize is adapted to shorter photoperiods and warmer nighttime temperatures than those that occur in the growing season of temperate U.S. production environments. Flowering time in

<footnote n="1">¹Department of Crop and Soil Sciences, North Carolina State University, Raleigh, NC 27695, USA</footnote>

<footnote n="2">²Department of Plant and Soil Sciences, University of Delaware, Newark, DE 19716, USA</footnote>

<footnote n="3">³USDA-ARS Plant Science Research Unit, Raleigh, NC 27695, USA</footnote>

<footnote n="4">⁴Laboratoire d'Écophysiologie des Plantes sous Stress Environmentaux, INRAE, University of Montpellier, L'Institut Agro, Montpellier, FR 34000, USA</footnote>

<footnote>*Corresponding author: USDA-ARS Plant Science Research Unit, Raleigh, NC 27695, USA. Email: jim.holland@usda.gov</footnote>
most tropical maize is quite sensitive to photoperiod, and correct timing of flowering is crucial for adaptation and fitness, as it controls the switch from vegetative growth to reproduction (Gouesnard et al. 2002). Therefore, tropical maize varieties tend to be maladapted to U.S. temperate environments—they flower too late in the season to achieve acceptable yields (Teixeira et al. 2015; Choquette et al. 2023). Phenological adaptation of tropical populations can be achieved through phenotypic recurrent selection, which also results in indirect improvements to seed yield and other agronomic traits (Choquette et al. 2023). For populations with low linkage disequilibrium and using selection procedures that minimize genetic drift, most of the standing variation can be maintained while the population is phenologically adapted (Wisser et al. 2019). However, phenotypic recurrent selection may require a decade to change the mean flowering time of a tropical population to an adapted state (Teixeira et al. 2015; Wisser et al. 2019). New approaches that reduce this timeline could facilitate the diversification of temperate maize.

Genomic selection offers the potential to reduce the time for adaptational breeding. Genomic selection leverages available genomic and phenotypic information on a training set of breeding families to predict breeding values for individuals with only genomic information (Meuwissen et al. 2001; Lorenz et al. 2011). Importantly, genomic selection can increase rates of genetic gain per unit time by allowing selection to be practiced in off-target environments (Heffner et al. 2009). However, accuracy of genomic selection depends on trait heritability, training population size, number of markers, linkage disequilibrium, and the relationship between the training and testing population sets (Heffner et al. 2009; Lorenzana and Bernardo 2009; Jannink et al. 2010; Crossa et al. 2017). As genomic technologies improve, the relative cost of genotyping to phenotyping continues to decrease, making genomic selection a more feasible option for breeders to achieve greater genetic gain per year.

Genomic selection could be exploited for adaptation of tropical maize to temperate environments. Maize flowering time has high heritability and a complex genetic architecture comprised of mostly additive effects (Chardon et al. 2004; Salvi et al. 2007; Buckler et al. 2009), conditions that are ideal for genomic selection. Previous research has indicated that phenotypic selection for early flowering time in low-latitude environments, which are employed by temperate breeding programs as "off-season" (winter) nurseries, show a limited response to selection and do not relieve other maladaptation effects (Choquette et al. 2023). Off-season nurseries typically have short daylengths that mask the expression of photoperiod-sensitivity genes, in addition to warm temperatures that compress differences in flowering time between individuals. Together, this makes off-season selection for earliness ineffective for tropical maize being adapted to the northern temperate zone. Here, we propose that genomic selection trained in target temperate environments can be used for selection in off season nurseries. This would provide an additional generation of selection per year compared to only selecting each summer in the target environment. Alternating generations of phenotypic selection in the target environment and genomic selection in the off-season nursery could lead to greater genetic gain per year.

In this study, we used two tropical synthetic maize populations, each created from two previous generations of phenotypic selection, to evaluate the use of genomic selection for earlier flowering in different temperate locations. We also compared genomic prediction models encompassing a range of genetic architectures (purely polygenic or including larger-effect variants; purely additive or additive and dominance variances) and non-genetic field

spatial effects. Results from this study can guide the use of genomic selection for geographical adaptation in maize and other crops.

## Materials and methods

### Phenotypic data and selection scheme

This study used populations derived from a tropical synthetic (TropicS) population created from seven tropical inbred lines (Weldekidan et al. 2022) that have good yield combining ability with temperate testers (Nelson and Goodman 2008). A common base population of the TropicS was derived by multiple generations of intercrossing and random intermating among the parental lines. In a separate study, this was followed by phenotypic selection for early flowering for two generations independently at eight locations (Weldekidan et al. 2022; Fig. 1). Here, we started with remnant non-inbred ($S_0$ generation) seeds of the second generation (G2) populations selected for early flowering at Delaware (TropicS G2-DE; latitude: 39.67 °N) and North Carolina (TropicS G2-NC; latitude: 35.67 °N). In the summer of 2019, 5,000 seeds of the TropicS-G2-DE $S_0$ population were planted in Newark, DE and 5,000 seeds of the TropicS-G2-NC population were planted in Clayton, NC. In both locations, plant stands were established with approximately 25.4 cm spacing between plants, and each experimental field was partitioned into four equally sized blocks. Within each block, 250 random plants were tagged to track their individual plant phenotypes for days to anthesis (DTA), days to silking (DTS), plant height, and ear height. To avoid edge effects on plant phenotypes, plants at the borders of each field were excluded from sampling. The 1,000 plants within each environment constituted location-specific training populations for genomic prediction. Leaf samples from each plant were collected for genotyping.

In both locations, phenotypic mass selection for early flowering time was also performed by intermating the 256 earliest-flowering plants (32 pollinations per block; half of which were among the earliest to shed pollen, and half among the earliest to produce silks). This 5% proportion of selected plants corresponds to a selection intensity of about 2.0 standard deviations (Falconer and Mackay 1996). A modified pollination scheme used by Weldekidan et al. (2022) was used to intermate selected plants. Briefly, a modified chain sib pollination method was used wherein pollen from five selected males was bulked together to pollinate five selected females, using each plant only once as a male or female parent, without selfing. To generate G3, which would be used for genomic selection, balanced seed bulks were created by combining 30 kernels from each ear of the selected plants.

In an off-season nursery in Homestead, FL (latitude: 25.50 °N), 3,150 seeds of the TropicS-G3-DE and TropicS-G3-NC populations were planted in October in 2019. Plant tissue was collected for genotyping from 768 random individuals in G3-NC and 1,536 individuals in G3-DE. All genotyped individuals were self-pollinated to create $S_{0.1}$ families and their dates of pollination were recorded.

### Genotyping

All plants were genotyped using a genotyping-by-sequencing (GBS) protocol optimized for heterozygous plants (Manching et al. 2017). The sequence data was processed using the RedRep pipeline (https://github.com/UD-CBCB/RedRep) with version 4 of the B73 reference genome. RedRep uses GATK genomicsDB for variant calling, which was applied to over 4,000 samples, including 3,682 samples designated for the current study. To reduce processing time, variant calling was run in parallel for 10 Mb windows
most tropical maize is quite sensitive to photoperiod, and correct timing of flowering is crucial for adaptation and fitness, as it controls the switch from vegetative growth to reproduction (Gouesnard et al. 2002). Therefore, tropical maize varieties tend to be maladapted to U.S. temperate environments—they flower too late in the season to achieve acceptable yields (Teixeira et al. 2015; Choquette et al. 2023). Phenological adaptation of tropical populations can be achieved through phenotypic recurrent selection, which also results in indirect improvements to seed yield and other agronomic traits (Choquette et al. 2023). For populations with low linkage disequilibrium and using selection procedures that minimize genetic drift, most of the standing variation can be maintained while the population is phenologically adapted (Wisser et al. 2019). However, phenotypic recurrent selection may require a decade to change the mean flowering time of a tropical population to an adapted state (Teixeira et al. 2015; Wisser et al. 2019). New approaches that reduce this timeline could facilitate the diversification of temperate maize.

Genomic selection offers the potential to reduce the time for adaptational breeding. Genomic selection leverages available genomic and phenotypic information on a training set of breeding families to predict breeding values for individuals with only genomic information (Meuwissen et al. 2001; Lorenz et al. 2011). Importantly, genomic selection can increase rates of genetic gain per unit time by allowing selection to be practiced in off-target environments (Heffner et al. 2009). However, accuracy of genomic selection depends on trait heritability, training population size, number of markers, linkage disequilibrium, and the relationship between the training and testing population sets (Heffner et al. 2009; Lorenzana and Bernardo 2009; Jannink et al. 2010; Crossa et al. 2017). As genomic technologies improve, the relative cost of genotyping to phenotyping continues to decrease, making genomic selection a more feasible option for breeders to achieve greater genetic gain per year.

Genomic selection could be exploited for adaptation of tropical maize to temperate environments. Maize flowering time has high heritability and a complex genetic architecture comprised of mostly additive effects (Chardon et al. 2004; Salvi et al. 2007; Buckler et al. 2009), conditions that are ideal for genomic selection. Previous research has indicated that phenotypic selection for early flowering time in low-latitude environments, which are employed by temperate breeding programs as “off-season” (winter) nurseries, show a limited response to selection and do not relieve other maladaptation effects (Choquette et al. 2023). Off-season nurseries typically have short daylengths that mask the expression of photoperiod-sensitivity genes, in addition to warm temperatures that compress differences in flowering time between individuals. Together, this makes off-season selection for earliness ineffective for tropical maize being adapted to the northern temperate zone. Here, we propose that genomic selection trained in target temperate environments can be used for selection in off season nurseries. This would provide an additional generation of selection per year compared to only selecting each summer in the target environment. Alternating generations of phenotypic selection in the target environment and genomic selection in the off-season nursery could lead to greater genetic gain per year.

In this study, we used two tropical synthetic maize populations, each created from two previous generations of phenotypic selection, to evaluate the use of genomic selection for earlier flowering in different temperate locations. We also compared genomic prediction models encompassing a range of genetic architectures (purely polygenic or including larger-effect variants; purely additive or additive and dominance variances) and non-genetic field

spatial effects. Results from this study can guide the use of genomic selection for geographical adaptation in maize and other crops.

## Materials and methods

### Phenotypic data and selection scheme

This study used populations derived from a tropical synthetic (TropicS) population created from seven tropical inbred lines (Weldekidan et al. 2022) that have good yield combining ability with temperate testers (Nelson and Goodman 2008). A common base population of the TropicS was derived by multiple generations of intercrossing and random intermating among the parental lines. In a separate study, this was followed by phenotypic selection for early flowering for two generations independently at eight locations (Weldekidan et al. 2022; Fig. 1). Here, we started with remnant non-inbred ($S_0$ generation) seeds of the second generation (G2) populations selected for early flowering at Delaware (TropicS G2-DE; latitude: 39.67 °N) and North Carolina (TropicS G2-NC; latitude: 35.67 °N). In the summer of 2019, 5,000 seeds of the TropicS-G2-DE $S_0$ population were planted in Newark, DE and 5,000 seeds of the TropicS-G2-NC population were planted in Clayton, NC. In both locations, plant stands were established with approximately 25.4 cm spacing between plants, and each experimental field was partitioned into four equally sized blocks. Within each block, 250 random plants were tagged to track their individual plant phenotypes for days to anthesis (DTA), days to silking (DTS), plant height, and ear height. To avoid edge effects on plant phenotypes, plants at the borders of each field were excluded from sampling. The 1,000 plants within each environment constituted location-specific training populations for genomic prediction. Leaf samples from each plant were collected for genotyping.

In both locations, phenotypic mass selection for early flowering time was also performed by intermating the 256 earliest-flowering plants (32 pollinations per block; half of which were among the earliest to shed pollen, and half among the earliest to produce silks). This 5% proportion of selected plants corresponds to a selection intensity of about 2.0 standard deviations (Falconer and Mackay 1996). A modified pollination scheme used by Weldekidan et al. (2022) was used to intermate selected plants. Briefly, a modified chain sib pollination method was used wherein pollen from five selected males was bulked together to pollinate five selected females, using each plant only once as a male or female parent, without selfing. To generate G3, which would be used for genomic selection, balanced seed bulks were created by combining 30 kernels from each ear of the selected plants.

In an off-season nursery in Homestead, FL (latitude: 25.50 °N), 3,150 seeds of the TropicS-G3-DE and TropicS-G3-NC populations were planted in October in 2019. Plant tissue was collected for genotyping from 768 random individuals in G3-NC and 1,536 individuals in G3-DE. All genotyped individuals were self-pollinated to create $S_{0.1}$ families and their dates of pollination were recorded.

### Genotyping

All plants were genotyped using a genotyping-by-sequencing (GBS) protocol optimized for heterozygous plants (Manching et al. 2017). The sequence data was processed using the RedRep pipeline (https://github.com/UD-CBCB/RedRep) with version 4 of the B73 reference genome. RedRep uses GATK genomicsDB for variant calling, which was applied to over 4,000 samples, including 3,682 samples designated for the current study. To reduce processing time, variant calling was run in parallel for 10 Mb windows
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Outline of the selection study. Phenotypic selection and genomic prediction model training were conducted in the target environment, genomic predictions were applied in the off-season nursery, and genomic predictions were evaluated on progeny lines in the target environment.</figcaption></figure>

of the maize genome, which were combined into a single vcf file for downstream analysis. Initially, several technical QC steps were sequentially applied: (i) genotype scores based on <12X sequence depth were set to NA; (ii) individual samples with >97.5% missing data were filtered; (iii) markers with greater than two alleles were filtered; and (iv) markers with a minor allele frequency (MAF) <0.4% were filtered. Because nearly all samples were derived from the same set of parents, the low MAF filter was intended to remove false variant calls. After technical QC filtering, 9,554 bi-allelic markers (SNPs + Indels) were retained and missing genotype calls were imputed with MaCH (Li et al. 2010). Afterward, a more stringent missing data filter was applied to remove markers imputed from <20% raw data, leaving 5,718 markers. Additionally, within each population, markers that deviated from Hardy-Weinberg expectations with P-values less than 0.05/m, where m is the number of markers, were excluded. After these quality control measures were performed, linkage disequilibrium was measured within populations as the squared correlation for all pairs of SNP loci (r²) using TASSEL (Bradbury et al. 2007). Finally, the intersection of remaining markers (4,340 SNPs) across populations was used for all subsequent analyses.

An additive genomic relationship matrix (G) was calculated using the R package AGHmatrix (Amadeu et al. 2016) using the VanRaden parameterization (VanRaden 2008; Isik et al. 2017):

$$ \frac{ZZ'}{2\Sigma p_i(1 - p_i)} $$

where ZZ' is the marker-based covariance matrix that reflects identity-by-descent relationships among the genotyped individuals and pᵢ is the minor allele frequency at each marker.

Genomic inbreeding coefficients for each individual were estimated as Gᵢᵢ - 1, where Gᵢᵢ is the diagonal element of G

corresponding to individual i. Individuals with inbreeding coefficients less than -0.2 or greater than 0.2 were removed. Additionally, closely related individuals within each population (Gᵢᵢ > 0.6) as well as individuals with suspected genotype sampling errors were identified and removed. After quality control based on inbreeding and pairwise relatedness, 665 G2-DE, 1,201 G3-DE, 698 G2-NC, and 613 G3-NC individuals remained. The additive genomic relationship matrix (G) was recalculated for this subset of lines. In addition, the matrix of realized dominance genomic relationships (D) was also calculated using the AGHmatrix package (Vitezica et al. 2013; Amadeu et al. 2016). Genetic relationship matrices were calculated for each population (DE or NC) separately as well as for both populations combined.

The genomic fixation index (Fst) was estimated between all pairs of the G2-DE and G3-DE and G2-NC and G3-NC populations using the R package hierstat (Goudet and Jombart 2015). Additionally, the expected heterozygosity was calculated for all of the populations using the adegenet R package (Jombart 2008).

## Initial genomic prediction to select G3 individuals for evaluation

The effectiveness of genomic selection in the off-season nursery was assessed using a subset of S₀.₁ families of G3 predicted plants. These were evaluated in 2020 at the original locations of phenotypic selection (target environments). Field space constraints limited the total number of progenies that could be evaluated in summer 2020. Therefore, we chose progenies from both the earliest and latest tails of genomic predictions as well as randomly sampled individuals for evaluation. Time constraints between processing the genotypic data (extracting DNA, sequencing, calling SNPs, quality control the data, and predicting GBLUPs) and planting time for the summer of 2020 limited
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Outline of the selection study. Phenotypic selection and genomic prediction model training were conducted in the target environment, genomic predictions were applied in the off-season nursery, and genomic predictions were evaluated on progeny lines in the target environment.</figcaption></figure>

of the maize genome, which were combined into a single vcf file for downstream analysis. Initially, several technical QC steps were sequentially applied: (i) genotype scores based on <12X sequence depth were set to NA; (ii) individual samples with >97.5% missing data were filtered; (iii) markers with greater than two alleles were filtered; and (iv) markers with a minor allele frequency (MAF) <0.4% were filtered. Because nearly all samples were derived from the same set of parents, the low MAF filter was intended to remove false variant calls. After technical QC filtering, 9,554 bi-allelic markers (SNPs + Indels) were retained and missing genotype calls were imputed with MaCH (Li et al. 2010). Afterward, a more stringent missing data filter was applied to remove markers imputed from <20% raw data, leaving 5,718 markers. Additionally, within each population, markers that deviated from Hardy-Weinberg expectations with P-values less than 0.05/m, where m is the number of markers, were excluded. After these quality control measures were performed, linkage disequilibrium was measured within populations as the squared correlation for all pairs of SNP loci (r²) using TASSEL (Bradbury et al. 2007). Finally, the intersection of remaining markers (4,340 SNPs) across populations was used for all subsequent analyses.

An additive genomic relationship matrix (G) was calculated using the R package AGHmatrix (Amadeu et al. 2016) using the VanRaden parameterization (VanRaden 2008; Isik et al. 2017):

$$ \frac{ZZ'}{2\Sigma p_i(1 - p_i)} $$

where ZZ' is the marker-based covariance matrix that reflects identity-by-descent relationships among the genotyped individuals and pᵢ is the minor allele frequency at each marker.

Genomic inbreeding coefficients for each individual were estimated as Gᵢᵢ − 1, where Gᵢᵢ is the diagonal element of G

corresponding to individual i. Individuals with inbreeding coefficients less than −0.2 or greater than 0.2 were removed. Additionally, closely related individuals within each population (Gᵢᵢ > 0.6) as well as individuals with suspected genotype sampling errors were identified and removed. After quality control based on inbreeding and pairwise relatedness, 665 G2-DE, 1,201 G3-DE, 698 G2-NC, and 613 G3-NC individuals remained. The additive genomic relationship matrix (G) was recalculated for this subset of lines. In addition, the matrix of realized dominance genomic relationships (D) was also calculated using the AGHmatrix package (Vitezica et al. 2013; Amadeu et al. 2016). Genetic relationship matrices were calculated for each population (DE or NC) separately as well as for both populations combined.

The genomic fixation index (Fst) was estimated between all pairs of the G2-DE and G3-DE and G2-NC and G3-NC populations using the R package hierstat (Goudet and Jombart 2015). Additionally, the expected heterozygosity was calculated for all of the populations using the adegenet R package (Jombart 2008).

## Initial genomic prediction to select G3 individuals for evaluation

The effectiveness of genomic selection in the off-season nursery was assessed using a subset of S₀.₁ families of G3 predicted plants. These were evaluated in 2020 at the original locations of phenotypic selection (target environments). Field space constraints limited the total number of progenies that could be evaluated in summer 2020. Therefore, we chose progenies from both the earliest and latest tails of genomic predictions as well as randomly sampled individuals for evaluation. Time constraints between processing the genotypic data (extracting DNA, sequencing, calling SNPs, quality control the data, and predicting GBLUPs) and planting time for the summer of 2020 limited
the range of initial prediction models that could be fit to the training data to select tails of the prediction distributions. Because previous research indicates that the genetic architecture of maize flowering time appears mostly additive with a distribution of allelic effects ranging from very small to moderate (effects up to 5% of total phenotypic variation in maize) (Buckler et al. 2009; Wisser et al. 2019), we fit two prediction models with distinct prior distributions of allelic effects. Both models have the form:

$$Y = \mu + M\beta + \varepsilon; \text{ [GS Model RR or BC]}$$

where Y is the vector of observations (data on individual plants of the G2 training sets), M is a matrix with columns corresponding to minor allele counts at each marker and rows corresponding to individuals, β is the vector of marker effects, and ε is the residual variance distributed as N(0, σ² residual). The marker effects were fit under Bayesian Ridge Regression (RR), which assumes a common distribution of polygenic marker effects, or BayesC (BC), which shrinks some marker effects to zero and allows other markers to have larger effects (de los Campos et al. 2013), using the default priors in the BGLR package (Pérez et al. 2014). The genomic prediction for each individual is the sum of the products of its minor allele counts times the marker effects over the total number of markers.

Predictions for DTA and DTS were standardized, and the mean standardized male and female flowering time was computed for the separate models trained on each G2 population. For each population in the off-season nursery (TropicS-DE-G3 and TropicS-NC-G3), the predictions were used to select individuals in the earliest and latest tails of genomic predictions as well as randomly sampled individuals. In total, 443 S0.1 families derived from each G3 population (DE and NC) were chosen based on ridge regression and BayesC predictions. The testing set based on G3 predictions from G3-DE included 42 individuals in the early tail predicted by both ridge regression and BayesC models, eight individuals in the early tail predicted by one of the two models, 35 individuals in the late tail predicted by both models, and 15 in the late tail predicted by one of the two models. We also chose 320 individuals at random from the remaining candidates. The same procedure was replicated for G3-NC.

To measure genomic prediction ability in a new year and generation, the testing set descended by self-fertilization from the chosen individuals within each of the G3-DE and G3-NC populations was planted in a common garden study during the summer of 2020 in Newark, DE and Clayton, NC. In addition, bulk population samples of the G2 and G3 generations of both DE and NC populations were planted to estimate genetic gain from phenotypic selection performed in summer 2019. The experimental design was an alpha design with two replicates of 40 incomplete blocks with 25 plots in each block in each location. Row and column coordinates of each plot in the field layout were used for spatial analysis.

DTA and DTS were measured as the number of days between planting and when at least 50% of the plants in a plot had visible anthers or silks, respectively. Plant and ear heights were measured on three plants within each family. Plant height was measured as the distance from the flag leaf node to the soil and ear height was measured as the distance from the primary ear-bearing node to the soil. For the G2 and G3 population entries, the traits were measured on about 12 individual plants per plot.

To estimate best linear unbiased estimators (BLUEs) for the lines in each location, a linear mixed model was fitted using ASReml-R as (Butler et al. 2017):

$$Y_{ijkl} = \mu + L_i + LR_{ij} + LRB_{ijk} + G_l + LG_{il} + \varepsilon; \text{ [EV Model 1]}$$

where Yijkl is the vector of observations, L_i is the random effect of the ith location, LR_{ij} is the random effect of the jth rep nested in the ith location, LRB_{ijk} is the random effect of the kth block nested in the jth rep in the Ith location; G_l is the fixed effect of the lth line, LG_{il} is the random interaction between line l and location i, and ε is the experimental error. The models were fit with heterogeneous residual variances among the locations.

Additionally, a spatial analysis was fitted using the same model as above except with spatially correlated residuals in row and columns directions (AR1×AR1; EV model 2). Residual effects ε within location i were distributed as σ² ei ∑ c (ρci) ⊗ ∑ r (ρri), where ∑ c (ρci) is the correlation matrix for columns and ρci is the auto-correlation parameter in the column direction and ∑ r (ρri) is the correlation matrix for rows and ρri is the auto-correlation parameter in the row direction within location i. The Bayesian Information Criterion (Schwarz 1978) was used to choose between EV models 1 and 2, and BLUEs for each line at each location were estimated from the chosen model.

## Retroactive assessment of alternative genomic prediction models

After the testing set was already defined for subsequent field evaluations based on the initial genomic prediction models (described above), we used the training sets to fit additional prediction models to assess how modeling field spatial effects and dominance effects would have changed the predictions. These models were fit using relationship matrices instead of marker incidence matrices to obtain individual genomic predictions (GBLUPs). A model accounting for field block effects was fit to the training data:

$$Y_{klm} = \mu + B_k + A_l + \varepsilon; \text{ [GS Model 2]}$$

where Yklm is the vector of observations, B_k is the random effect of the kth block with distribution N(0, σ²), A_l is the genomic breeding value for the lth individual with distribution N(0, σ²), where G is the estimated additive genomic relationship matrix and σ² is the additive genetic variance component, and ε is the residual error distributed as N(0, σ² residual).

GBLUPs from model 2 were predicted as:

$$GBLUP_l = \mu + \hat{A}_l$$

In both locations of the 2019 training trials, heterogeneity within blocks was observed visually, motivating the use of a spatial analysis to potentially improve the estimation of genetic effects. Plants were tagged in a serpentine pattern within blocks in specific rows and columns, so maps of spatial coordinates of plants were created. The spatial analysis model was fitted like the models above except with spatially autocorrelated residuals in row and column directions (AR1×AR1) only:

$$Y_{kl} = \mu + B_k + A_l + \varepsilon; \text{ [GS Model 3]}$$

where Ykl is the vector of observations, B_k is the random effect of the kth block, A_l is the genomic breeding value for the lth
the range of initial prediction models that could be fit to the training data to select tails of the prediction distributions. Because previous research indicates that the genetic architecture of maize flowering time appears mostly additive with a distribution of allelic effects ranging from very small to moderate (effects up to 5% of total phenotypic variation in maize) (Buckler et al. 2009; Wisser et al. 2019), we fit two prediction models with distinct prior distributions of allelic effects. Both models have the form:

$$Y = \mu + M\beta + \varepsilon; \quad [\text{GS Model RR or BC}]$$

where Y is the vector of observations (data on individual plants of the G2 training sets), M is a matrix with columns corresponding to minor allele counts at each marker and rows corresponding to individuals, β is the vector of marker effects, and ε is the residual variance distributed as N(0, σ² residual). The marker effects were fit under Bayesian Ridge Regression (RR), which assumes a common distribution of polygenic marker effects, or BayesC (BC), which shrinks some marker effects to zero and allows other markers to have larger effects (de los Campos et al. 2013), using the default priors in the BGLR package (Pérez et al. 2014). The genomic prediction for each individual is the sum of the products of its minor allele counts times the marker effects over the total number of markers.

Predictions for DTA and DTS were standardized, and the mean standardized male and female flowering time was computed for the separate models trained on each G2 population. For each population in the off-season nursery (TropicS-DE-G3 and TropicS-NC-G3), the predictions were used to select individuals in the earliest and latest tails of genomic predictions as well as randomly sampled individuals. In total, 443 S0.1 families derived from each G3 population (DE and NC) were chosen based on ridge regression and BayesC predictions. The testing set based on G3 predictions from G3-DE included 42 individuals in the early tail predicted by both ridge regression and BayesC models, eight individuals in the early tail predicted by one of the two models, 35 individuals in the late tail predicted by both models, and 15 in the late tail predicted by one of the two models. We also chose 320 individuals at random from the remaining candidates. The same procedure was replicated for G3-NC.

To measure genomic prediction ability in a new year and generation, the testing set descended by self-fertilization from the chosen individuals within each of the G3-DE and G3-NC populations was planted in a common garden study during the summer of 2020 in Newark, DE and Clayton, NC. In addition, bulk population samples of the G2 and G3 generations of both DE and NC populations were planted to estimate genetic gain from phenotypic selection performed in summer 2019. The experimental design was an alpha design with two replicates of 40 incomplete blocks with 25 plots in each block in each location. Row and column coordinates of each plot in the field layout were used for spatial analysis.

DTA and DTS were measured as the number of days between planting and when at least 50% of the plants in a plot had visible anthers or silks, respectively. Plant and ear heights were measured on three plants within each family. Plant height was measured as the distance from the flag leaf node to the soil and ear height was measured as the distance from the primary ear-bearing node to the soil. For the G2 and G3 population entries, the traits were measured on about 12 individual plants per plot.

To estimate best linear unbiased estimators (BLUEs) for the lines in each location, a linear mixed model was fitted using ASReml-R as (Butler et al. 2017):

$$Y_{ijkl} = \mu + L_i + LR_{ij} + LRB_{ijk} + G_l + LG_{il} + \varepsilon; \quad [\text{EV Model 1}]$$

where Yijkl is the vector of observations, L_i is the random effect of the ith location, LR_{ij} is the random effect of the jth rep nested in the ith location, LRB_{ijk} is the random effect of the kth block nested in the jth rep in the Ith location; G_l is the fixed effect of the lth line, LG_{il} is the random interaction between line l and location i, and ε is the experimental error. The models were fit with heterogeneous residual variances among the locations.

Additionally, a spatial analysis was fitted using the same model as above except with spatially correlated residuals in row and columns directions (AR1×AR1; EV model 2). Residual effects ε within location i were distributed as σ² ei ∑ c (ρci) ⊗ ∑ r (ρri), where ∑ c (ρci) is the correlation matrix for columns and ρci is the auto-correlation parameter in the column direction and ∑ r (ρri) is the correlation matrix for rows and ρri is the auto-correlation parameter in the row direction within location i. The Bayesian Information Criterion (Schwarz 1978) was used to choose between EV models 1 and 2, and BLUEs for each line at each location were estimated from the chosen model.

## Retroactive assessment of alternative genomic prediction models

After the testing set was already defined for subsequent field evaluations based on the initial genomic prediction models (described above), we used the training sets to fit additional prediction models to assess how modeling field spatial effects and dominance effects would have changed the predictions. These models were fit using relationship matrices instead of marker incidence matrices to obtain individual genomic predictions (GBLUPs). A model accounting for field block effects was fit to the training data:

$$Y_{klm} = \mu + B_k + A_l + \varepsilon; \quad [\text{GS Model 2}]$$

where Yklm is the vector of observations, B_k is the random effect of the kth block with distribution N(0, σ²_B), A_l is the genomic breeding value for the lth individual with distribution N(0, σ²_A), where G is the estimated additive genomic relationship matrix and σ²_A is the additive genetic variance component, and ε is the residual error distributed as N(0, σ² residual).

GBLUPs from model 2 were predicted as:

$$GBLUP_l = \mu + \hat{A}_l$$

In both locations of the 2019 training trials, heterogeneity within blocks was observed visually, motivating the use of a spatial analysis to potentially improve the estimation of genetic effects. Plants were tagged in a serpentine pattern within blocks in specific rows and columns, so maps of spatial coordinates of plants were created. The spatial analysis model was fitted like the models above except with spatially autocorrelated residuals in row and column directions (AR1×AR1) only:

$$Y_{kl} = \mu + B_k + A_l + \varepsilon; \quad [\text{GS Model 3}]$$

where Ykl is the vector of observations, B_k is the random effect of the kth block, A_l is the genomic breeding value for the lth
individual with distribution $A_l \sim N(0, G\sigma_A^2)$, where G is the estimate additive genomic relationship matrix for the population and $\sigma_A^2$ is the additive genetic variance component, and $\varepsilon$ is the residual error with distribution $\varepsilon \sim N(0, \sigma_e^2 \Sigma_c (\rho_c) \otimes \Sigma_r (\rho_r))$, where $\Sigma_c (\rho_c)$ is the correlation matrix for residual effects across columns and $\rho_c$ is the auto-correlation parameter in the column direction and $\Sigma_r (\rho_r)$ is the correlation matrix for residual effects across rows and $\rho_r$ is the auto-correlation parameter in the row direction (Isik et al. 2017).

A fourth model was fit using the dominance relationship matrix (D) to evaluate if modeling dominance effects improves prediction ability. The following model was fit to the data for each population separately:

$$Y_{kl} = \mu + B_k + A_l + D_l + \varepsilon; \quad [\text{GS Model 4}]$$

where $Y_{kl}$ is the vector of observations, $B_k$ is the effect of the kth block, $A_l$ is the genomic breeding value for the lth individual with distribution $A_l \sim N(0, G\sigma_A^2)$, where G is the estimated additive genomic relationship matrix for the population and $\sigma_A^2$ is the additive genetic variance component, $D_l$ is the estimated dominance effect for the lth individual with distribution $D_l \sim N(0, D\sigma_D^2)$ and $\sigma_D^2$ is the dominance genetic variance component, and $\varepsilon$ is the residual error with distribution $\varepsilon \sim N(0, \sigma_{\text{residual}}^2)$.

GBLUPs were predicted from Model 4 in two ways:

i) GBLUP$_l = \mu + \hat{A}_l$, the breeding value prediction, and

ii) GBLUP$_l = \mu + \hat{A}_l + \frac{1}{2}\hat{D}_l$, the predicted genotypic value of the $S_{0.1}$ progeny resulting from self-fertilization of individual l.

## Evaluation of genomic prediction ability

Genomic prediction ability for the different models was estimated with Pearson correlations between genomic predictions made from the 2019 training sets (individual plants) and the EV Model 1 BLUEs from the 2020 testing set ($S_{0.1}$ families). Recall that the training and testing data are similar in terms of the locations used but differ by the nature of their genetic constitution (segregating individuals vs selfed progeny rows) and year of evaluation (2019 vs 2020). In addition, training sets were specific to the population and corresponding location of selection, but the testing set included samples from all populations evaluated together at both locations in 2020. Therefore, models trained on TropicS-G2-DE and TropicS-G2-NC populations from 2019 were evaluated based on TropicS-G3-DE and TropicS-G3-NC $S_{0.1}$ families measured in environments of 2020 for each of the corresponding locations. In addition, cross-environment and cross-population predictions were evaluated. For cross-environment analysis, for example, prediction ability for GBLUPs trained on the 2019 TropicS-G2-DE data was determined using TropicS-G3-DE family data from the 2020 NC location, and vice versa for the NC sets (i.e. prediction abilities for models trained on each population were tested across locations). Cross-population predictions involved predicting the G3 generation of one population based on the G2 generation of the other lineage. For example, models trained on TropicS-G2-DE were evaluated on TropicS-G3-NC individuals used as the prediction set.

## Evaluation of phenotypic prediction ability

We also estimated the correlations between pollination date of G3 $S_0$ plants in the off-season nursery and the mean flowering times of their G3 $S_{0.1}$ progeny lines in the summer 2020 environments. Pollination dates serve as proxies for DTS of the parental plants

and permit estimation of the effectiveness of phenotypic selection for earlier plants in the off-season nursery.

## Measuring prediction ability with randomly sampled subset of markers

We selected random subsets of markers to determine if the number of makers was limiting prediction ability. Subsets of 10%, 25%, 50%, 75%, and 90% of the markers were randomly sampled 10 times and GS Model 2 was refit each time to generate predictions. Prediction ability was measured as the pairwise correlation between the predictions and the BLUEs for each population representing a combination of population and evaluation location.

## Reliability of selection candidates

Reliability is the squared correlation between true and predicted effects ($\hat{r}_{g,\hat{g}}^2$). The reliability of predictions for the selection candidates in the off-season nursery was estimated as:

$$\hat{r}_{g,\hat{g}}^2 = 1 - \frac{\text{PEV}}{\hat{\sigma}_g^2}; \quad (1)$$

where $\hat{r}_{g,\hat{g}}^2$ is the prediction error variance and $\hat{\sigma}_g^2$ is the estimated genetic variance (Isik et al. 2017).

This estimate of reliability is related to heritability by:

$$E(\hat{r}_{g,\hat{g}}^2) = \frac{\sigma_g^2}{\sigma_P^2} = h^2; \quad (2)$$

where $\sigma_g^2$ is the phenotypic variance due to the genotype, $\sigma_P^2$ is the total phenotypic variance, and $h^2$ is heritability. By averaging reliability across the selected BLUPS, a generalized heritability is obtained and applicable for predicting responses to selection based on those predictions (Isik et al. 2017; Schmidt et al. 2019). Here, GS Model 2 was used to estimate generalized heritabilities for selection candidates. Additionally, the reliability of the test set $S_{0.1}$ families was calculated using EV Model 2 with the lines fit as random effects.

## Estimating theoretical maximum prediction ability

The estimated theoretical maximum prediction ability for the prediction models was calculated as:

$$\hat{r}_{G(19,20)};$$

where $\hat{r}_{G(19,20)}$ is the genetic correlation between 2019 and 2020 and h is the square root of the mean reliability of the selection candidates in the target environment described above. Essentially, $\hat{r}_{G(19,20)}$ is an estimate for the true correlation of the genetic values across years. However, the true genetic values of 2020 are unknown but can be estimated with BLUEs. These BLUEs are related to their true breeding values through the square root of the heritability.

To estimate the genetic correlation between 2019 and 2020, a separate set of inbred lines derived without intentional selection from the TropicS-NC-G1 population was planted in 2019 and 2020 in both locations. In 2019, 205 inbreds were planted in an alpha lattice design with two replicates with 18 blocks, containing 12 entries per block. In 2020, a subset of 100 of the inbreds was planted within the alpha lattice design used for the testing set of
individual with distribution $A_l \sim N(0, G\sigma_A^2)$, where G is the estimate additive genomic relationship matrix for the population and $\sigma_A^2$ is the additive genetic variance component, and $\epsilon$ is the residual error with distribution $\epsilon \sim N(0, \sigma_\epsilon^2 \Sigma_c(\rho_c) \otimes \Sigma_r(\rho_r))$, where $\Sigma_c(\rho_c)$ is the correlation matrix for residual effects across columns and $\rho_c$ is the auto-correlation parameter in the column direction and $\Sigma_r(\rho_r)$ is the correlation matrix for residual effects across rows and $\rho_r$ is the auto-correlation parameter in the row direction (Isik et al. 2017).

A fourth model was fit using the dominance relationship matrix (D) to evaluate if modeling dominance effects improves prediction ability. The following model was fit to the data for each population separately:

$$Y_{kl} = \mu + B_k + A_l + D_l + \epsilon; \quad [\text{GS Model 4}]$$

where $Y_{kl}$ is the vector of observations, $B_k$ is the effect of the kth block, $A_l$ is the genomic breeding value for the lth individual with distribution $A_l \sim N(0, G\sigma_A^2)$, where G is the estimated additive genomic relationship matrix for the population and $\sigma_A^2$ is the additive genetic variance component, $D_l$ is the estimated dominance effect for the lth individual with distribution $D_l \sim N(0, D\sigma_D^2)$ and $\sigma_D^2$ is the dominance genetic variance component, and $\epsilon$ is the residual error with distribution $\epsilon \sim N(0, \sigma_{\text{residual}}^2)$.

GBLUPs were predicted from Model 4 in two ways:

i) GBLUP$_l = \mu + \hat{A}_l$, the breeding value prediction, and

ii) GBLUP$_l = \mu + \hat{A}_l + \frac{1}{2}\hat{D}_l$, the predicted genotypic value of the $S_{0.1}$ progeny resulting from self-fertilization of individual l.

## Evaluation of genomic prediction ability

Genomic prediction ability for the different models was estimated with Pearson correlations between genomic predictions made from the 2019 training sets (individual plants) and the EV Model 1 BLUEs from the 2020 testing set ($S_{0.1}$ families). Recall that the training and testing data are similar in terms of the locations used but differ by the nature of their genetic constitution (segregating individuals vs selfed progeny rows) and year of evaluation (2019 vs 2020). In addition, training sets were specific to the population and corresponding location of selection, but the testing set included samples from all populations evaluated together at both locations in 2020. Therefore, models trained on TropicS-G2-DE and TropicS-G2-NC populations from 2019 were evaluated based on TropicS-G3-DE and TropicS-G3-NC $S_{0.1}$ families measured in environments of 2020 for each of the corresponding locations. In addition, cross-environment and cross-population predictions were evaluated. For cross-environment analysis, for example, prediction ability for GBLUPs trained on the 2019 TropicS-G2-DE data was determined using TropicS-G3-DE family data from the 2020 NC location, and vice versa for the NC sets (i.e. prediction abilities for models trained on each population were tested across locations). Cross-population predictions involved predicting the G3 generation of one population based on the G2 generation of the other lineage. For example, models trained on TropicS-G2-DE were evaluated on TropicS-G3-NC individuals used as the prediction set.

## Evaluation of phenotypic prediction ability

We also estimated the correlations between pollination date of G3 $S_0$ plants in the off-season nursery and the mean flowering times of their G3 $S_{0.1}$ progeny lines in the summer 2020 environments. Pollination dates serve as proxies for DTS of the parental plants

and permit estimation of the effectiveness of phenotypic selection for earlier plants in the off-season nursery.

## Measuring prediction ability with randomly sampled subset of markers

We selected random subsets of markers to determine if the number of makers was limiting prediction ability. Subsets of 10%, 25%, 50%, 75%, and 90% of the markers were randomly sampled 10 times and GS Model 2 was refit each time to generate predictions. Prediction ability was measured as the pairwise correlation between the predictions and the BLUEs for each population representing a combination of population and evaluation location.

## Reliability of selection candidates

Reliability is the squared correlation between true and predicted effects ($\hat{r}_{g,\hat{g}}^2$). The reliability of predictions for the selection candidates in the off-season nursery was estimated as:

$$\hat{r}_{g,\hat{g}}^2 = 1 - \frac{\text{PEV}}{\hat{\sigma}_g^2}; \quad (1)$$

where $\hat{r}_{g,\hat{g}}^2$ is the prediction error variance and $\hat{\sigma}_g^2$ is the estimated genetic variance (Isik et al. 2017).

This estimate of reliability is related to heritability by:

$$E(\hat{r}_{g,\hat{g}}^2) = \frac{\sigma_g^2}{\sigma_P^2} = h^2; \quad (2)$$

where $\sigma_g^2$ is the phenotypic variance due to the genotype, $\sigma_P^2$ is the total phenotypic variance, and $h^2$ is heritability. By averaging reliability across the selected BLUPS, a generalized heritability is obtained and applicable for predicting responses to selection based on those predictions (Isik et al. 2017; Schmidt et al. 2019). Here, GS Model 2 was used to estimate generalized heritabilities for selection candidates. Additionally, the reliability of the test set $S_{0.1}$ families was calculated using EV Model 2 with the lines fit as random effects.

## Estimating theoretical maximum prediction ability

The estimated theoretical maximum prediction ability for the prediction models was calculated as:

$$\hat{r}_{G(19,20)}h;$$

where $\hat{r}_{G(19,20)}$ is the genetic correlation between 2019 and 2020 and h is the square root of the mean reliability of the selection candidates in the target environment described above. Essentially, $\hat{r}_{G(19,20)}$ is an estimate for the true correlation of the genetic values across years. However, the true genetic values of 2020 are unknown but can be estimated with BLUEs. These BLUEs are related to their true breeding values through the square root of the heritability.

To estimate the genetic correlation between 2019 and 2020, a separate set of inbred lines derived without intentional selection from the TropicS-NC-G1 population was planted in 2019 and 2020 in both locations. In 2019, 205 inbreds were planted in an alpha lattice design with two replicates with 18 blocks, containing 12 entries per block. In 2020, a subset of 100 of the inbreds was planted within the alpha lattice design used for the testing set of
S0.1 families. In both years, DTA, DTS, plant height, and ear height were measured.

The following model was used to estimate the genetic correlations between environments based on the TropicS-NC-G1 derived inbred lines. The model was fit as:

$$Y_{ijkl} = \mu + E_i + ER_{ij} + ERB_{ijk} + EG_{il} + \varepsilon;$$

where $Y_{ijkl}$ is the vector of observations, $E_i$ is the fixed effect of the ith environment (location by year combination), $ER_{ij}$ is the random effect of the jth replication nested in the ith environment, $ERB_{ijk}$ is the random effect of the kth block nested in the jth replication in the ith environment; $EG_{il}$ is the random effect of the lth inbred line in the ith environment, and $\varepsilon$ is the residual error. $EG_{il}$ effects were modeled with an unstructured covariance matrix that included separate genetic variances within each of the four environments and separate genetic covariances for all six pairs of environments, resulting in six estimates of the genetic correlation ($\hat{r}_g$) between environment pairs.

## Estimating genetic gain from phenotypic and genomic selection

Genetic gain for the phenotypic selection that occurred in the summer of 2019 was estimated as the difference between BLUEs for the G3 and G2 populations. This response was estimated within the same location where selection occurred (i.e. TropicS-G2-DE and TropicS-G3-DE grown in DE and TropicS-G2-NC and TropicS-G3-NC grown in NC). Genetic gain from genomic selection was estimated as the difference between the mean flowering times of the lines derived from individuals in the earliest 5th percentile of G3 based on initial genomic prediction models and the mean of the randomly chosen individuals. The earliest 5th percentile among 768 genotyped plants in the G3-NC population was represented by 38 individuals in the earliest 5th percentile in both RR and BayesC models. The earliest tail for the DE population was represented by the 42 individuals in the earliest tail in both models plus 15 additional individuals in the earliest tail in each of the models individually. This set of 72 is approximately 5% of the total sample size of 1,536 individuals genotyped in the G3-DE population.

## Results and discussion

### Population differentiation and expected heterozygosity

Population differentiation, measured by the fixation index ($F_{st}$) between the G2-DE and G2-NC populations used for selection, was 0.003. Previous phenotypic selection for early flowering within DE and NC independently appeared to have changed genome-wide allele frequencies very slowly or in largely the same direction. $F_{st}$ between generations of the two lineages also showed low differentiation ($F_{st} = 0.003$ and 0.002 between G2 and G3 for the DE and NC lineages, respectively), reflecting minimal impacts of selection or drift on allele frequency changes within and between lineages. This indicates that neither drift nor selection separated the populations very much during phenotypic selection performed prior to the current study.

The expected heterozygosity ($H_e$) was calculated within populations of each lineage to determine if selection changed the diversity in the populations used for the current study. For G2-DE and G3-DE, $H_e$ was 0.34 and 0.33, respectively. Similarly, $H_e$ was 0.34 for G2-NC and 0.33 for G3-NC. In both populations, diversity

decreased slightly due to a generation of selection. This is consistent with the finding that genome-wide diversity during recurrent phenotypic selection can be largely maintained in the short term (Wisser et al. 2019).

### Prediction ability for genomic selection

We used different prior distributions of marker effects (Bayesian Ridge Regression or BayesC) for the initial genomic selection models in case genomic predictions were sensitive to the genetic architecture of flowering time in these populations. We also compared the GBLUPS for individuals from the G3 populations with the BLUES for their selfed progenies (S0.1 families) to measure the prediction ability for both models. For both G3-DE and G3-NC, prediction ability was not affected by the choice of model (Supplementary Figs. 1 and 2). Prediction abilities for DTA were 0.39 and 0.29 for G3-DE and G3-NC, respectively, for both Ridge Regression and BayesC. These results suggest that, in these populations, genetic variation for flowering time is largely due to polygenic effects with limited influence from loci with strong effects.

Additional prediction models were fit retroactively to the G2 training sets to determine if field spatial effects with varying levels of complexity could have improved prediction ability. Genomic selection Model 2 including blocks as the only spatial effect had the best fit (lower BIC) for DTS, plant height, and ear height measured in the G2-DE training set, whereas a more complicated spatial analysis with residual correlations due to proximity in both row and column directions (GS Model 3) had a better fit for DTA (Supplementary Table 1). For the G2-NC training set, Model 2 was better for DTA, DTS, and ear height, while Model 3 had a better fit statistic for plant height (Supplementary Table 1). Models including with spatially correlated residuals (EV Model 2) had better fit than models with uncorrelated residuals but block main effects (EV Model 1) in the evaluation trials (Supplementary Table 2). For consistency in calculating GBLUPS retroactively, GS Model 2 with random block effects was used for both training sets and all traits. Prediction abilities were higher by only 0.01 to 0.02 when blocks were added to the prediction model.

For flowering time, there was a significant positive correlation (ranging from 0.30 to 0.40) between the Model 2 S0 GBLUPS and the S1 BLUEs in both populations in both target environments (Fig. 2, Supplementary Fig. 3). Prediction abilities were higher for G3-DE compared to G3-NC for both DTA and DTS. This pattern held for plant and ear height, for which prediction ability was higher in the DE population compared to the NC population (Supplementary Figs. 4 and 5). These results suggest that within-field heterogeneity was less of a problem in the 2019 DE than in the 2019 NC training experiments, resulting in better modeling of genomic breeding values in the 2019 DE environment.

We tested if models trained on one lineage could predict individuals from the other lineage. For this, G2 training data from one of the reference populations (e.g. G2-DE) was used to predict next generation progeny from the other population (i.e. G3-NC). Prediction abilities for flowering time were consistently lower using the training population from the alternative lineage. Specifically, prediction abilities for the DE lineage decreased by 0.1 and 0.07 for DTA and DTS, respectively, when the models were trained on G2-NC. Similarly, prediction abilities decreased for the NC lineage by 0.05 (DTA and DTS) for models trained on G2-DE. The decreases in prediction ability for cross-lineage predictions were much greater than those for cross-environment predictions within each lineage (Table 1). Although the $F_{st}$ estimates indicated that the populations are not very diverged after two generations of local phenotypic selection, even a small amount of
S0.1 families. In both years, DTA, DTS, plant height, and ear height were measured.

The following model was used to estimate the genetic correlations between environments based on the TropicS-NC-G1 derived inbred lines. The model was fit as:

$$Y_{ijkl} = \mu + E_i + ER_{ij} + ERB_{ijk} + EG_{il} + \varepsilon;$$

where $Y_{ijkl}$ is the vector of observations, $E_i$ is the fixed effect of the ith environment (location by year combination), $ER_{ij}$ is the random effect of the jth replication nested in the ith environment, $ERB_{ijk}$ is the random effect of the kth block nested in the jth replication in the ith environment; $EG_{il}$ is the random effect of the lth inbred line in the ith environment, and $\varepsilon$ is the residual error. $EG_{il}$ effects were modeled with an unstructured covariance matrix that included separate genetic variances within each of the four environments and separate genetic covariances for all six pairs of environments, resulting in six estimates of the genetic correlation ($\hat{r}_g$) between environment pairs.

## Estimating genetic gain from phenotypic and genomic selection

Genetic gain for the phenotypic selection that occurred in the summer of 2019 was estimated as the difference between BLUEs for the G3 and G2 populations. This response was estimated within the same location where selection occurred (i.e. TropicS-G2-DE and TropicS-G3-DE grown in DE and TropicS-G2-NC and TropicS-G3-NC grown in NC). Genetic gain from genomic selection was estimated as the difference between the mean flowering times of the lines derived from individuals in the earliest 5th percentile of G3 based on initial genomic prediction models and the mean of the randomly chosen individuals. The earliest 5th percentile among 768 genotyped plants in the G3-NC population was represented by 38 individuals in the earliest 5th percentile in both RR and BayesC models. The earliest tail for the DE population was represented by the 42 individuals in the earliest tail in both models plus 15 additional individuals in the earliest tail in each of the models individually. This set of 72 is approximately 5% of the total sample size of 1,536 individuals genotyped in the G3-DE population.

## Results and discussion

### Population differentiation and expected heterozygosity

Population differentiation, measured by the fixation index ($F_{st}$) between the G2-DE and G2-NC populations used for selection, was 0.003. Previous phenotypic selection for early flowering within DE and NC independently appeared to have changed genome-wide allele frequencies very slowly or in largely the same direction. $F_{st}$ between generations of the two lineages also showed low differentiation ($F_{st} = 0.003$ and 0.002 between G2 and G3 for the DE and NC lineages, respectively), reflecting minimal impacts of selection or drift on allele frequency changes within and between lineages. This indicates that neither drift nor selection separated the populations very much during phenotypic selection performed prior to the current study.

The expected heterozygosity ($H_e$) was calculated within populations of each lineage to determine if selection changed the diversity in the populations used for the current study. For G2-DE and G3-DE, $H_e$ was 0.34 and 0.33, respectively. Similarly, $H_e$ was 0.34 for G2-NC and 0.33 for G3-NC. In both populations, diversity

decreased slightly due to a generation of selection. This is consistent with the finding that genome-wide diversity during recurrent phenotypic selection can be largely maintained in the short term (Wisser et al. 2019).

### Prediction ability for genomic selection

We used different prior distributions of marker effects (Bayesian Ridge Regression or BayesC) for the initial genomic selection models in case genomic predictions were sensitive to the genetic architecture of flowering time in these populations. We also compared the GBLUPS for individuals from the G3 populations with the BLUES for their selfed progenies (S0.1 families) to measure the prediction ability for both models. For both G3-DE and G3-NC, prediction ability was not affected by the choice of model (Supplementary Figs. 1 and 2). Prediction abilities for DTA were 0.39 and 0.29 for G3-DE and G3-NC, respectively, for both Ridge Regression and BayesC. These results suggest that, in these populations, genetic variation for flowering time is largely due to polygenic effects with limited influence from loci with strong effects.

Additional prediction models were fit retroactively to the G2 training sets to determine if field spatial effects with varying levels of complexity could have improved prediction ability. Genomic selection Model 2 including blocks as the only spatial effect had the best fit (lower BIC) for DTS, plant height, and ear height measured in the G2-DE training set, whereas a more complicated spatial analysis with residual correlations due to proximity in both row and column directions (GS Model 3) had a better fit for DTA (Supplementary Table 1). For the G2-NC training set, Model 2 was better for DTA, DTS, and ear height, while Model 3 had a better fit statistic for plant height (Supplementary Table 1). Models including with spatially correlated residuals (EV Model 2) had better fit than models with uncorrelated residuals but block main effects (EV Model 1) in the evaluation trials (Supplementary Table 2). For consistency in calculating GBLUPS retroactively, GS Model 2 with random block effects was used for both training sets and all traits. Prediction abilities were higher by only 0.01 to 0.02 when blocks were added to the prediction model.

For flowering time, there was a significant positive correlation (ranging from 0.30 to 0.40) between the Model 2 S0 GBLUPS and the S1 BLUEs in both populations in both target environments (Fig. 2, Supplementary Fig. 3). Prediction abilities were higher for G3-DE compared to G3-NC for both DTA and DTS. This pattern held for plant and ear height, for which prediction ability was higher in the DE population compared to the NC population (Supplementary Figs. 4 and 5). These results suggest that within-field heterogeneity was less of a problem in the 2019 DE than in the 2019 NC training experiments, resulting in better modeling of genomic breeding values in the 2019 DE environment.

We tested if models trained on one lineage could predict individuals from the other lineage. For this, G2 training data from one of the reference populations (e.g. G2-DE) was used to predict next generation progeny from the other population (i.e. G3-NC). Prediction abilities for flowering time were consistently lower using the training population from the alternative lineage. Specifically, prediction abilities for the DE lineage decreased by 0.1 and 0.07 for DTA and DTS, respectively, when the models were trained on G2-NC. Similarly, prediction abilities decreased for the NC lineage by 0.05 (DTA and DTS) for models trained on G2-DE. The decreases in prediction ability for cross-lineage predictions were much greater than those for cross-environment predictions within each lineage (Table 1). Although the $F_{st}$ estimates indicated that the populations are not very diverged after two generations of local phenotypic selection, even a small amount of
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Validation of genomic prediction for days to anthesis in the DE-G3 and NC-G3 populations grown in DE and NC. For each population, S0 Genomic Best Linear Unbiased Predictors (GBLUPs) from GS Model 2 are plotted against S1 Best Linear Unbiased Estimators (BLUEs) for each population grown in 2020 in DE and NC target locations.</figcaption></figure>

<table><caption>**Table 1.** Prediction abilities for models using G2-DE and G2-NC as training population and either G3-DE or G3-NC as the testing population.</caption><thead><tr><th rowspan="3">Evaluation site</th><th colspan="4">DTA</th><th colspan="4">DTS</th></tr><tr><th colspan="2">G2-DE</th><th colspan="2">G2-NC</th><th colspan="2">G2-DE</th><th colspan="2">G2-NC</th></tr><tr><th colspan="4">Evaluation population</th><th colspan="2">Evaluation population</th><th colspan="2">Evaluation population</th></tr><tr><th></th><th>G3-DE</th><th>G3-NC</th><th>G3-DE</th><th>G3-NC</th><th>G3-DE</th><th>G3-NC</th><th>G3-DE</th><th>G3-NC</th></tr></thead><tbody><tr><td>DE</td><td>0.40</td><td>0.26</td><td>0.30</td><td>0.31</td><td>0.36</td><td>0.29</td><td>0.29</td><td>0.35</td></tr><tr><td>NC</td><td>0.38</td><td>0.25</td><td>0.28</td><td>0.30</td><td>0.33</td><td>0.30</td><td>0.26</td><td>0.35</td></tr></tbody></table>

G3 S0 BLUPs from each model were validated against S1 BLUEs grown in either DE or NC.

genetic differentiation combined with relatively small genotype-by-environment interaction variance had an observable impact on genomic prediction ability.

We also tested if modeling both additive and dominance effects improved prediction ability. In general, dominance variance was negligible for flowering time, but for plant height $\hat{\sigma}_D^2$ ranged from 12 to 25% of the magnitude of additive variance. However, its inclusion for either flowering time or plant and ear height had little impact on estimates of $\hat{\sigma}_A^2$ (Table 2). Consequently, including the dominance relationship matrix did not improve prediction ability for any trait in either population regardless of the method of including dominance effects (Table 2). Although inclusion of dominance effects improved height prediction ability in maize hybrids (Ramstein et al. 2020), the smaller contribution of dominance effects to partially inbred S0.1 lines compared to hybrid genetic variances limits its utility for predicting selfed progeny means.

## Limits to genomic prediction ability

For flowering time, we compared the realized prediction abilities for each location to an estimate of the theoretical maximum based on the between-environment genetic correlation adjusted by within-environment heritability. These genetic correlations

ranged from 0.79-0.98 (mean: 0.90; Supplementary Table 3), indicating modest to low genetic-by-environment interactions between training and testing environments. The theoretical maxima for prediction abilities were: 0.87 for DTA and DTS in G3-DE; 0.83 for DTA in G3-NC, and 0.88 for DTS in G3-NC. This is two to three times greater than estimates of the realized prediction abilities (Table 1), a disparity that could be due to prediction models, number of markers, composition of training population, and heritability of the traits.

Prediction accuracy was not limited by prediction models, as evidenced by the very small impact of different assumptions about genetic architecture (ridge regression assuming a normal distribution of effects, BayesC assuming a mixture distribution of effects, or modeling of both additive and dominance effects). We tested if the number of markers available limited our prediction ability in these populations by taking random subsets of 10% - 90% of markers (Yu et al. 2009) (Fig. 3). In most cases for both populations, prediction ability using only 10 to 50% of the markers was much less than when using 100% of the markers. For both populations, however, prediction ability plateaued as the proportion of markers reached ~90% of the 4,340 markers in the full data set (Fig. 3). The difference in mean prediction ability
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Validation of genomic prediction for days to anthesis in the DE-G3 and NC-G3 populations grown in DE and NC. For each population, S0 Genomic Best Linear Unbiased Predictors (GBLUPs) from GS Model 2 are plotted against S1 Best Linear Unbiased Estimators (BLUEs) for each population grown in 2020 in DE and NC target locations.</figcaption></figure>

<table><caption>**Table 1.** Prediction abilities for models using G2-DE and G2-NC as training population and either G3-DE or G3-NC as the testing population.</caption><thead><tr><th rowspan="3">Evaluation site</th><th colspan="4">DTA</th><th colspan="4">DTS</th></tr><tr><th colspan="2">G2-DE</th><th colspan="2">G2-NC</th><th colspan="2">G2-DE</th><th colspan="2">G2-NC</th></tr><tr><th colspan="4">Evaluation population</th><th colspan="2">Evaluation population</th><th colspan="2">Evaluation population</th></tr><tr><th></th><th>G3-DE</th><th>G3-NC</th><th>G3-DE</th><th>G3-NC</th><th>G3-DE</th><th>G3-NC</th><th>G3-DE</th><th>G3-NC</th></tr></thead><tbody><tr><td>DE</td><td>0.40</td><td>0.26</td><td>0.30</td><td>0.31</td><td>0.36</td><td>0.29</td><td>0.29</td><td>0.35</td></tr><tr><td>NC</td><td>0.38</td><td>0.25</td><td>0.28</td><td>0.30</td><td>0.33</td><td>0.30</td><td>0.26</td><td>0.35</td></tr></tbody></table>

G3 S0 BLUPs from each model were validated against S1 BLUEs grown in either DE or NC.

genetic differentiation combined with relatively small genotype-by-environment interaction variance had an observable impact on genomic prediction ability.

We also tested if modeling both additive and dominance effects improved prediction ability. In general, dominance variance was negligible for flowering time, but for plant height $\hat{\sigma}_D^2$ ranged from 12 to 25% of the magnitude of additive variance. However, its inclusion for either flowering time or plant and ear height had little impact on estimates of $\hat{\sigma}_A^2$ (Table 2). Consequently, including the dominance relationship matrix did not improve prediction ability for any trait in either population regardless of the method of including dominance effects (Table 2). Although inclusion of dominance effects improved height prediction ability in maize hybrids (Ramstein et al. 2020), the smaller contribution of dominance effects to partially inbred S0.1 lines compared to hybrid genetic variances limits its utility for predicting selfed progeny means.

## Limits to genomic prediction ability

For flowering time, we compared the realized prediction abilities for each location to an estimate of the theoretical maximum based on the between-environment genetic correlation adjusted by within-environment heritability. These genetic correlations

ranged from 0.79-0.98 (mean: 0.90; Supplementary Table 3), indicating modest to low genetic-by-environment interactions between training and testing environments. The theoretical maxima for prediction abilities were: 0.87 for DTA and DTS in G3-DE; 0.83 for DTA in G3-NC, and 0.88 for DTS in G3-NC. This is two to three times greater than estimates of the realized prediction abilities (Table 1), a disparity that could be due to prediction models, number of markers, composition of training population, and heritability of the traits.

Prediction accuracy was not limited by prediction models, as evidenced by the very small impact of different assumptions about genetic architecture (ridge regression assuming a normal distribution of effects, BayesC assuming a mixture distribution of effects, or modeling of both additive and dominance effects). We tested if the number of markers available limited our prediction ability in these populations by taking random subsets of 10% - 90% of markers (Yu et al. 2009) (Fig. 3). In most cases for both populations, prediction ability using only 10 to 50% of the markers was much less than when using 100% of the markers. For both populations, however, prediction ability plateaued as the proportion of markers reached ~90% of the 4,340 markers in the full data set (Fig. 3). The difference in mean prediction ability
<table><caption>**Table 2.** Variance component ($\hat{\sigma}^2$) and prediction ability estimates from GS model 2 with additive genetic effects ($r_A$) and GS Model 4 with additive and dominance genetic effects ($r_{A+D}$) based on G2 training data on individual plants grown in summer, 2019 in DE and NC.</caption><thead><tr><th rowspan="2">Parameter estimate</th><th colspan="4">G2-DE</th><th colspan="4">G2-NC</th></tr><tr><th>DTA</th><th>DTS</th><th>Plant height</th><th>Ear height</th><th>DTA</th><th>DTS</th><th>Plant height</th><th>Ear height</th></tr></thead><tbody><tr><td colspan="9"><strong>GS Model 2 with additive genetic effects only</strong></td></tr><tr><td>σ<sup>2</sup><sub>Block</sub></td><td>0.4</td><td>0.6</td><td>26.3</td><td>11.2</td><td>0.2</td><td>1.1</td><td>296.4</td><td>51.5</td></tr><tr><td>σ<sup>2</sup><sub>A</sub></td><td>3.1</td><td>5.7</td><td>236.5</td><td>133.1</td><td>2.5</td><td>3.7</td><td>128.5</td><td>81.2</td></tr><tr><td>σ<sup>2</sup><sub>e</sub></td><td>6.4</td><td>12.4</td><td>251.6</td><td>159.0</td><td>6.3</td><td>10.0</td><td>285.3</td><td>153.9</td></tr><tr><td>Prediction ability (r<sub>A</sub>)</td><td>0.40</td><td>0.36</td><td>0.47</td><td>0.50</td><td>0.30</td><td>0.35</td><td>0.42</td><td>0.35</td></tr><tr><td colspan="9"><strong>GS Model 4 with additive and dominance genetic effects</strong></td></tr><tr><td>σ<sup>2</sup><sub>Block</sub></td><td>0.4</td><td>0.6</td><td>26.4</td><td>11.4</td><td>0.2</td><td>1.1</td><td>295.5</td><td>51.4</td></tr><tr><td>σ<sup>2</sup><sub>A</sub></td><td>3.1</td><td>5.7</td><td>238.2</td><td>130.0</td><td>2.5</td><td>3.7</td><td>126.6</td><td>79.3</td></tr><tr><td>σ<sup>2</sup><sub>D</sub></td><td>0</td><td>0.2</td><td>30.2</td><td>16.0</td><td>0.1</td><td>0.4</td><td>31.7</td><td>15.6</td></tr><tr><td>σ<sup>2</sup><sub>e</sub></td><td>6.4</td><td>12.2</td><td>220.0</td><td>144.4</td><td>6.3</td><td>9.6</td><td>255.0</td><td>139.3</td></tr><tr><td>Prediction ability (r<sub>A+D</sub>)</td><td>0.40</td><td>0.36</td><td>0.47</td><td>0.50</td><td>0.30</td><td>0.35</td><td>0.42</td><td>0.35</td></tr></tbody></table>

Traits measured include days to anthesis (DTA), days to silking (DTS), plant height, and ear height.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Prediction abilities for days to anthesis (DTA) and days to silk (DTS) in the G3-DE and G3-NC populations based on models using random subsets of 10 to 100% of the full set of markers.</figcaption></figure>

between 90 and 100% of the full marker set was never greater than 0.01. These results indicate that prediction abilities are not limited by the available number of markers. Curiously, prediction abilities for both populations using only 10% of the total markers remained greater than $r=0.20$. This suggests that the disequilibrium between markers and causal alleles was sufficiently extensive in these populations so that relatively few markers can lead to a substantial positive prediction ability. Direct estimation of the linkage disequilibrium within each chromosome showed that it extends over Mbp distances in both populations (Supplementary Figs. 6, 7 and 8). Although most of the founder linkage disequilibrium was reduced during development of the TropicS-G0 base population, some large blocks remained in the population (Weldekidan et al. 2022). Moreover, the current study was initiated with populations that had undergone two generations of selection, which is also expected to increase linkage disequilibrium.

Genetic distance between training and testing sets was quite small, as indicated above by the values of genomic differentiation

near zero, and as expected for testing sets comprising selfed grandchildren (G3 $S_{0.1}$ families) of individuals in the training set (G2 individuals). Therefore, training set composition is not likely to be hindering prediction ability substantially. The use of selfed families to validate parental GBLUPs could also introduce some limits on our prediction ability due to inbreeding depression. Self-fertilization typically delays flowering time in maize (Flint-Garcia et al. 2009; Li et al. 2018), and variation in inbreeding effects on flowering time could introduce additional variation in the progenies that does not exist in the parents (Samayoa et al. 2021). Nevertheless, because variation in inbreeding depression for flowering time is relatively low compared to additive genetic variation (Samayoa et al. 2021), we speculate that inbreeding depression was unlikely to have strongly affected prediction ability.

The most likely factor limiting prediction ability was the heritability within the training set (see equations (1) and (2)). The genomic selection models were trained on individual plants in 2019, whereas the evaluations were conducted on $S_{0.1}$ families
<table><caption>**Table 2.** Variance component ($\hat{\sigma}^2$) and prediction ability estimates from GS model 2 with additive genetic effects ($r_A$) and GS Model 4 with additive and dominance genetic effects ($r_{A+D}$) based on G2 training data on individual plants grown in summer, 2019 in DE and NC.</caption><thead><tr><th rowspan="2">Parameter estimate</th><th colspan="4">G2-DE</th><th colspan="4">G2-NC</th></tr><tr><th>DTA</th><th>DTS</th><th>Plant height</th><th>Ear height</th><th>DTA</th><th>DTS</th><th>Plant height</th><th>Ear height</th></tr></thead><tbody><tr><td colspan="9"><strong>GS Model 2 with additive genetic effects only</strong></td></tr><tr><td>σ<sup>2</sup><sub>Block</sub></td><td>0.4</td><td>0.6</td><td>26.3</td><td>11.2</td><td>0.2</td><td>1.1</td><td>296.4</td><td>51.5</td></tr><tr><td>σ<sup>2</sup><sub>A</sub></td><td>3.1</td><td>5.7</td><td>236.5</td><td>133.1</td><td>2.5</td><td>3.7</td><td>128.5</td><td>81.2</td></tr><tr><td>σ<sup>2</sup><sub>e</sub></td><td>6.4</td><td>12.4</td><td>251.6</td><td>159.0</td><td>6.3</td><td>10.0</td><td>285.3</td><td>153.9</td></tr><tr><td>Prediction ability (r<sub>A</sub>)</td><td>0.40</td><td>0.36</td><td>0.47</td><td>0.50</td><td>0.30</td><td>0.35</td><td>0.42</td><td>0.35</td></tr><tr><td colspan="9"><strong>GS Model 4 with additive and dominance genetic effects</strong></td></tr><tr><td>σ<sup>2</sup><sub>Block</sub></td><td>0.4</td><td>0.6</td><td>26.4</td><td>11.4</td><td>0.2</td><td>1.1</td><td>295.5</td><td>51.4</td></tr><tr><td>σ<sup>2</sup><sub>A</sub></td><td>3.1</td><td>5.7</td><td>238.2</td><td>130.0</td><td>2.5</td><td>3.7</td><td>126.6</td><td>79.3</td></tr><tr><td>σ<sup>2</sup><sub>D</sub></td><td>0</td><td>0.2</td><td>30.2</td><td>16.0</td><td>0.1</td><td>0.4</td><td>31.7</td><td>15.6</td></tr><tr><td>σ<sup>2</sup><sub>e</sub></td><td>6.4</td><td>12.2</td><td>220.0</td><td>144.4</td><td>6.3</td><td>9.6</td><td>255.0</td><td>139.3</td></tr><tr><td>Prediction ability (r<sub>A+D</sub>)</td><td>0.40</td><td>0.36</td><td>0.47</td><td>0.50</td><td>0.30</td><td>0.35</td><td>0.42</td><td>0.35</td></tr></tbody></table>

Traits measured include days to anthesis (DTA), days to silking (DTS), plant height, and ear height.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Prediction abilities for days to anthesis (DTA) and days to silk (DTS) in the G3-DE and G3-NC populations based on models using random subsets of 10 to 100% of the full set of markers.</figcaption></figure>

between 90 and 100% of the full marker set was never greater than 0.01. These results indicate that prediction abilities are not limited by the available number of markers. Curiously, prediction abilities for both populations using only 10% of the total markers remained greater than $r=0.20$. This suggests that the disequilibrium between markers and causal alleles was sufficiently extensive in these populations so that relatively few markers can lead to a substantial positive prediction ability. Direct estimation of the linkage disequilibrium within each chromosome showed that it extends over Mbp distances in both populations (Supplementary Figs. 6, 7 and 8). Although most of the founder linkage disequilibrium was reduced during development of the TropicS-G0 base population, some large blocks remained in the population (Weldekidan et al. 2022). Moreover, the current study was initiated with populations that had undergone two generations of selection, which is also expected to increase linkage disequilibrium.

Genetic distance between training and testing sets was quite small, as indicated above by the values of genomic differentiation

near zero, and as expected for testing sets comprising selfed grandchildren (G3 $S_{0.1}$ families) of individuals in the training set (G2 individuals). Therefore, training set composition is not likely to be hindering prediction ability substantially. The use of selfed families to validate parental GBLUPs could also introduce some limits on our prediction ability due to inbreeding depression. Self-fertilization typically delays flowering time in maize (Flint-Garcia et al. 2009; Li et al. 2018), and variation in inbreeding effects on flowering time could introduce additional variation in the progenies that does not exist in the parents (Samayoa et al. 2021). Nevertheless, because variation in inbreeding depression for flowering time is relatively low compared to additive genetic variation (Samayoa et al. 2021), we speculate that inbreeding depression was unlikely to have strongly affected prediction ability.

The most likely factor limiting prediction ability was the heritability within the training set (see equations (1) and (2)). The genomic selection models were trained on individual plants in 2019, whereas the evaluations were conducted on $S_{0.1}$ families
<table><caption>**Table 3.** Genomic heritability ($h^2$) estimates for DTA and DTS in G2-DE and G2-NC S0 individual plants, and for G3-DE and G3-NC S1 lines, and mean reliability (rel) estimates for G3-DE and G3-NC S0 individual plants.</caption><thead><tr><th rowspan="2">Generation</th><th colspan="2">DE</th><th colspan="2">NC</th></tr><tr><th>DTA</th><th>DTS</th><th>DTA</th><th>DTS</th></tr></thead><tbody><tr><td>h<sup>2</sup> G2 S0s</td><td>0.47</td><td>0.47</td><td>0.46</td><td>0.39</td></tr><tr><td>rel G3 S0s</td><td>0.30</td><td>0.31</td><td>0.31</td><td>0.24</td></tr><tr><td>h<sup>2</sup> G3 S1s</td><td>0.76</td><td>0.75</td><td>0.69</td><td>0.77</td></tr></tbody></table>

G2 generations were grown in summer 2019 and heritability estimates correspond to individual S0 plants with both phenotype and genotype data. G3 S0s were grown in winter 2019–2020 and mean reliability estimates correspond to GBLUPs of individual S0 plants with genotype but no phenotype data. G3 S1s were grown in replicated trials in 2020 and heritability estimates correspond to S0:1 family means.

composed of 50 plants replicated within plots and across two field blocks. The heritability for DTA based on the G3-DE family means in DE was 0.76, whereas the genomic heritability based on the G2-DE grandparental individuals in the training set was 0.47 (Table 3). Similarly, the heritability for DTA based on G3-NC family means was 0.69 in NC, whereas the genomic heritability of the G2-NC grandparental training set was 0.46 (Table 3). Heritability of family means is expected to be substantially greater than heritability of individual plants because the phenotypic variance of means is reduced by averaging across the effects of random micro-environmental influences in the field. Thus, substantial increases in the genomic prediction ability would require improved trait heritabilities in the training set. Unfortunately, this would require at least one additional mating generation to permit progeny testing in replicated trials to increase heritability. This additional investment in time for training might be advantageous if the genomic prediction models maintained effectiveness over multiple generations of selection before requiring retraining, and if multiple generations of genomic selection can be executed in the off-season.

## Genetic gain comparisons for phenotypic vs genomic selection

Our results permit comparison of genomic selection conducted in a winter nursery to phenotypic selection conducted either in target environments or an off-season nursery. First, to estimate the genetic gain from genomic selection in the off-season nursery, we compared the mean flowering times of the lines predicted to be in the earliest 5% tail of the population distribution based on initial genomic selection models to the mean of the randomly sampled lines within each population. This value was compared to the mean difference between the population bulk means of G3 and G2 (which differ due to a single generation of phenotypic selection in the target nurseries). In the DE population, phenotypic selection resulted in gains of about -2 d for anthesis and silk dates; genomic selection in the off-season nursery achieved relative gains of about 50 to 60% (Table 4). In the NC population, phenotypic selection gains were closer to -1 d and genomic selection achieved gains of about 100 to 150% relative to phenotypic selection gains (Table 4). Our estimates of the response to genomic selection may be biased downward, as the genomic selections for extreme tails were based on initial models that did not incorporate the adjustments for field blocks and were implemented in much smaller population samples than the phenotypic selections. These results suggest that genomic selection in the off-season nursery could boost annual gains in selection response by at least

<table><caption>**Table 4.** Estimated gains from phenotypic selection in TropicS G2 in summer, 2019, and evaluated in summer, 2020 in the same locations.</caption><thead><tr><th rowspan="2"></th><th colspan="2">DE</th><th colspan="2">NC</th></tr><tr><th>DTA</th><th>DTS</th><th>DTA</th><th>DTS</th></tr></thead><tbody><tr><td>PS</td><td>-2.3 d</td><td>-1.9 d</td><td>-1.2 d</td><td>-1.4 d</td></tr><tr><td>GS</td><td>-1.2 d</td><td>-1.1 d</td><td>-1.7 d</td><td>-1.6 d</td></tr><tr><td>GS/PS</td><td>52%</td><td>58%</td><td>140%</td><td>114%</td></tr></tbody></table>

Estimated gains from genomic selection in TropicS G3 in winter, 2019–2020, based on prediction models trained on summer 2019 phenotypes and evaluated in summer 2020 in the same locations. Genetic gains were measured in days to anthesis (DTA) and days to silk (DTS) and as proportions of phenotypic selection gains.

50%, and likely closer to 100%, thus reducing the time required to change the population mean flowering time to an acceptably adapted state.

Second, to compare the relative effectiveness of genomic and phenotypic selection for earlier flowering in the off-season nursery, we estimated the correlation between the pollination dates for G3 individuals recorded in the off-season nursery to flowering times of their S<sub>0.1</sub> progeny family BLUEs in target temperate environments (Supplementary Figs. 9 and 10) and compared it to genomic prediction ability. The phenotypic correlations were lower by 0.03 to 0.12 than the genomic prediction ability. This is likely due to the suppression of flowering time photoperiod response and reduced genetic variation that is observed in the tropical off-season environment (Teixeira et al. 2015; Choquette et al. 2023). The positive phenotypic correlations between off-season and target location environments observed here may not occur for all off-season nurseries. Previous evaluations of this same population using the same locations demonstrated that selection for flowering time in the off-season nursery had a very small impact on the phenotypic response measured in the DE and NC target environments (Choquette et al. 2023). Therefore, our specific results likely underestimate the relative effectiveness of genomic to phenotypic selection in the off-season nursery.

For the current study, genomic selection models were trained on individuals evaluated in the target environment where photoperiod sensitivity is induced. These models should select on photoperiod sensitivity in addition to early flowering time per se, whereas phenotypic selection in off-season nurseries with short daylengths acts only on flowering time per se (Teixeira et al. 2015; Choquette et al. 2023). However, in this population, photoperiod sensitivity was dramatically reduced after the first two generations of phenotypic selection in the target environments before we trained genomic selection models (Choquette et al. 2023). It is likely that implementing genomic selection in the off-season from the beginning of the adaptation breeding program would be more advantageous and lead to relatively greater gains in selection for photoperiod sensitivity compared to off-season nursery selection. In the off-season nursery, phenotypic selection is more effective on shortening the basic vegetative phase rather than photoperiod sensitivity. With a better understanding of prediction ability between summer and winter phenotypes, breeders could utilize both genomic selection models and assign weights accordingly to improve prediction models. Minimally, breeders should initially use genomic selection in the off-season nursery in addition to direct phenotypic selection in target environments when breeding for adaptation of a tropical population to temperate environments.
<table><caption>**Table 3.** Genomic heritability ($h^2$) estimates for DTA and DTS in G2-DE and G2-NC S0 individual plants, and for G3-DE and G3-NC S1 lines, and mean reliability (rel) estimates for G3-DE and G3-NC S0 individual plants.</caption><thead><tr><th rowspan="2">Generation</th><th colspan="2">DE</th><th colspan="2">NC</th></tr><tr><th>DTA</th><th>DTS</th><th>DTA</th><th>DTS</th></tr></thead><tbody><tr><td>h<sup>2</sup> G2 S0s</td><td>0.47</td><td>0.47</td><td>0.46</td><td>0.39</td></tr><tr><td>rel G3 S0s</td><td>0.30</td><td>0.31</td><td>0.31</td><td>0.24</td></tr><tr><td>h<sup>2</sup> G3 S1s</td><td>0.76</td><td>0.75</td><td>0.69</td><td>0.77</td></tr></tbody></table>

G2 generations were grown in summer 2019 and heritability estimates correspond to individual S0 plants with both phenotype and genotype data. G3 S0s were grown in winter 2019–2020 and mean reliability estimates correspond to GBLUPs of individual S0 plants with genotype but no phenotype data. G3 S1s were grown in replicated trials in 2020 and heritability estimates correspond to S0:1 family means.

composed of 50 plants replicated within plots and across two field blocks. The heritability for DTA based on the G3-DE family means in DE was 0.76, whereas the genomic heritability based on the G2-DE grandparental individuals in the training set was 0.47 (Table 3). Similarly, the heritability for DTA based on G3-NC family means was 0.69 in NC, whereas the genomic heritability of the G2-NC grandparental training set was 0.46 (Table 3). Heritability of family means is expected to be substantially greater than heritability of individual plants because the phenotypic variance of means is reduced by averaging across the effects of random micro-environmental influences in the field. Thus, substantial increases in the genomic prediction ability would require improved trait heritabilities in the training set. Unfortunately, this would require at least one additional mating generation to permit progeny testing in replicated trials to increase heritability. This additional investment in time for training might be advantageous if the genomic prediction models maintained effectiveness over multiple generations of selection before requiring retraining, and if multiple generations of genomic selection can be executed in the off-season.

## Genetic gain comparisons for phenotypic vs genomic selection

Our results permit comparison of genomic selection conducted in a winter nursery to phenotypic selection conducted either in target environments or an off-season nursery. First, to estimate the genetic gain from genomic selection in the off-season nursery, we compared the mean flowering times of the lines predicted to be in the earliest 5% tail of the population distribution based on initial genomic selection models to the mean of the randomly sampled lines within each population. This value was compared to the mean difference between the population bulk means of G3 and G2 (which differ due to a single generation of phenotypic selection in the target nurseries). In the DE population, phenotypic selection resulted in gains of about -2 d for anthesis and silk dates; genomic selection in the off-season nursery achieved relative gains of about 50 to 60% (Table 4). In the NC population, phenotypic selection gains were closer to -1 d and genomic selection achieved gains of about 100 to 150% relative to phenotypic selection gains (Table 4). Our estimates of the response to genomic selection may be biased downward, as the genomic selections for extreme tails were based on initial models that did not incorporate the adjustments for field blocks and were implemented in much smaller population samples than the phenotypic selections. These results suggest that genomic selection in the off-season nursery could boost annual gains in selection response by at least

<table><caption>**Table 4.** Estimated gains from phenotypic selection in TropicS G2 in summer, 2019, and evaluated in summer, 2020 in the same locations.</caption><thead><tr><th rowspan="2"></th><th colspan="2">DE</th><th colspan="2">NC</th></tr><tr><th>DTA</th><th>DTS</th><th>DTA</th><th>DTS</th></tr></thead><tbody><tr><td>PS</td><td>-2.3 d</td><td>-1.9 d</td><td>-1.2 d</td><td>-1.4 d</td></tr><tr><td>GS</td><td>-1.2 d</td><td>-1.1 d</td><td>-1.7 d</td><td>-1.6 d</td></tr><tr><td>GS/PS</td><td>52%</td><td>58%</td><td>140%</td><td>114%</td></tr></tbody></table>

Estimated gains from genomic selection in TropicS G3 in winter, 2019–2020, based on prediction models trained on summer 2019 phenotypes and evaluated in summer 2020 in the same locations. Genetic gains were measured in days to anthesis (DTA) and days to silk (DTS) and as proportions of phenotypic selection gains.

50%, and likely closer to 100%, thus reducing the time required to change the population mean flowering time to an acceptably adapted state.

Second, to compare the relative effectiveness of genomic and phenotypic selection for earlier flowering in the off-season nursery, we estimated the correlation between the pollination dates for G3 individuals recorded in the off-season nursery to flowering times of their S<sub>0.1</sub> progeny family BLUEs in target temperate environments (Supplementary Figs. 9 and 10) and compared it to genomic prediction ability. The phenotypic correlations were lower by 0.03 to 0.12 than the genomic prediction ability. This is likely due to the suppression of flowering time photoperiod response and reduced genetic variation that is observed in the tropical off-season environment (Teixeira et al. 2015; Choquette et al. 2023). The positive phenotypic correlations between off-season and target location environments observed here may not occur for all off-season nurseries. Previous evaluations of this same population using the same locations demonstrated that selection for flowering time in the off-season nursery had a very small impact on the phenotypic response measured in the DE and NC target environments (Choquette et al. 2023). Therefore, our specific results likely underestimate the relative effectiveness of genomic to phenotypic selection in the off-season nursery.

For the current study, genomic selection models were trained on individuals evaluated in the target environment where photoperiod sensitivity is induced. These models should select on photoperiod sensitivity in addition to early flowering time per se, whereas phenotypic selection in off-season nurseries with short daylengths acts only on flowering time per se (Teixeira et al. 2015; Choquette et al. 2023). However, in this population, photoperiod sensitivity was dramatically reduced after the first two generations of phenotypic selection in the target environments before we trained genomic selection models (Choquette et al. 2023). It is likely that implementing genomic selection in the off-season from the beginning of the adaptation breeding program would be more advantageous and lead to relatively greater gains in selection for photoperiod sensitivity compared to off-season nursery selection. In the off-season nursery, phenotypic selection is more effective on shortening the basic vegetative phase rather than photoperiod sensitivity. With a better understanding of prediction ability between summer and winter phenotypes, breeders could utilize both genomic selection models and assign weights accordingly to improve prediction models. Minimally, breeders should initially use genomic selection in the off-season nursery in addition to direct phenotypic selection in target environments when breeding for adaptation of a tropical population to temperate environments.
## Data availability

Data and R scripts are available on figshare at https://doi.org/10.25387/g3.22191277 for access along with the paper.

Supplemental material available at G3 online.

## Acknowledgements

We are grateful for the staff and students who helped with this project. We would especially like to thank Eric Butoto, Danielle Darden, Kendra Jones, Will Pequigney, Anna Rogers, Ryan Sartor, and Will Voelker for aiding in tissue sampling and conducting many pollinations as well as Greg Marshall and Shannon Sermons for assisting with harvest and shelling.

## Funding and competing interests

This work was supported by the Agriculture and Food Research Initiative, USDA National Institute of Food and Agriculture, Grant nos. 2011-67003-30342 and 2019-67013-29170, in addition to INRAE through the project EXPOSE (PAF_18). The research funders had no role in the research design, execution, analysis, interpretation, and reporting.

## Literature cited

Amadeu RR, Cellon C, Olmstead JW, Garcia AAF, Resende MFR, et al. AGHmatrix: R package to construct relationship matrices for autotetraploid and diploid species: a blueberry example. Plant Genome. 2016;9(3):1-10. doi:10.3835/plantgenome2016.01.0009.

Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, et al. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633-2635. doi:10.1093/bioinformatics/btm308.

Buckler ES, Holland JB, Bradbury PJ, Acharya CB, Brown PJ, et al. The genetic architecture of maize flowering time. Science. 2009; 325(5941):714-718. doi:10.1126/science.1174276.

Butler DG, Cullis BR, Gilmour AR, Gogel BG, Thompson R. ASReml-R Reference Manual Version 4. Hemel Hempstead (UK): VSN International Ltd; 2017.

Chardon F, Virlon B, Moreau L, Falque M, Joets J, et al. Genetic architecture of flowering time in maize as inferred from quantitative trait loci meta-analysis and synteny conservation with the rice genome. Genetics. 2004;168(4):2169-2185. doi:10.1534/genetics.104.032375.

Choquette NE, Holland JB, Weldekidan T, Drouault J, de Leon N, et al. Environment-specific selection alters flowering-time plasticity and results in pervasive pleiotropic responses in maize. New Phytol. 2023;238(2):737-749. doi:10.1111/nph.18769.

Cortés AJ, López-Hernández F. Harnessing crop wild diversity for climate change adaptation. Genes (Basel). 2021;12(5):783. doi:10.3390/genes12050783.

Crossa J, Pérez-Rodríguez P, Cuevas J, Montesinos-López O, Jarquín D, et al. Genomic selection in plant breeding: methods, models, and perspectives. Trends Plant Sci. 2017;22(11):961-975. doi:10.1016/j.trends.2017.08.011.

de los Campos G, Hickey JM, Pong-Wong R, Daetwyler HD, Calus MPL. Whole-genome regression and prediction methods applied to plant and animal breeding. Genetics. 2013;193(2):327-345. doi: 10.1534/genetics.112.143313.

Duvick DN. The contribution of breeding to yield advances in maize (Zea mays L.). Adv Agron. 2005;86:83-145. doi:10.1016/S0065-2113(05)86002-X.

Falconer DS, Mackay TFC. Introduction to Quantitative Genetics. Harlow (UK): Longman; 1996.

Flint-Garcia SA, Buckler ES, Tiffin P, Ersoz E, Springer NM. Heterosis is prevalent for multiple traits in diverse maize germplasm. PLoS One. 2009;4(10):e7433. doi:10.1371/journal.pone.0007433.

Frey TJ, Weldekidan T, Colbert T, Wolters PJCC, Hawk JA. Fitness evaluation of Rcg1, a locus that confers resistance to Colletotrichum graminicola (Ces.) G.W. Wils. Using near-isogenic maize hybrids. Crop Sci. 2011;51(4):1551-1563. doi:10.2135/cropsci2010.10.0613.

Goodman MM. Developing temperate inbreds using tropical maize germplasm: rationale, results, conclusions. Maydica. 2004;49(3): 209-219.

Goudet J, Jombart T. hierfstat: estimation and tests of hierarchical F-statistics. R Packag. version 0.04-22, 2015.

Gouesnard B, Rebourg C, Welcker C, Charcosset A. Analysis of photo-period sensitivity within a collection of tropical maize populations. Genet Resour Crop Evol. 2002;49(5):471-481. doi:10.1023/A:1020982827604.

Heffner EL, Sorrells ME, Jannink J-L. Genomic selection for crop improvement. Crop Sci. 2009;49(1):1-12. doi:10.2135/cropsci2008.08.0512.

Holland JB, Goodman MM, Castillo-Gonzalez F. Identification of agronomically superior latin American maize accessions via multi-stage evaluations. Crop Sci. 1996;36(3):778-784. doi:10.2135/cropsci1996.0011183X003600030041x.

Hufford MB, Xu X, Van Heerwaarden J, Pyhäjärvi T, Chia JM, et al. Comparative population genomics of maize domestication and improvement. Nat Genet. 2012;44(7):808-811. doi:10.1038/ng.2309.

Isik F, Holland J, Maltecca C. Genetic Data Analysis for Plant and Animal Breeding. New York (NY): Springer; 2017.

Jannink J-L, Lorenz AJ, Iwata H. Genomic selection in plant breeding: from theory to practice. Brief Funct Genomics. 2010;9(2):166-177. doi:10.1093/bfgp/elq001.

Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24(11):1403-1405. doi:10.1093/bioinformatics/btn129.

Li Y, Willer CJ, Ding J, Scheet P, Abecasis GR. MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes. Genet Epidemiol. 2010;34(8):816-834. doi:10.1002/GEPI.20533.

Li Z, Coffey L, Garfin J, Miller ND, White MR, et al. Genotype-by-environment interactions affecting heterosis in maize, (L. Lukens, ed.). PLoS One. 2018;13(1):e0191321. doi:10.1371/journal.pone.0191321.

Liu K, Goodman M, Muse S, Smith JS, Buckler E, et al. Genetic structure and diversity among maize inbred lines as inferred from DNA microsatellites. Genetics. 2003;165(4):2117-2128. doi:10.1093/genetics/165.4.2117.

Lorenz AJ, Chao S, Asoro FGFG, Heffner EL, Hayashi T, et al. Genomic selection in plant breeding: knowledge and prospects. Adv Agron. 2011;110:77-123. doi:10.1016/B978-0-12-385531-2.00002-5.

Lorenzana RE, Bernardo R. Accuracy of genotypic value predictions for marker-based selection in biparental plant populations. Theor Appl Genet. 2009;120:151-161. doi:10.1007/s00122-009-1166-3.

Manching H, Sengupta S, Hopper KR, Polson SW, Ji Y, et al. Phased genotyping-by-sequencing enhances analysis of genetic diversity and reveals divergent copy number variants in maize. G3 (Bethesda). 2017;7(7):2161-2170. doi:10.1534/g3.117.042036.

Meuwissen THE, Hayes BJ, Goddard ME. Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2001; 157(4):1819-1829. doi:10.1093/genetics/157.4.1819.
## Data availability

Data and R scripts are available on figshare at https://doi.org/10.25387/g3.22191277 for access along with the paper.

Supplemental material available at G3 online.

## Acknowledgements

We are grateful for the staff and students who helped with this project. We would especially like to thank Eric Butoto, Danielle Darden, Kendra Jones, Will Pequigney, Anna Rogers, Ryan Sartor, and Will Voelker for aiding in tissue sampling and conducting many pollinations as well as Greg Marshall and Shannon Sermons for assisting with harvest and shelling.

## Funding and competing interests

This work was supported by the Agriculture and Food Research Initiative, USDA National Institute of Food and Agriculture, Grant nos. 2011-67003-30342 and 2019-67013-29170, in addition to INRAE through the project EXPOSE (PAF_18). The research funders had no role in the research design, execution, analysis, interpretation, and reporting.

## Literature cited

Amadeu RR, Cellon C, Olmstead JW, Garcia AAF, Resende MFR, et al. AGHmatrix: R package to construct relationship matrices for autotetraploid and diploid species: a blueberry example. Plant Genome. 2016;9(3):1-10. doi:10.3835/plantgenome2016.01.0009.

Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, et al. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633-2635. doi:10.1093/bioinformatics/btm308.

Buckler ES, Holland JB, Bradbury PJ, Acharya CB, Brown PJ, et al. The genetic architecture of maize flowering time. Science. 2009; 325(5941):714-718. doi:10.1126/science.1174276.

Butler DG, Cullis BR, Gilmour AR, Gogel BG, Thompson R. ASReml-R Reference Manual Version 4. Hemel Hempstead (UK): VSN International Ltd; 2017.

Chardon F, Virlon B, Moreau L, Falque M, Joets J, et al. Genetic architecture of flowering time in maize as inferred from quantitative trait loci meta-analysis and synteny conservation with the rice genome. Genetics. 2004;168(4):2169-2185. doi:10.1534/genetics.104.032375.

Choquette NE, Holland JB, Weldekidan T, Drouault J, de Leon N, et al. Environment-specific selection alters flowering-time plasticity and results in pervasive pleiotropic responses in maize. New Phytol. 2023;238(2):737-749. doi:10.1111/nph.18769.

Cortés AJ, López-Hernández F. Harnessing crop wild diversity for climate change adaptation. Genes (Basel). 2021;12(5):783. doi:10.3390/genes12050783.

Crossa J, Pérez-Rodríguez P, Cuevas J, Montesinos-López O, Jarquín D, et al. Genomic selection in plant breeding: methods, models, and perspectives. Trends Plant Sci. 2017;22(11):961-975. doi:10.1016/j.trends.2017.08.011.

de los Campos G, Hickey JM, Pong-Wong R, Daetwyler HD, Calus MPL. Whole-genome regression and prediction methods applied to plant and animal breeding. Genetics. 2013;193(2):327-345. doi: 10.1534/genetics.112.143313.

Duvick DN. The contribution of breeding to yield advances in maize (Zea mays L.). Adv Agron. 2005;86:83-145. doi:10.1016/S0065-2113(05)86002-X.

Falconer DS, Mackay TFC. Introduction to Quantitative Genetics. Harlow (UK): Longman; 1996.

Flint-Garcia SA, Buckler ES, Tiffin P, Ersoz E, Springer NM. Heterosis is prevalent for multiple traits in diverse maize germplasm. PLoS One. 2009;4(10):e7433. doi:10.1371/journal.pone.0007433.

Frey TJ, Weldekidan T, Colbert T, Wolters PJCC, Hawk JA. Fitness evaluation of Rcg1, a locus that confers resistance to Colletotrichum graminicola (Ces.) G.W. Wils. Using near-isogenic maize hybrids. Crop Sci. 2011;51(4):1551-1563. doi:10.2135/cropsci2010.10.0613.

Goodman MM. Developing temperate inbreds using tropical maize germplasm: rationale, results, conclusions. Maydica. 2004;49(3): 209-219.

Goudet J, Jombart T. hierfstat: estimation and tests of hierarchical F-statistics. R Packag. version 0.04-22, 2015.

Gouesnard B, Rebourg C, Welcker C, Charcosset A. Analysis of photo-period sensitivity within a collection of tropical maize populations. Genet Resour Crop Evol. 2002;49(5):471-481. doi:10.1023/A:1020982827604.

Heffner EL, Sorrells ME, Jannink J-L. Genomic selection for crop improvement. Crop Sci. 2009;49(1):1-12. doi:10.2135/cropsci2008.08.0512.

Holland JB, Goodman MM, Castillo-Gonzalez F. Identification of agronomically superior latin American maize accessions via multi-stage evaluations. Crop Sci. 1996;36(3):778-784. doi:10.2135/cropsci1996.0011183X003600030041x.

Hufford MB, Xu X, Van Heerwaarden J, Pyhäjärvi T, Chia JM, et al. Comparative population genomics of maize domestication and improvement. Nat Genet. 2012;44(7):808-811. doi:10.1038/ng.2309.

Isik F, Holland J, Maltecca C. Genetic Data Analysis for Plant and Animal Breeding. New York (NY): Springer; 2017.

Jannink J-L, Lorenz AJ, Iwata H. Genomic selection in plant breeding: from theory to practice. Brief Funct Genomics. 2010;9(2):166-177. doi:10.1093/bfgp/elq001.

Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24(11):1403-1405. doi:10.1093/bioinformatics/btn129.

Li Y, Willer CJ, Ding J, Scheet P, Abecasis GR. MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes. Genet Epidemiol. 2010;34(8):816-834. doi:10.1002/GEPI.20533.

Li Z, Coffey L, Garfin J, Miller ND, White MR, et al. Genotype-by-environment interactions affecting heterosis in maize, (L. Lukens, ed.). PLoS One. 2018;13(1):e0191321. doi:10.1371/journal.pone.0191321.

Liu K, Goodman M, Muse S, Smith JS, Buckler E, et al. Genetic structure and diversity among maize inbred lines as inferred from DNA microsatellites. Genetics. 2003;165(4):2117-2128. doi:10.1093/genetics/165.4.2117.

Lorenz AJ, Chao S, Asoro FGFG, Heffner EL, Hayashi T, et al. Genomic selection in plant breeding: knowledge and prospects. Adv Agron. 2011;110:77-123. doi:10.1016/B978-0-12-385531-2.00002-5.

Lorenzana RE, Bernardo R. Accuracy of genotypic value predictions for marker-based selection in biparental plant populations. Theor Appl Genet. 2009;120:151-161. doi:10.1007/s00122-009-1166-3.

Manching H, Sengupta S, Hopper KR, Polson SW, Ji Y, et al. Phased genotyping-by-sequencing enhances analysis of genetic diversity and reveals divergent copy number variants in maize. G3 (Bethesda). 2017;7(7):2161-2170. doi:10.1534/g3.117.042036.

Meuwissen THE, Hayes BJ, Goddard ME. Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2001; 157(4):1819-1829. doi:10.1093/genetics/157.4.1819.
Mikel MA, Dudley JW. Evolution of North American dent corn from public to proprietary germplasm. Crop Sci. 2006;46(3): 1193-1205. doi:10.2135/cropsci2005.10-0371.

Nelson PT, Goodman MM. Evaluation of elite exotic maize inbreds for use in temperate breeding. Crop Sci. 2008;48(1):85-92. doi:10.2135/cropsci2007.05.0287.

Pérez P, de los Campos G, Goddard ME. Genome-wide regression and prediction with the BGLR statistical package. Genetics. 2014; 198(2):483-495. doi:10.1534/genetics.114.164442.

Pollak LM. The history and success of the public-private project on germplasm enhancement of maize (GEM). Adv Agron. 2003;78: 45-87. doi:10.1016/S0065-2113(02)78002-4.

Ramstein GP, Larsson SJ, Cook JP, Edwards JW, Ersoz ES, et al. Dominance effects and functional enrichments improve prediction of agronomic traits in hybrid maize. Genetics. 2020;215(1): 215-230. doi:10.1534/genetics.120.303025.

Rogers AR, Bian Y, Krakowsky M, Peters D, Turnbull C, et al. Genomic prediction for the germplasm enhancement of maize project. Plant Genome. 2022;15(4):e20267. doi:10.1002/tpg2.20267.

Romay MC, Millard MJ, Glaubitz JC, Peiffer JA, Swarts KL, et al. Comprehensive genotyping of the USA national maize inbred seed bank. Genome Biol. 2013;14:R55. doi:10.1186/gb-2013-14-6-r55.

Salvi S, Sponza G, Morgante M, Torres D, Niu X, et al. Conserved non-coding genomic sequences associated with a flowering-time quantitative trait locus in maize. Proc Natl Acad Sci U S A. 2007; 104(27):11376-11381. doi:10.1073/pnas.0704145104.

Samayoa LF, Olukolu BA, Yang CJ, Chen Q, Stetter MG, et al. Domestication reshaped the genetic basis of inbreeding depression in a maize landrace compared to its wild relative, teosinte. PLoS Genet. 2021;17(12):e1009797. doi:10.1371/journal.pgen.1009797.

Schmidt P, Hartung J, Bennewitz J, Hans-Peter P. Heritability in plant breeding on a genotype-difference basis. Genetics. 2019;212(4): 991-1008. doi:10.1534/genetics.119.302134.

Schwarz G. Estimating the dimension of a model. Ann Stat. 1978;6(2): 461-464. doi:10.1214/aos/1176344136.

Smith JS, Trevisan W, McCunn A, Huffman WE. Global dependence upon corn belt dent maize germplasm: challenges and opportunities. Crop Sci. 2022;62(6):2039-2066. doi:10.1002/csc2.20802.

Teixeira JEC, Weldekidan T, De Leon N, Flint-Garcia S, Holland JB, et al. Hallauer's tusón: a decade of selection for tropical-to-temperate phenological adaptation in maize. Heredity (Edinb). 2015;114(2):229-240. doi:10.1038/hdy.2014.90.

VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91(11):4414-4423. doi:10.3168/jds.2007-0980.

Vitezica ZG, Varona L, Legarra A. On the additive and dominant variance and covariance of individuals within the genomic selection scope. Genetics. 2013;195(4):1223-1230. doi:10.1534/genetics.113.155176.

Weldekidan T, Manching H, Choquette N, de Leon N, Flint-Garcia S, et al. Registration of tropical populations of maize selected in parallel for early flowering time across the United States. J Plant Regist. 2022;16(1):100-108. doi:10.1002/plr2.20181.

White MR, Mikel MA, De Leon N, Kaepler SM. Diversity and heterotic patterns in North American proprietary dent maize germplasm. Crop Sci. 2020;60(1):100-114. doi:10.1002/csc2.20050.

Wisser RJ, Fang Z, Holland JB, Teixeira JEC, Dougherty J, et al. The genomic basis for short-term evolution of environmental adaptation in maize. Genetics. 2019;213(4):1479-1494. doi:10.1534/genetics.119.302780.

Yu J, Zhang Z, Zhu C, Tabanao DA, Pressoir G, et al. Simulation appraisal of the adequacy of number of background markers for relationship estimation in association mapping. Plant Genome. 2009;2(1):63-77. doi:10.3835/plantgenome2008.09.0009.

Editor: A. Lipka
Mikel MA, Dudley JW. Evolution of North American dent corn from public to proprietary germplasm. Crop Sci. 2006;46(3): 1193-1205. doi:10.2135/cropsci2005.10-0371.

Nelson PT, Goodman MM. Evaluation of elite exotic maize inbreds for use in temperate breeding. Crop Sci. 2008;48(1):85-92. doi:10.2135/cropsci2007.05.0287.

Pérez P, de los Campos G, Goddard ME. Genome-wide regression and prediction with the BGLR statistical package. Genetics. 2014; 198(2):483-495. doi:10.1534/genetics.114.164442.

Pollak LM. The history and success of the public-private project on germplasm enhancement of maize (GEM). Adv Agron. 2003;78: 45-87. doi:10.1016/S0065-2113(02)78002-4.

Ramstein GP, Larsson SJ, Cook JP, Edwards JW, Ersoz ES, et al. Dominance effects and functional enrichments improve prediction of agronomic traits in hybrid maize. Genetics. 2020;215(1): 215-230. doi:10.1534/genetics.120.303025.

Rogers AR, Bian Y, Krakowsky M, Peters D, Turnbull C, et al. Genomic prediction for the germplasm enhancement of maize project. Plant Genome. 2022;15(4):e20267. doi:10.1002/tpg2.20267.

Romay MC, Millard MJ, Glaubitz JC, Peiffer JA, Swarts KL, et al. Comprehensive genotyping of the USA national maize inbred seed bank. Genome Biol. 2013;14:R55. doi:10.1186/gb-2013-14-6-r55.

Salvi S, Sponza G, Morgante M, Tomes D, Niu X, et al. Conserved non-coding genomic sequences associated with a flowering-time quantitative trait locus in maize. Proc Natl Acad Sci U S A. 2007; 104(27):11376-11381. doi:10.1073/pnas.0704145104.

Samayoa LF, Olukolu BA, Yang CJ, Chen Q, Stetter MG, et al. Domestication reshaped the genetic basis of inbreeding depression in a maize landrace compared to its wild relative, teosinte. PLoS Genet. 2021;17(12):e1009797. doi:10.1371/journal.pgen.1009797.

Schmidt P, Hartung J, Bennewitz J, Hans-Peter P. Heritability in plant breeding on a genotype-difference basis. Genetics. 2019;212(4): 991-1008. doi:10.1534/genetics.119.302134.

Schwarz G. Estimating the dimension of a model. Ann Stat. 1978;6(2): 461-464. doi:10.1214/aos/1176344136.

Smith JS, Trevisan W, McCunn A, Huffman WE. Global dependence upon corn belt dent maize germplasm: challenges and opportunities. Crop Sci. 2022;62(6):2039-2066. doi:10.1002/csc2.20802.

Teixeira JEC, Weldekidan T, De Leon N, Flint-Garcia S, Holland JB, et al. Hallauer's tusón: a decade of selection for tropical-to-temperate phenological adaptation in maize. Heredity (Edinb). 2015;114(2):229-240. doi:10.1038/hdy.2014.90.

VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91(11):4414-4423. doi:10.3168/jds.2007-0980.

Vitezica ZG, Varona L, Legarra A. On the additive and dominant variance and covariance of individuals within the genomic selection scope. Genetics. 2013;195(4):1223-1230. doi:10.1534/genetics.113.155176.

Weldekidan T, Manching H, Choquette N, de Leon N, Flint-Garcia S, et al. Registration of tropical populations of maize selected in parallel for early flowering time across the United States. J Plant Regist. 2022;16(1):100-108. doi:10.1002/plr2.20181.

White MR, Mikel MA, De Leon N, Kaepler SM. Diversity and heterotic patterns in North American proprietary dent maize germplasm. Crop Sci. 2020;60(1):100-114. doi:10.1002/csc2.20050.

Wisser RJ, Fang Z, Holland JB, Teixeira JEC, Dougherty J, et al. The genomic basis for short-term evolution of environmental adaptation in maize. Genetics. 2019;213(4):1479-1494. doi:10.1534/genetics.119.302780.

Yu J, Zhang Z, Zhu C, Tabanao DA, Pressoir G, et al. Simulation appraisal of the adequacy of number of background markers for relationship estimation in association mapping. Plant Genome. 2009;2(1):63-77. doi:10.3835/plantgenome2008.09.0009.

Editor: A. Lipka