The goal of wheat breeding is the development of superior cultivars tailored to specific environments, and the identification of promising crosses is crucial for the success of breeding programs. Although genomic estimated breeding values were developed to estimate additive effects of genotypes before testing as parents, application has focused on predicting performance of candidate lines, ignoring non-additive genetic effects. However, non-additive genetic effects are hypothesized to be especially importance in allopolyploid species due to the interaction between homeologous genes. The objectives of this study were to model additive and additive-by-additive epistatic effects to better delineate the genetic architecture of grain yield in wheat and to the improve accuracy of genomewide predictions. The dataset utilized consisted of 3740 F5:6experimental lines tested in the K-State wheat breeding program across the years 2016 and 2018. Covariance matrices were calculated based on whole and sub-genome marker data and the natural and orthogonal interaction approach (NOIA) was used to estimate variance components for additive and additive-by-additive epistatic effects. Incorporating epistatic effects in additive models resulted in non-orthogonal partitioning of genetic effects but increased total genetic variance and reduced DIC. Estimation of sub-genome effects indicated that genotypes with the greatest whole genome effects often combine sub-genomes with intermediate to high effects, suggesting potential for crossing parental lines which have complementary sub-genome effects. Modeling epistasis in either whole-genome or sub-genome models led to a marginal (3%) but significant improvement in genomic prediction accuracy, which could result in significant genetic gains across multiple cycles of breeding.