Genome-wide QTL mapping of yield and agronomic traits in two widely adapted winter wheat cultivars from multiple mega-environments

Plant materials, field trials and phenotyping

A population of 124 F_5:7 RILs derived from the cross ‘TAM 112’/‘TAM 111’ was used to map grain yield and related agronomic traits. Both TAM 111 and TAM 112 were developed by Texas A&M AgriLife Research (Lazar et al., 2004; Rudd et al., 2014) and well-adapted hard red winter wheat (HRWW) cultivars in the Southern Great Plains of the United States. TAM 112 showed better adaption to low input environments, whereas TAM 111 typically produced higher yield in moderate to high input environments. The RILs along with the two parents were planted in five locations including Bushland, TX (35°06′N, 102°27′W), Chillicothe, TX (34°15′N, 99°30′W), Clovis, NM (34°24′N, 103°12′W), Etter, TX (35°51′N, 101°58′W), and Uvalde, TX (29°21′N, 99°75′W) during five years (2011 to 2014, and 2017) (Table 1). Yield data were collected from 28 year-location-management combinations (environments) which included Bushland dryland (BD) in 2011 (11BD), 2012 (12BD) and 2017 (17BD), Bushland irrigated (BI) in 2017 (17BI), Chillicothe dryland (CH) in 2011 (11CH), 2012 (12CH), and 2014 (14CH), Clovis irrigated (CVI) in 2017 (17CVI), Etter with linear irrigation in 2017 (17EI), and various irrigation levels in Etter and Uvalde. Etter with five irrigation levels as 40% of evapotranspiration (ET) demand (EP1), 50% of ET demand (EP2), 65% of ET demand (EP3), 75% of ET demand (EP4), 100% of ET demand (EP5) during 2011–2014 included 14 environments (11EP1, 11EP2, 11EP3, 11EP4, 11EP5, 12EP1, 12EP2, 12EP3, 13EP2, 13EP3, 13EP4, 13EP5, 14EP4, and 14EP5). Uvalde included dryland (UVD) and Uvalde at three irrigation levels at 50% (UV5), 70% (UV7), and 100% of ET demand (UVL) in 2012 and 2013 (12UVD, 12UV5, 12UV7, 12UVL, and 13UVL). Overall, sixteen environments (11BD, 12BD, 17BD, 11CH, 12CH, 14CH, 11EP1, 11EP2, 11EP3, 12EP1, 12EP2, 12EP3, 13EP2, 13EP3, 12UVD, and 12UV5) were considered as dryland condition and twelve others (17BI, 17CVI, 11EP4, 11EP5, 13EP4, 13EP5, 14EP4, 14EP5, 17EI, 12UV7, 12UVL, and 13UVL) were considered as irrigated conditions. All trials were replicated twice in an alpha-lattice design with an incomplete block size of five plots and each parent occurring three times in each replication. The plot dimension was 6.09 m × 1.52 m in the dryland environments and 4.57 m × 1.52 m in the irrigated ones with a 0.3-m space between plots. Standard agronomic practices were performed in each trial (Dhakal et al., 2021; Yang et al., 2020b).

All 28 environments were harvested using a combine harvester and the total plot weight was used to calculate grain yield (YLD). Traits of test weight (TW), heading date (HD), and plant height (HT) were measured in a subset of environments. TW was measured using Seedburo equipment (http://www.seedburo.com, Des Plaines, IL, USA) from 19 environments (11BD, 17BD, 17BI, 12CH, 14CH, 11EP5, 12EP1, 12EP2, 12EP3, 13EP2, 13EP3, 13EP4, 13EP5, 14EP4, 14EP5, 17EP5, 12UV5, 12UV7, and 12UVL). HD was recorded at Feekes growth stage 10.1 when half of the plants were fully visible on heads from 11 environments (11BD, 12BD, 17BD, 11EP1, 11EP2, 11EP3, 11EP4, 11EP5, 12EP1, 12EP2, and 12EP3). Plant height was measured in centimeters (cm) from representative plants in each plot as the distance from the base of the stem to the top of the spike excluding awns at maturity in 11 environments (11BD, 12BD, 17BD, 12CH, 17CVI, 11EP1, 11EP2, 11EP3, 11EP4, 11EP5, and 17EI).

Statistical analysis

Descriptive statistics were calculated using PROC UNIVARIATE in SAS 9.4 (SAS Institute, Cary, NC, USA). The histograms of the residuals for all traits were approximately normal. Best linear unbiased prediction (BLUP) means of lines were calculated for the agronomic traits from the individual environment (IE) as well as multi-environment trials (MET). BLUP was calculated using a restricted maximum likelihood (REML) approach implemented in the ‘lme4’ package through META-R (Alvarado et al., 2018).The analysis of variance (ANOVA) for IE followed a linear statistical model of individual environment analysis with replication and incomplete block as random effects, while for MET linear statistical model of combined environment analysis was run using PROC MIXED in SAS 9.4 with the environment, replication and incomplete block as random effects. The entry-mean heritabilities, pearson’s correlations, biplot analyses, mega-environments classifications followed the same procedures (Dhakal et al., 2021).

Genotyping and linkage map construction

DNA extraction, Illumina Infinium iSelect 90K array SNP genotyping and Genome Studio clustering followed the procedures described in the literature (Assanga et al., 2017; Dhakal et al., 2018; Liu et al., 2016; Yang et al., 2019). This population was also genotyped with double digest restriction-site associated DNA sequencing (ddRADSeq) method developed by Peterson et al. (2012) on an Illumina HiSeq 2500 platform (2  × 125 bp paired-end) following the standard procedures (Yang et al., 2020b). JoinMap v4.0 software (Van Ooijen, 2006) was used to construct the genetic maps using the standard procedures (Dhakal et al., 2021; Yang et al., 2020b). The SNP sequences were listed in Table S1.

QTL analysis

The additive effects, epistasis, additive-by-environment, and epistasis-by-environment interactions were analyzed using the integrated composite interval mapping (ICIM) function implemented in QTL IciMapping software (Meng et al., 2015). The threshold for declaring a QTL was determined through a permutation test (n = 1, 000) for a single environment to obtain a 0.05 genome-wide probability level of Type I error. QTL analyses were run for individual environment (IE) for additive effects (ADD), multiple environment (EPI), and within and across MEs following the standard procedures (Dhakal et al., 2021). QTL designation followed McIntosh et al. (2003) guidelines with a slight change using the format Qtrait.tamu.chrom.Mb, where trait represents a trait name, tamu indicates Texas A&M University, chrom means the chromosome harboring the QTL, and Mb indicates the Mega base pair (Mb) position of the peak SNP within a QTL according to sequence alignment using the IWGSC RefSeq v1.0 reference genome (International Wheat Genome Sequencing Consortium, 2014).

Phenotypic variability, heritability, and correlations

Each environment was unique and different from the others across years. In general, Chillicothe and Uvalde, TX were warmer and Etter, TX located in the northern Texas High Plains tends to be cooler. Bushland received below average rainfall in 2011 and almost double in 2012 and 2017; however, the rainfalls during the wheat growing seasons were low during 2011-2014, which were all considered drought years for wheat (Table 1). Temperature was within the optimum range required for winter wheat growth. Combined ANOVA showed highly significant genotype differences (p < 0.001) among RILs for all the traits (p < 0.01) (Table 2). All the traits were highly heritable (0.77–0.96). Means for YLD, TW, HD, and HT were 287.7 g m^-2, 760.6 kg m⁻³, 115.5 Julian days, and 73.8 cm, respectively. Transgressive segregation was observed for all the traits. Based on the means across all environments, TAM 112 had higher yield while TAM 111 had higher TW, HD and HT based on the overall BLUP means from all tested environments. However, the two parents did not differ significantly for all traits (Table 2).

The phenotypic relationships between grain yield and agronomic traits in individual environment was determined using Pearson’s correlation coefficients (Table S2). In general, negative correlations were found between HD and YLD in the dryland environments with high correlation values at some dryland environments. Most correlations between YLD and HT were positive and significant in the dryland environments, especially in 2011. Most correlations between YLD and TW were positive and significant except for a few dryland environments. HD and HT were significantly and negatively correlated in dryland environments except for 17BD as 2017 was a high rainfall year while they were positively correlated in highly irrigated environments. TW had negative correlations with HD while it had positive correlations with HT in dryland environments but negatively correlated in irrigated environments with limited amount of data (Table S2).

Boxplot, biplot, and mega-environment

Heading date was later for the year 2011 except the 11BD environment. Earlier HD was seen in 11BD, 12BD, and 17BD (Fig. S1). Plant height was severely affected by drought, as shown in the Bushland dryland and less affected under irrigated environments in Etter. Plant height was less affected by the year with high rainfall such as 2017. Plants were taller in irrigated and higher rainfall environments (17CVI and 17EI), almost three times taller in these optimum environments. Test weight was also affected by drought. The lowest TW was observed in 12EP1 and 12EP3, while the highest TW was in 17BI. Grain yield was lowest in the driest year and dry environments. However, the top 20 high yielding lines across individual environments were distributed within the 106 lines out of 124 lines (Table S3), showing the strong genotype-by-environment interactions and the necessary for mega-environement classification.

Mega-environment analysis categorized all the environments into two or four different mega-environments for different traits (Fig. 1). Grain yield had ME1 (11CH, 14EP4, 14EP5, 17BI, and 17CVI), ME2 (11BD, 11EP5, 12BD, 12EP1, 12EP2, and 17EI), ME3 (12CH, 12UVD, 13UVL, and 17BD), and ME4 (12UV5, 12UV7, 12UVL, and 14CH); Test weight had ME1 (12BD, 12EP2, 12UV7, 12UVL, 12UV5, 13EP3, and 14CH) and ME2 (11EP5, 12CH, 12EP1, 12EP3, 13EP2, 13EP4, 13EP5, 14EP4, 14EP5, 17BD, 17BI, and 17EI). Plant height had ME1 (11EP3, 11BD, 12BD, 17CVI, and 12CH) and ME2 (11EP1, 11EP2, 11EP4, 11EP5, 17BD, and 17EI). Heading date had ME1 (12EP2, 11EP4, 11BD, 12EP1, and 17BD) and ME2 (12EP3, 12BD, 11EP5, 11EP1, and 11EP3). 11EP2 was excluded since it was far deviated from the rest (Fig. 1). QTL was analyzed within and across ME’s for yield and test weight.

Linkage map and QTL analyses

We used JoinMap 4.0 and QTL ICIMapping 4.2 to construct genetic maps and QTL analyses following similar procedures as Yang et al. (2019).

A set of 115 unique QTL regions significantly associated with YLD and related agronomic traits across 28 environments over five years was identified through the analyses of data from IE, MET, ME (Table 3, Table S4). Among them, 51 unique consistent QTL were associated with a single trait from at least two out of IE, MET, ME analyses (Table 3, Fig. 2, Figs. S2, and S3). These consistent QTL were identified for all the four traits analyzed on 15 chromosomes including 1A, 1D, 2B, 2D, 3B, 3D, 4B, 4D, 5A, 5B, 6A, 6B, 6D, 7B, and 7D.

Grain yield

A total of 18 QTL for YLD were detected. Nine common consistent and pleiotropic QTL were located on chromosomes 1D, 4B, 4D, 6A, 7B, and 7D. The other nine consistent QTL were mapped on chromosomes 1A, 3B, 4B, 4D, and 7D (Table 3, Table S4, and Fig. 2). Favorable alleles from TAM 111 for QTL on the chromosomes 1A (411.7 and 585.6 Mb), 3B (48.6 Mb), 4B (226.8 Mb), 4D (445.5 Mb), 6A (12.4 and 608.5 Mb), 7B (15.6 Mb), and 7D (84.3 Mb) explained up to 41.3% of total phenotypic variations and increased yield by up to 37.41 g m⁻² from Qyld.tamu.4B.267 in 17CVI while the remaining QTL with favorable alleles from TAM 112 increased yield by up to 13.54 g m^-2 and explained a total of 25% of yield variations at ME3 including 12CH, 12UVD, 13UVL, and 17BD from Qyld.tamu.1A.587 (Table 3, Table S4, and Fig. S3). Favorable allele switched between two parents for QTL on the chromosome 4B at 266.8 Mb and 4D at 445.5 Mb. TAM 111 favorable allele of Qyld.tamu.4B.267 increased yield by 37.41 g m^-2 in 17CVI and increased yield by 3.38 g m^-2 across ME1, including 11CH, 14EP4, 14EP5, 17BI, and 17CVI while TAM 112 allele only increased yield by 0.42 g m⁻² when analyzed across all 28 environments. However, additive-by-environment contributed very large part of PVE. The additive-by-17CVI interaction explained 63.3% out of the total 63.5% PVE (Table S4). The QTL Qyld.tamu.4D.446 had favorable alleles from TAM 111 that increased yield by 8.59 and 3.98 g m⁻² when analyzed for individual environment in 14CH and across ME4, including 12UV5, 12UV7, 12UVL, and 14CH while it had favorable allele from TAM 112 that increased yield by 0.47 g m^-2 when analyzed across all 28 environments (Table 3, Table S4).

Among the major consistent QTL for yield, those QTL that had larger additive effect contributions from across environment analyses are of interest to breeders for yield improvement. Two QTL, Qyld.tamu.1A.587 and Qyld.tamu.7D.61, increased yield by 4.43 and 6.57 g m^-2 with TAM 112 allele contributing 14.6% and 13.3% additive effect for percentage of variations explained (PVE) across the four mega-environments and ME3 (including 12CH, 12UVD, 13UVL, and 17BD) analyses, respectively (Table S4). A third QTL at 659.2 Mb on 4B increased yield by 6.43 g m^-2 with TAM 112 allele contributing 4.9% additive PVE. A set of four additional QTL, at 20.6 and 109.8 Mb on 4D, 15.6 Mb on 7B, and 84.3 Mb on 7D, increased yield by 2.55 to 8.16 g m^-2 with TAM 111 alleles contributing additive PVE ranging from 4.9% to 7.9%. These seven yield QTL will be our focus for further application.

The additive-by-environment interaction increased yield by up to 37.83 g m^-2 from Qyld.tamu.4B.267 under 17CVI with TAM 111 allele while the same QTL had favorable allele from TAM 112 that only increased 11.13 g m⁻² from interaction under 17BI when it was analyzed for across all 28 environments. Furthermore, the same QTL had interactional effect that increased yield by 34.03 g m^-2 under 17CVI with TAM 111 allele and increased yield by 14.92 g m^-2 under 17BI with TAM 112 allele from analyses of across ME1 including 11CH, 14EP4, 14EP5, 17BI, and 17CVI. Similarly, Qyld.tamu.7D.84 had interactional effect increased yield by 20.3 g m^-2 under 17BI with TAM 111 allele and increased by 11.12 g m^-2 under 12UV5 with TAM 112 allele when analyzed across all 28 environments. The same QTL had an interactional effect of increasing yield by 13.53 g m^-2 with TAM 111 allele under 17BI when it was analyzed across ME1 including 11CH, 14EP4, 14EP5, 17BI, and 17CVI.

Five major QTL had QTL-by-environment interactions from TAM 111 alleles. Qyld.tamu.3B.49 increased yield by 10.8 g m^-2 under 12UVL when it was analyzed across ME4 including 12UV5, 12UV7, 12UVL, and 14CH. The effect of increasing yield by 11.01 and 26.24 g m⁻² were observed when all 28 environments were analyzed for QTL Qyld.tamu.4D.446 under 11EP4, and Qyld.tamu.4D.110 under 17BI, respectively. The two QTL, Qyld.tamu.7D.84 and Qyld.tamu.4D.110 had interactional effect of increasing yield by 13.53 and 19.75 g m⁻² with TAM 111 allele under 17BI from across ME1 analyses (Table S4).

Test weight

A total of 32 unique QTL for test weight were identified. Six consistent and pleiotropic QTL were mapped on chromosomes 1A at 11.8 Mb, 1D at 421.8 Mb, 4D at 109.8 Mb, 6A at 12.4 Mb and 608.5 Mb, and 7D at 60.6 Mb. Additional 26 consistent QTL were located on 14 chromosomes (Table 3, Fig. 2, and Fig. S2). The favorable alleles for TW were from TAM 111 on chromosome 1A at 11.8, 380.7, and 485.2 Mb, 3B at 507.0 Mb, 3D at 548.6, 554.7, and 562.7 Mb, 5B at 589.4 and 646.0 Mb, 6A at 603.3, 608.5 and 611.6 Mb, 6D at 459.2 Mb, 7B at 8.5 and 64.5 Mb, and 7D at 604 Mb that explained up to 11.3% phenotypic variation and increased TW by 7.23 kg m⁻³ with Qtw.tamu.6D.459 in 12UV5 (Table 3, Table S4). The rest of the QTL had favorable allele from TAM 112 that increased TW by up to 10.35 kg m⁻³ with Qtw.tamu.4D.63 in 13EP2. The QTL on the chromosome 6A at 603.3 Mb received favorable alleles from both parents depending on the particular environments. Qtw.tamu.6A.603 explained up to 21.4% of the phenotypic variation and increased TW by 16.58 kg m⁻³ with alleles from TAM 112 in 12UV7 but it only increased 0.08 kg m⁻³ with TAM 111 allele in the analyses across all 19 environments.

Based on PVE of additive effects, only Qtw.tamu.1D.422 along with other four QTL on 4D at 109.8 Mb, 6A at 12.4 Mb, and on 6B at 466.0 and 559.4 Mb, had values ranging from 3.8% to 10.2% with TAM 112 alleles that increased test weight by 2.26 to 3.47 kg m⁻³ (Table S4). One QTL Qtw.tamu.6D.459 increased test weight by 3.66 kg m⁻³ with TAM 111 allele.

For the QTL-by-environmental interactions based on across all individual environment and across ME2 analyses, Qtw.tamu.4D.63 had additive effects of 7.34 and 7.05 kg m⁻³ in environment 13EP2, respectively while their QTL additive effects only increased 2.89 and 3.18 kg/m3 (Table S4). Qtw.tamu.6A.603 had interactional effect of increasing TW by 16.65 and 15.14 kg m⁻³ under the environment 12UV7 from the analyses of across all environments and ME1, respectively (Table S2). All these favorable alleles were from TAM 112.

Pleiotropic QTL to grain yield and test weight

From the multiple trait QTL analyses, nine unique consistent QTL regions for grain yield are also linked to agronomic traits (Table 3, Table S4, and Fig. 2). They were detected on chromosome 1A, 1D, 4B, 4D, 6A, 7B, and 7D. QTL linked to HD and TW on chromosomes 1A at 11.8 Mb increased HD by 0.28 day with TAM 112 allele while the TAM 111 alleles increased TW by 1.87 kg m⁻³ (Table 3, Table S4). QTL on 1D at 421.8 Mb had a favorable allele from TAM 112 that increased TW by 6.41 kg m⁻³ and YLD by 6.87 g m^-2. QTL on chromosome 4B at 659.2 Mb was associated with HD and YLD. Allele from TAM 112 increased HD by 1.53 days and YLD by 9.88 g m^-2. Plant height and YLD were associated with QTL on chromosome 4D at 20.6 Mb. TAM 112 allele increased HT by 0.47 cm under combined DRY while TAM 111 allele increased YLD by 17.32 g m^-2 under 17CVI. Another QTL on chromosome 4D at 109.8 Mb was associated with TW and YLD. TAM 112 allele increased TW by up to 8.01 kg m⁻³ in 13EP3 while TAM 111 allele increased YLD by 28.06 g m⁻² under 17BI. The two QTL on chromosome 6A at 12.4 Mb and 608.5 Mb were associated with both TW and YLD. TAM 112 allele of QTL at 12.4 Mb increased TW by 3.58 kg m⁻³ while TAM 111 allele increased YLD by 19.58 g m^-2. For the QTL at 608.5 Mb, TAM 111 allele increased TW by 5.19 kg m⁻³ under 13EP5 and YLD by 3.18 g m^-2 under ME3. Heading date and YLD were co-located on chromosome 7B at 15.6 Mb. TAM 112 allele increased HD by 0.25 day while TAM 111 allele increased YLD by up to 7.23 g m^-2 under 12CH. QTL on chromosome 7D at 60.6 Mb was associated with HT, TW, and YLD. The favorable allele from TAM 111 increased HT by 1.05 cm under 17EI while the favorable allele from TAM 112 increased TW and YLD by up to 8.48 kg m⁻³ under 11EP5 and 8.73 g m^-2 under 12UVD, respectively. Another QTL on chromosome 7D at 64.3 Mb was associated with HT and yield. TAM 111 allele increased HT by 1.63 cm under 11EP5 and TAM 112 allele increased yield by 6.45 g m^-2 under 12CH. Among the five pleiotropic QTL associated with both TW and yield, the two QTL on 1D and 7D have favorable alleles increasing both traits from TAM 112 while the QTL on 4D at 109.8 and 6A at 12.4 Mb had TAM 112 allele for TW and TAM 111 allele for yield. The QTL on 6A at 608.5 Mb had TAM 111 allele increasing both traits (Table 3, Table S4).

Interactions of epistasis, epistasis-by-environment, and additive-by-environment

There were 359 interactions of additive-by-additive, additive-by-environment, and epistasis-by-environment with a total LOD ≥ 12 for all traits (Table S5, Fig. S4). Only 139 out of 359 interactions had LOD(AA) >10.0. There were only two interactions for HD, one for yield, and the rest of the 136 for TW. None of the interactions was the same as any of those major consistent and pleiotropic QTL for heading date and yield while there were five major consistent QTL associated with the epistasis interactions for TW (Table S5). They were a marker linked to Qtw.tamu.5A.702 interacting with two other markers on 1A at 465.5 Mb and on 5A at 584.4 Mb, a marker linked to Qtw.tamu.5B.646 interacting with IWB5813 on 2B at 25.2 Mb, a marker linked to Qtw.tamu.6A.12 interacting with a marker on 1B at 466 Mb, IWB38972 linked to Qtw.tamu.6B.130 interacting with a marker on 1A at 544.6 Mb, IWB6455 linked to Qtw.tamu.7B.9 interacting with a marker on 1D at 418.5 Mb. However, they only explained TW variations by 1.3% to 2.8% with additive effects from epistasis less than 0.25 kg m⁻³. Furthermore, neither of any interactions from epistasis, the interactions between either marker from the epistasis with the environment, or the epistasis-by-environment interactions had effects that increased TW by more than 1 kg m⁻³.

For the interaction of yield between a marker on 5B at 655.5 Mb and a marker on 7B at 740.1 Mb, its epistasis increased yield by 33.74 g m^-2. The marker on 5B at 655.5 Mb increased yield by more than 30 g m^-2 with TAM 111 allele under drier environments 11CH, 12EP3, and 12UV5, and with TAM 112 allele under higher irrigated environments 11EP5, 12UV7, and 17BI. The marker on 7B at 740.1 Mb increased yield by more than 30 g m^-2 with TAM 111 alleles under environments 11EP5, 12UV5, and 13UVL, and with TAM 112 alleles under environments 12CH and 17EI. The epistasis-by-environment interaction effects increased yield by more than 30 g m⁻² with TAM 111 alleles under less irrigated environments 11EP2, 12UV7, 13EP2, and 13EP3, but with TAM 112 alleles under highly irrigated environments 11EP3, 11EP4, 11EP5, 17BI and 17CVI (Table S5). Neither marker regions were involved with major QTL for yield.

Among the other 220 interactions with LOD(AA) <10.0, six interactions for HD, eight interactions for HT, 53 interactions for TW, and 153 interactions for yield. For TW, there were five interactions that increased TW by more than 0.3 kg m⁻³ with TAM 111 alleles and 14 interactions with TAM 112 alleles but all of them had effect less than 0.5 kg m⁻³. One major QTL Qtw.tamu.4D.63 interacting with 12EP1, 12EP3, 17BD increased TW by 0.38–0.44 kg m⁻³ with TAM 112 allele while the other two major QTL Qtw.tamu.2D.487 and Qtw.tamu.7B.9 had interactional effects less than 0.3 kg m⁻³ (Table S4).

For yield, a total of 1,092 interactions increased yield by more than 50 g m^-2 based on the 153 pairs of QTL across 28 environments with 513 from TAM 112 alleles and 579 from TAM 111 alleles. A subset of 87 interactions increased yield by more than 100 g m^-2. At the first locus of the epistasis, environments 11EP4, 11EP5, 12UV5, and 17BI interacted with a set of 40, 55, 43, and 41 markers increased YLD more than 50 g m^-2 while the corresponding subsets of four, three, three, and one marker increased yield by more than 100 g m^-2. At the second locus, environments 11EP4, 11EP5, 12UV7, 12UVL, 17BI interacted with a set of 35, 40, 29, 17, and 23 markers increased yield by more than 50 g m^-2 while the corresponding subsets of one, five, two, three, and two markers increased yield by more than 100 g m^-2. For epistasis-by-environment interactions, environments 11EP3, 11EP4, 11EP5, 12UV5, 12UV7, 12UVL, 17BI, and 17CVI interacted with a set of 38, 52, 55, 75, 46, 45, 75, and 29 epistasis increased yield by more than 50 g m⁻² while the corresponding subset of four, seven, 14, 14, two, three, 18, and one marker increased yield by more than 100 g m^-2 with three interactions. The top four marker pairs increased yield by more than 100 g m⁻² with additive-by-environment or epistasis-by-environment interactions. A marker on 6A at 613.8 Mb interacted with 12UV5 increased yield by 127.2 g m^-2 with TAM 111 allele. The epistasis between 6A marker with a marker on 7D at 585.6 Mb increased yield by 105.7 g m⁻² under 12UV5 with TAM 112 allele while it increased yield by 106.2 g m^-2 under 11EP4 with TAM 111 alleles. The 2nd set of three interactions were as follow: a marker on 4A at 378.3 Mb and IWA5751 on 4D at 408.8 Mb interacted with 11EP4 increased yield by 104.19 g m⁻² with TAM 112 allele and by 110.72 g m⁻² with TAM 111 allele, respectively; the epistasis between these two markers interacted with 12UV5 increased yield by 102.07 g m^-2 with TAM 111 allele. The 3rd set of three interactions were: marker IWB52359 on 7D at 40.1 Mb interacted with 12UVL and 12UV7 increased yield by 128.98 and 109.82 g m^-2 with TAM 112 allele, respectively; its epistasis with a marker on 7B at 6.8 Mb interacting with 17BI increased yield by 146.04 g m^-2 with TAM 111 allele. The fourth set of three interactions were all epistasis-by-environment effects between IWB73713 on 1B at 675.6 Mb and IWA1924 on 6D at 461.4 Mb that increased yield by 105.73, 106.52, and 106.25 g m^-2, respectively, under 11EP5, 12UV5, and 12UVL with the first interactional allele from TAM 111 and the rest two from TAM 112 (Table S5). The highest effect from TAM 111 allele increased yield by 184.15 g m^-2 that was from a major QTL Qyld.tamu.6A.12 under 17BI while the highest effect from TAM 112 alleles increased yield by 155.18 g m^-2 was from epistasis-by-environment effect between IWA4746 on 2D at 14.4 Mb and a marker on 3A at 7.6 Mb under 11EP5. Two other major QTL, Qyld.tamu.4D.21 and Qyld.tamu.4B.659, were involved with epistasis-by-environment interactions but most interactional effect only increased yield by less than 100 g m^-2 (Table S5).