Genome wide identification of QTL associated with yield and yield components in two popular wheat cultivars TAM 111 and TAM 112

Abstract

Two drought-tolerant wheat cultivars, ‘TAM 111’ and ‘TAM 112’, have been widely grown in the Southern Great Plains of the U.S. and used as parents in many wheat breeding programs worldwide. This study aimed to reveal genetic control of yield and yield components in the two cultivars under both dryland and irrigated conditions. A mapping population containing 124 F5:7 recombinant inbred lines (RILs) was developed from the cross of TAM 112/TAM 111. A set of 5,948 SNPs from the wheat 90K iSelect array and double digest restriction-site associated DNA sequencing was used to construct high-density genetic maps. Data for yield and yield components were obtained from 11 environments. QTL analyses were performed based on 11 individual environments, across all environments, within and across mega-environments. Thirty-six unique consistent QTL regions were distributed on 13 chromosomes including 1A, 1B, 1D, 2A, 2D, 3D, 4B, 4D, 6A, 6B, 6D, 7B, and 7D. Ten unique QTL with pleiotropic effects were identified on four chromosomes and eight were in common with the consistent QTL. These QTL increased dry biomass grain yield by 16.3 g m-2, plot yield by 28.1 g m-2, kernels spike-1 by 0.7, spikes m-2 by 14.8, thousand kernel weight by 0.9 g with favorable alleles from either parent. TAM 112 alleles mainly increased spikes m-2 and thousand kernel weight while TMA 111 alleles increased kernels spike-1, harvest index and grain yield. The saturated genetic map and markers linked to significant QTL from this study will be very useful in developing high throughput genotyping markers for tracking the desirable haplotypes of these important yield-related traits in popular parental cultivars.

About the Authors:

Yan Yang

Roles Data curation, Formal analysis, Funding acquisition, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

Affiliations Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America, Department of Soil and Crop Science, Texas A&M University, College Station, Texas, United States of America

Smit Dhakal

Roles Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Validation, Visualization, Writing – review & editing

Affiliations Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America, Department of Soil and Crop Science, Texas A&M University, College Station, Texas, United States of America

ORCID logo https://orcid.org/0000-0002-1166-1168

Chenggen Chu

Roles Data curation, Formal analysis, Funding acquisition, Methodology, Software, Writing – review & editing

Current address: USDA-ARS, Sugarbeet & Potato Research Unit, Edward T. Schafer Agricultural Research Center, Fargo, ND, United States of America

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

ORCID logo https://orcid.org/0000-0001-6189-5985

Shichen Wang

Roles Data curation, Formal analysis, Funding acquisition, Methodology, Software, Writing – review & editing

Affiliation: Genomics and Bioinformatics Service Center, Texas A&M AgriLife Research, College Station, TX, United States of America

Qingwu Xue

Roles Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Project administration, Resources, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Jackie C. Rudd

Roles Conceptualization, Data curation, Funding acquisition, Investigation, Project administration, Resources, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Amir M. H. Ibrahim

Roles Conceptualization, Data curation, Funding acquisition, Resources, Writing – review & editing

Affiliation: Department of Soil and Crop Science, Texas A&M University, College Station, Texas, United States of America

Kirk Jessup

Roles Data curation, Investigation, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Jason Baker

Roles Data curation, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Maria Pilar Fuentealba

Roles Data curation, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Ravindra Devkota

Roles Data curation, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Shannon Baker

Roles Data curation, Writing – review & editing

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Charles D. Johnson

Roles Data curation, Investigation, Resources, Writing – review & editing

Affiliation: Genomics and Bioinformatics Service Center, Texas A&M AgriLife Research, College Station, TX, United States of America

Richard Metz

Roles Data curation, Writing – review & editing

Affiliation: Genomics and Bioinformatics Service Center, Texas A&M AgriLife Research, College Station, TX, United States of America

Shuyu Liu

Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Texas A&M AgriLife Research and Extension Center, Amarillo, Texas, United States of America

Introduction

Wheat (Triticum aestivum L.) is one of the most important food crops worldwide. The significance of wheat lies on its physical and chemical properties of grains, which provide over 20% of the calories and protein requirements for human nutrition. Yield is a polygenic complex trait and the most important to breeders and farmers. However, environmental conditions and the genetic-by-environmental interactions throughout all processes of vegetative and reproductive growth and development could seriously affect yield [1]. In general, grain yield can be broken into three major components as number of spikes m-2 (SPM), kernels spike-1 (KPS), and thousand kernel weight (TKW) with each controlled by multiple genes or quantitative trait loci (QTL). Interactions among QTL and between QTL and environments also modify the expression of the QTL in different genetic backgrounds [2]. Typically, a QTL detected in one environment but not in another might be a indication of QTL × environment interaction (QEI). However, assessing the effects of such interactions is difficult due to the unpredictable random change of environments. Goldringer et al. [3] proposed the additive and epistatic genetic variances for agronomic traits in a doubled haploid population and demonstrated that yield and its components showed either additive or additive plus epistatic effects. Significant epistasis and QEI for yield were identified subsequently in other researches [4–8]. Thus, dissection of QTL effects and their interactions may facilitate better understanding of the genetic control of the complex yield traits [9].

Saturated genetic linkage maps play a crucial role in QTL identification for providing measurements of the relative effects of alleles in a mapped chromosomal region as well as selectable DNA markers for breeders to integrate the traits through marker-assisted selection (MAS) [10]. More recently, single nucleotide polymorphisms (SNPs) as the common source of genetic variation among individuals of any species and the smallest unit of genetic variation with virtually unlimited numbers, were used to develop high-density linkage maps and QTL identification in many crops. The availability of diverse SNP genotyping platforms, particularly genotyping-by-sequencing (GBS) based on the next-generation sequencing, were facilitated in genetic dissection, marker discovery, and genomic selection of complex traits [11, 12]. However, the extensive abundance of conserved repetitive element nature of the hexaploid wheat genome (~80%) has slowed the progress in SNP discovery and detection [13]. Cavanagh et al. [14] developed 9K SNP assays and constructed the first high-density wheat consensus SNP map containing 7,504 polymorphic loci. A set of 40K out of 90K SNP assay from wheat was mapped onto chromosomes [15], thus provides a powerful resource for genome-wide dissecting traits of interests and developing new tools for efficient selection in breeding. Liu et al. [16] mapped 4k to 8k array SNPs in three wheat bi-parental mapping populations.

In this study, the highly-saturated genetic maps constructed with SNPs from 90K iSelect array and double digest restriction-site associated DNA sequencing (ddRADseq) were used to dissect QTL associated with yield, yield components, and related traits in popular cultivars TAM 111 and TAM 112. Additionally, through extensive analysis of additive-by-environment interactions, epistasis, and epistasis-by-environment interactions in individual and mega environments, the consistent and pleiotropic QTL were identified and summarized.

Materials and methods

Plant material and phenotyping

A population of 124 F5:7 recombinant inbred lines (RILs) was derived from the cross between TAM 112 and TAM 111. Both the parents are hard red winter wheat (HRWW) released by Texas A&M AgriLife Research, and they are the top-ranked cultivars grown in the U.S. Great Plains. TAM 111 has the pedigree of ‘TAM 107’//TX78V3630/ ‘Centurk78’/3/TX87V1233 with excellent performance under both drought and irrigated conditions, whereas TAM 112 has the pedigree of U1254–7–9–2-1/TXGH10440 and is highly adapted to drought condition [17, 18]. Genetic analysis of the population thus can detect the favorable alleles from the two parents.

The 124 recombinant inbred line (RILs) of TAM 112/TAM 111 along with their parents were evaluated for yield and yield component traits in field experiments across 11 environments during five crop years harvested in 2011, 2012, 2013, 2014, and 2017. The combination of the location-year-irrigation level is an environment. Field locations used in this study included Texas AgriLife Research stations in Bushland (35° 06' N, 102° 27' W) in 2011, 2012 and 2017 (designated as 11BD, 12BD for dryland and 17BI for irrigated, respectively), Chillicothe (34° 07' N, 99° 18' W) in 2012 and 2014 (designated as 12CH and 14CH, respectively), two irrigation levels (75% and 100%) in Etter (35° 59' N, 101° 59' W), TX in 2013 and 2014 (designated as 13EP4, 13EP5, 14EP4 and 14EP5, accordingly), irrigated plots at Clovis (34° 24' N, 103° 12' W), NM (designated as 17CVI), and Dumas (35° 51' N, 101° 58' W) (designated as 17DMS), TX in 2017. All of these locations are belonging to the Texas A&M AgriLife Research or New Mexico State University experimental stations so no specific permission was required. It is confirmed that the field studies did not involve endangered or protected species. All trials were planted using alpha lattice design with an incomplete block size of five plots, and each trial has two replications in every environment. Standard agronomic practices were carried out for each environment. The data collection followed similar procedures as outlines by Assanga et al. [19]. Plot grain yield from combine harvester (YLD) was recorded. Biomass sample harvested from a random half-meter inner row with representing uniform plant performance of each plot was oven-dried for 72 hrs at 60°C and used to measure total dry biomass (BM), grain weight from the biomass sample (BMYLD), and yield components. Thousand kernel weight (TKW) was calculated by the weight of 200 seeds and scaling to 1000 seeds from biomass sample; the number of spikes m-2 (SPM) was calculated from the plot sample by counting the number of heads. Kernels spike-1 (KPS) was calculated using BMYLD, TKW and SPM; the harvest index (HI) was calculated as grain weight (BMYLD) divided by total weight of biomass sample (BM) from each plot. Single head dry weight (SHDW) was calculated through dividing the total dry head weight including glumes and awns per plot biomass sample by the number of heads. Single head grain weight (SHGW) was calculated by dividing the total BMYLD by the number of heads.

DNA extraction and genotyping

Whole genomic DNA was extracted from leaf samples of parents and 124 RILs using the CTAB method with minor modification as described by Liu et al. [20]. SNP genotyping with Infinium iSelect assays containing 90K SNPs was performed in USDA Small Grains Genotyping Laboratory at Fargo, ND according to manufacturer’s protocol (Illumina Inc., San Diego, CA, USA), and the assay was designed under the International Wheat SNP Consortium protocols [14]. The fluorescence signal was captured by Illumina scanner and analyzed using GenomeStudio software (www.illumina.com). More details for polymorphic SNP sorting and conversion in this population were outlined in Liu et al. [16] and Dhakal et al. [21].

The 124 RILs and two parents were also genotyped in Texas AgriLife Research Genomics & Bioinformatics Services at College Station, TX (http://www.txgen.tamu.edu/) following the Double digest restriction-site associated DNA (ddRADSeq) method with some noted modifications. The libraries were constructed using a 96-plex plate with single random blank well included for quality control. DNA was co-digested with the restriction enzymes PstI (CTGCAG) and MspI (CCGG), and barcoded adapters were ligated to individual samples. SNP calling was processed as described by Yang et al. [22]. More details are described below with some modification from Elshire et al. [23], Poland et al. [24], and Ballare et al. [25].

Adapters: All oligos were purchased from Integrated DNA Technologies (IDT), and were received as a 100 μM stock. Adapters were made by mixing equimolar amounts (30 μM of top and bottom oligos in 100 μl 1x annealing buffer (10 mM Tris-HCl, pH 8.0, 50 mM NaCl, 1 mM EDTA). The oligos were held at 95°C for 10 sec, then cooled to 12°C at a rate of 0.1°C per sec. P5-Index Adapters were made by annealing the following oligos (where XXXXXXXX represents 8-base i5 index sequences): Top (5’ to 3’): AAT GAT ACG GCG ACC ACC GAG ATC TAC ACX XXX XXX XTC TTT CCC T, Bottom (5’ to 3’): /5Phos/AXX XXX XXX GTG TAG ATC TCG GTG GTC GCC GTA TCA TT. The P5-PstI-Bridge adapters was made by annealing top (Pster_T, 5’ to 3’): /5Phos/ACA CGA CGC TCT TCC GAT CTT GCA and bottom (Pster_B, 5’ to 3’): AGA TCG GAA GAG CGT CGT GTA GGG AAA G oligos. P7-MluCI Adapter was made by annealing top (P7-MluCI_T, 5’ to 3’): AAT TAG ATC GGA AGA GCA CAC GTC TGA ACT CCA GTC AC and bottom (P7-MluCI_B, 5’ to 3’): GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATC T.

Dual ligation ddRAD: At the end of each step in this protocol, samples were quantified on a DeNovix DS-11 spectrophotometer. Sample DNA concentrations were determined using Quant-iT-Picogreen reagent (Thermo Fisher Scientific) and read on a SpectraMax Spectrofluorometric plate reader (Molecular Devices, LLC). One hundred nanograms of DNA per sample in 96 well plate format was digested in a final volume of 25 μl in 1x NEB Cut Smart Buffer and 200 U PstI-HF and 100 U MluCI (New England Biolabs) at 37°C for 4 hours. Following a 20 min 80°C enzyme inactivation, samples were held at 12°C until ligation. To each 25 μl digest, 5.7 μl of a master mix was added such that each well got the equivalent of 3.2 μl 10x Ligase buffer (NEB), 0.25 μl T4 DNA Ligase (New England Biolabs) and P5-PstI-Bridge and P7-MluCI adapters at a final concentration of 500 nM. In addition, each well got 1 of 48 unique P5-Index Adapters (400 nM final concentration) and were mixed by pipetting. Plates were spun down and incubated at 16°C for 8 hr followed by a 15 min heat inactivation at 65°C and held at 4°C. Pools were made by combining no more than 48 samples with unique P5 indexes. To each pool, EDTA was added to 0.25 mM to further inhibit ligase activity. Pools were precipitated by adding 1/10th volume of 3 M sodium acetate (pH 5.2), evenly dividing them into two or three 2.0 ml microcentrifuge tubes and adding 2 volumes 100% ethanol and placing them at -20°C for at least 1 hr. Tubes were spun at 20,000 xg for 10 min and supernatant poured out. Tubes were washed with 1 ml 70% ethanol, spun 5 min and supernatant removed. Pellets were resuspended in 200 μl EB and purified through QIAquick PCR Purification Kits as per manufacturer’s protocol (Qiagen) eluting twice with 50 μl EB. Combined elutants (100 μl total) were further purified to remove unligated adapters using 0.9x volume AMPure XP beads as per the manufacturer’s protocol (Beckman-Coulter) eluting in 35 μl EB. Up to 3000 ng of each pool was size selected at 390–610 bp (280–600 bp inserts plus 110 bp adapters) on Pippin prep 2% dye-free gels (Sage Science). Recovered DNA was purified with 0.9x AMPure XP beads as described earlier (Beckman-Coulter) and eluted in 32 μl EB.

Incorporation of a biotin moiety at the P5 side (for further purification—described later) and the addition of the i7 index was accomplished in a Pre-Selection PCR step. Using the primers P5_Select (5’-3’): /5BiotinTEG/AAT GAT ACG GCG ACC ACC GAG ATC TAC AC and one of twenty-four i7 indexed reverse primers (TDX 1–24, 5’-3’ where XXXXXXX represent bases used for indexes: CAAGCA GAA GAC GGC ATA CGA GAT XXX XXX XGTGAC TGG AGT TCA GAC GTG TGC). PreSelect PCR reactions (200 μl total volume, split into two reactions of 100 μl each) contained up to 150 ng size-selected DNA, 0.4 mM dNTPs, 1 μM each primer (P5-Select and a TDX-reverse index), 20 U of Q5 Polymerase in 1x Q5 High Fidelty DNA Polymerase Buffer. Reactions were denatured at 98°C for 30 sec, then subjected to 15 cycles of 98°C for 10 sec, 62°C for 20 sec and 72°C for 40 sec with a final elongation at 72°C for 5 min followed by a 10°C hold. Pre-Selection PCR reactions were cleaned up with QIA quick columns and AMPure XP beads as described above with a final elution in 50μl EB.

Selection of only PstI-MluCI fragments was accomplished using Dynabeads M-280 Streptavidin beads (Invitrogen) to capture fragments with biotin incorporated at their P5 ends during PreSelect PCR. Dynabeads (50 μl per pool) were captured on a magnet and washed twice in 300 μl 1x bead washing buffer (10 mM Tris, pH 7.5, 2 M NaCl, 1 mM EDTA) by resuspending beads in buffer, capturing on a magnet, and removing the supernatant with pipette tips. After the second wash, beads were resuspended in 100 μl 2x bead wash buffer per pool, and 100 μl washed beads was mixed with up to 3000 ng of Pre-Selected DNA (in 100 μl EB). Samples were incubated at RT for 20 min then captured on a magnet. MluCI-MluCI fragments lacking biotin were washed away as follows: beads were washed three times in 200 μl 1x bead washing buffer, twice in 200 μl nuclease free water and once in 100 μl 1x SSC (150 mM NaCl, 15 mM sodium citrate, pH 7.0, BIO-RAD Laboratories, Inc.) and finally resuspended in 50 μl 1x SSC. PstI-MluCI fragments were obtained by heating beads at 95°C for 5 min thus denaturing off the non-biotinylated strand (leaving both strands of PstI-PstI fragments attached to the beads). Following 5 min at 95°C, tubes were transferred immediately to a magnet and supernatant removed quickly to new tubes. This was repeated for a total of two heated elutions with a final volume of 100 μl. Elutants were purified with QIA quick PCR columns as described above, eluted in 40 μl EB and quantified.

Final libraries were produced in a PCR reaction of 50 μl containing 10 ng single-stranded, Dynabead-selected DNA, 0.4 μM dNTPS, 0.5 μM each final PCR primers (DuLig-F1, 5'-AAT GAT ACG GCG ACC ACC GAG ATC TAC AC-3' and DuLig-R1, 5'-CAA GCA GAA GAC GGC ATA CGA GAT-3') and 20 U/μl Q5 DNA Polymerase in 1x Q5 DNA polymerase reaction buffer (New England Biolabs). Reaction conditions were the same as Pre-Selection PCR, but with total eight cycles. Final PCR reactions were cleaned up with 0.9x AMPure XP beads, eluted in 35 μl EB, quantified and assessed for quality on a Fragment Analyzer (Agilent Technologies) diluted and quantified by qPCR (Kappa Biosystems). Libraries were sequenced on the Illumina NovaSeq 6000 system, S4XP flowcell running 2x 150 bp recipe.

Statistical analysis

The analysis of variance (ANOVA) from individual and across environment data was calculated to determine the significance of genetic (G), environment (E), and genetic-by-environment interaction (GEI) variances. Broad-sense heritability was calculated, and only single environment with heritability ≥ 0.05 were included into the analysis. Pearson's correlation coefficients among all variables were calculated. Best linear unbiased predictors (BLUP) and best linear unbiased estimator (BLUE) of individual environment and across all environments were computed using a restricted maximum likelihood (REML) approach based on META-R program with lme4 package in R software from Alvarado rt al. [26]. Mega-environments for each trait were classified according to the biplot clustering for the environments. In most case, the BLUP values were used for the QTL analyses. The BLUE values were used only if the BLUP values were the same for all the RILs in that environment.

Linkage map construction and QTL identification

Linkage map construction in this population has been described in Yang et al. [22]. Of the marker data generated, the false double crossovers were manually checked and removed according to the alignment of SNP orders between genetic maps and physical base pair location from IWGSC RefSeq v1.0 (https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies, accessed on February 8, 2020). QTL analysis was performed using QTL IciMapping software [27]. Individual environment QTL analysis was conducted using single trait from single individual environment and across multiple environments. The multi-environment QTL analyses were also performed using single trait across classified mega-environments and within each mega-environment. The genetic position and effects of individual environment QTL and multi-environment trait QTL were determined by integrated composite interval mapping (ICIM) function for additive effect (ICIM-ADD) and epistasis effect (ICIM-EPI). To identify an appropriate threshold likelihood of odd (LOD) score for declaring a significant QTL, permutation test was conducted for 1,000 times for ICIM-ADD for individual and across environments. Consistent QTL was determined if a QTL was significant at least from two individual environments or two out of the four analyses including individual environment, across all individual, within and across mega environments. Pleiotropic QTL was determined if a QTL was significantly associated with two or more traits that were not highly correlated to each other. For ICIM-EPI, since it is too long to run permutation, LOD = 5 was chosen but actual thresholds for each trait from ICIM-ADD as reference or 10 if too many interactions exist.

Identified QTL were designated in the format as Qtrait.tamu.chrom.Mb, where trait is an abbreviation for a trait name, tamu represents Texas A&M University, chrom is chromosome on which the QTL is located, and Mb is a physical position of the QTL peak linked marker according to alignment with the IWGSC RefSeq v1.0 reference genome.

Results

Trait analysis

The combined ANOVA across environments indicated a significant genetic variance for all traits except dry biomass and significant genetic by environment interactions (GEI) except spikes m-2 (P < 0.01) (S1 Table).

The entry-mean heritability ranged from moderate (0.4 to 0.6) to high (> 0.6). Yield and three yield components, including thousand kernel weight, kernels spike-1, spikes m-2, as well as single head dry weight and single head grain weight, exhibited higher heritability (0.75–0.90). Harvest index displayed moderate heritability (0.56), whereas dry biomass and biomass grain yield expressed relatively low heritability of 0.23–0.33. For the overall BLUP means of three duplicated sets across seven to 11 environments, TAM 112 had higher dry biomass and spikes m-2 while TAM 111 had higher kernels spike-1, single head dry weight and grain weight (S1 Table).

Based on the best linear unbiased prediction (BLUP) means, positive genetic correlations were found between yield and the three yield components expect thousand kernel weight (TKW) (S2 Table). Dry biomass had significant positive correlations with all tested traits except thousand kernel weight and spikes m-2. Harvest index had significant correlations with all traits except spike m-2. Thousand kernel weight only had a significant low positive correlation with harvest index but a significant negative correlation with spikes m-2. Spikes m-2 had a significant negative correlation with kernels spike-1. Thousand kernel weight had the least significant correlations with other traits related to yield indicated that it could be improved independently. However, kernels spike-1 can be increased along with thousand kernel weight for improved yield but not together with spikes m-2 due to the significant negative correlations. Across all individual 11 environments, harvest index were significantly and positively correlated with all the traits except thousand kernel weight in 12BD, 17CVI, and 17BI (significantly negative), spikes m-2 in 12BD, 12CH, 17DMS, 11BD and 13EP5 (last two significantly negative), kernels spike-1 in 14CH, and yield in 13EP5 (S2 Table); thousand kernel weight was significantly and negatively correlated with most of the rest traits except kernels spike-1 in 12BD, 12CH, 13EP5, 17CVI, 17DMS, 13EP4, 14CH, and 14EP5 (last three significantly positive), spikes m-2 in 14CH and 14EP5, yield in 13EP4, 14EP4, 17CVI, 12BD, 12CH, 14CH 14EP5, and 17DMS (last five significantly positive); spikes m-2 was significantly and negatively correlated with kernels spike-1 in 11BD,12BD, 12CH, and 13EP4, except in 14CH, 17DMS, and 13EP5 (last one significantly positive), significantly and positively correlated with yield in 11BD, 12BD, and 14CH, except in 12CH, 13EP4, 13EP5, and 14EP5 (last two significantly negative); kernels spike-1 was significantly and positively correlated with yield in all 11 environments (S2 Table).

The significant correlations between yield and its components implied that yield can be increased through the indirect selection of the higher component traits. Hence, mapping the QTL for yield and associated yield components could reveal significant QTL across environments and improve the indirect selections.

Boxplot and biplot across all environments, and mega environment classification

From the boxplot of all the traits across individual environments, it is easy to define that the lower yield environments were from the drought years (eight environments from 11BD to 14EP5) while the higher yield environments were from the irrigated location in a good rainfall year (three environments, 17BI, 17CVI, and 17DMS) with the latter had higher genetic variations (S1 Table). Similar trends were found on dry biomass and biomass grain yield; however, several other traits did not follow this trend. The three environments having higher harvest index were 11BD, 14EP4, and 14EP5, and all the environments had relatively larger variations. Kernels spike-1 had very similar means across all the environments except 17BI and 17DMS. Thousand kernel weight were classified into two groups and the higher median group included all the irrigated environments in 2014 and 2017. However, spikes m-2 did not have similar trends as any other traits. Its median and ranges were very similar in the two irrigated environments, 17BI and 17DMS (S1A Fig).

Biplot of all the environments for each trait could help us to classify those environments where the performance of individual lines had similar trends; therefore, we classified them as a mega environment (ME) (S1B Fig). Yield had ME1 (17BI, 17CVI, 14EP4, 14EP5); ME2 (12BD, 13EP4, 13EP5), and ME3 (12CH, 14CH); dry biomass had ME1 (11BD, 14CH, 14EP4) and ME2 (13EP4, 13EP5); biomass grain yield had ME1 (11BD, 14CH), ME2 (13EP4, 12BD, 12CH), ME3 (14EP4, 14EP5, 17CVI); harvest index had ME1 (14EP5, 17BI, 17CVI, 17DMS), ME2 (12CH, 14CH, 13EP4), ME3 (12BD, 11BD, 13EP5); kernels spike-1 had ME1 (17CVI, 17DMS, 14EP4, 14EP5), ME2 (11BD, 12BD, 13EP4, 14CH); spikes m-2 had ME1 (14EP5, 14CH), ME2 (12BD, 13EP4, 13EP5, 17CVI), ME3 (12CH, 17DMS); thousand kernel weight had ME1 (11BD, 12BD, 12CH, 13EP4, 13EP5, 14CH), ME2 (14EP4, 14EP5, 17BI, 17CVI, 17DMS). If an environment had a small vector length, RILs did not have larger genetic variations so it was omitted in ME classification. The mega environments allowed us to identify some consistent genetic factors across similar individual environments, within and across mega environments.

Linkage map

A set of 5,948 markers including 3,193 from ddRADseq and 2,740 from 90K iSelect array SNPs, and 15 microsatellites and kompetitive allele specific PCR (KASP) markers were used to construct 25 linkage groups covering all 21 chromosomes that were used for QTL analyses (S3 Table). The cumulative genetic map length is 2,703.9 cM with an average marker density of 0.6 SNP/cM or 2.8 SNP/Mb. The total covered physical base pair length is about 12.6 Gb with average length of 602.2 Mb per chromosome.

Consistent QTL identification for individual trait

A set of 87 unique QTL regions significantly associated with nine yield and yield related traits across 11 environments over five years were identified through the analyses of data from individual and mega-environments (S4 Table; S2 and S3 Figs). Among them, a set of 36 unique consistent QTL was identified to be associated with one trait but from at least two out of the analyses from individual, across all individual, within and across each defined mega environments based on biplot and overall best linear unbiased prediction (BLUP) or best linear unbiased estimation (BLUE) for each trait (Tables 1 and 2; Fig 1). A set of 10 unique pleiotropic QTL were found to be associated with at least two traits that were not highly correlated to each other (Table 2). Among the consistent and pleiotropic QTL, eight were in common (Tables 1 and 2; Fig 1).

[Figure omitted. See PDF.]

Fig 1. Consistent and pleiotropic QTL identified from individual and mega-environments for all traits.

Traits include 1) Yield from combine plots (YLD), 2) dry biomass from hand harvested 0.5 m long inner row sample from crown (BM), 3) grain weight from b) as hand harvested dry grain (BMYLD), 4) harvest index (HI), 5) kernels spike-1 (KPS), 6) spikes m-2 (SPM), 7) thousand kernel weight (TKW), 8) single head dry weight (SHDW), 9) single head grain weight (SHGW). Identified QTL were designated in the format as Qtrait.tamu.chrom.Mb.

https://doi.org/10.1371/journal.pone.0237293.g001

[Figure omitted. See PDF.]

Table 1. Consistent QTL associated with yield and yield components from at least two analyzed across individual or mega environments in TAM 112/TAM 111.

https://doi.org/10.1371/journal.pone.0237293.t001

[Figure omitted. See PDF.]

Table 2. Pleiotropic QTL associated with at least two different traits of yield and yield components in TAM 112/TAM 111.

https://doi.org/10.1371/journal.pone.0237293.t002

Yield

A set of 14 consistent QTL for yield was identified on chromosomes 1B, 1D, 2A, 4B, 4D, 6A, 6D, 7B and 7D (Table 1, S4 Table). There were four major QTL at 20.6 and 109.8 Mb on 4D, 12.4 Mb on 6A, and 90.2 Mb on 7D that increased yield up to 19.6–28.1 g m-2 from the analyses of individual environment 17BI or 17CVI, ME1 (including 17BI, 17CVI, 14EP4, and 14EP5) and across 11 environments and all had favorable alleles from TAM 111. Eight minor QTL with favorable alleles from TAM 112 that increased yield by 2.5–9.9 g m-2 were located at 376.1 Mb on 1B, 421.8 Mb on 1D, 659.2 Mb on 4B, 455.3 Mb on 4D, 19.6 Mb on 6D, 617.0 Mb on 7B, and 64.3 Mb and 591.2 Mb on 7D. From the LOD score and values of percentage of explained variations (PVE) of additive effects, only five out of the eight minor QTL, Qyld.tamu.1B.376, Qyld.tamu.1D.422, Qyld.tamu.7B.16, Qyld.tamu.7B.617, and Qyld.tamu.7D.64 had larger proportion of additive effects while the rest had larger additive-by-environment interactions than additive effects indicating the complex of yield inheritance (S3A Fig). Among the four major QTL that had larger additive effects, results from across individual environment analyses showed that the corresponding additive-by-environment interactions increased yield by 15.8–24.0 g m-2 at 17BI or 17CVI (S4 Table).

TKW

Eight QTL were identified for thousand kernel weight including one on chromosome 1D at 12.3 Mb, four on 2D at 15.7, 63.3, 486.8 and 531.4 Mb, one on 4D at 343.2 Mb, and two on 7D at 40.1 and 64.3 Mb (Table 1). Two QTL Qtkw.tamu.2D.16 and Qtkw.tamu.4D.343 had the favorable alleles from TAM 111 and increased TKW up to 0.3 g while the other six QTL had alleles from TAM 112 and increased TKW up to 0.9 g. All QTL appeared across 11 environments and ME2 (including 14EP4, 14EP5, 17BI, 17CVI, and17DMS) analyses except three QTL Qtkw.tamu.1D.12, Qtkw.tamu.2D.487, and Qtkw.tamu.7D.40. Four major QTL Qtkw.tamu.2D.487, Qtkw.tamu.2D.531, Qtkw.tamu.7D.40 and Qtkw.tamu.7D.64 increased thousand kernel weight from 0.6 to 0.9 g at 12CH, 17CVI, 11BD, and 17DMS, respectively (S4 Table). Their corresponding additive-by-environment interactions increased thousand kernel weight by 0.6, 0.3, 0.5, and 0.3 g, respectively (S4 Table and S3G Fig).

KPS

Only four QTL significantly associated with kernels spike-1 were identified on chromosomes 1A at 13.8 Mb and 411.7 Mb, 4D at 445.5 Mb, and 7B at 647.8 Mb (Table 1). All favorable alleles were from TAM 111 and increased kernels spike-1 up to 0.7 except the Qkps.tamu.4D.446. Two of the three QTL appeared in the analyses of across eight environments, or either 17CVI or 17DMS and ME1 (including 17CVI, 17DMS, 14EP4, and 14EP5). The corresponding additive-by-environment interactions of the three QTL increased kernels spike-1 by 0.3–0.4 at 17DMS or 17CVI (S4 Table and S3E Fig).

SPM

For spikes m-2, three QTL were detected on chromosomes 1A at 356.4 Mb, 4D at 484.7 Mb, and 6B at 673.8 Mb (Table 1). Qspm.tamu.1A.356 and Qspm.tamu.6B.674 had alleles from TAM 112 and increased spikes m-2 by 14.8 while Qspm.tamu.4D.459 had allele from TAM 111 and increased spikes m-2 by 11.5 at 17CVI. All three QTL appeared in the analyses of ME2 (including 12BD, 13EP4, 13EP5, 17CVI) and the two QTL on 1A and 6B appeared in the analyses of across seven environments. Qspm.tamu.6B.674 had the highest additive effects of 14.8 from 12BD and its additive-by-environment interactions increased spikes m-2 by 10.6 while Qspm.tamu.1A.356 increased spikes m-2 by 12.8 and its interactions at 13EP4 increased 8.3 (S4 Table and S3F Fig).

HI

Five QTL for harvest index were detected, in which the two QTL at 1.7 Mb on chromosome 2D and at 29 Mb on 4D had favorable alleles from TAM 111. Qhi.tamu.4D.29 increased harvest index by 1.1% and were consistent in two individual environments, 17CVI and 17DMS. The other three QTL on chromosomes 1B, 3D and 7D had favorable alleles from TAM 112 and increased harvest index by 0.6% at environments 17DMS or 12CH. These five QTL could increase harvest index by 0.25 to 0.87 from additive-by-environment interactions at their corresponding environments, 17DMS, 17CVI or 12CH (Table 1 and S4 Table, S3D Fig).

BMYLD

For biomass grain yield (BMYLD) collected from 0.5-m long in an inner row, three significant QTL were identified on chromosomes 4D, 6D, and 7D (Table 1). Qbmyld.tamu.4D.26 had favorable allele from TAM 111 and increased biomass yield by 16.3 g m-2 at two individual environments 17CVI and 14EP5. The other two QTL at 459.2 Mb on 6D and 64.3 Mb on 7D had favorable alleles from TAM 112 and increased biomass yield up to 13.3 g m-2 at 12CH. Only the major QTL Qbmyld.tamu.4D.26 had a larger additive LOD scores compared with those of corresponding additive-by-environment interactions that increased biomass yield by 7.6 and 11.9 g m-2 at environments 17BI and 17CVI, respectively (S4 Table and S3B Fig). On the other hand, the QTL Qbmyld.tamu.7D.64 had additive-by-environment interaction effect of 16.1 g m-2 with allele from TAM 111 at 17BI from the analyses of across eight environments (S4 Table). At another environment 12CH, the additive-by-environment interaction of the same QTL increased 7.7 g m-2 with allele from TAM 112. Only one significant QTL for dry biomass at 455.8 Mb on chromosome 4D had favorable allele from TAM 112 and increased biomass by 12.1 g m-2 (Table 1).

In general, we can see that TAM 111 favorable alleles mainly increased kernels spike-1 while TAM 112 favorable alleles mainly increased spikes m-2 and thousand kernel weight. For biomass yield, yield, and harvest index, almost half of the QTL had TAM 111 favorable alleles and half had TAM 112 alleles (Table 1). A major QTL had the highest additive effects for certain trait at a particular environment. In the meantime, it had a higher effect from additive-by-environment interactions at the same environment. The additive by environment effects for yield, biomass, biomass grain yield, and harvest index had higher values on explained phenotypic variations and LOD than those of their corresponding additive effects for most QTL while KPS and SPM had higher values from additive components with the TKW had the higher values on additive for all the QTL (S3A–S3G Fig).

Pleiotropic QTL

A set of ten unique significant QTL regions was found to affect more than one trait and thus considered having pleiotropic effects (Table 2 and Fig 1). Eight were in common with the 36 consistent QTL identified for all nine evaluated traits. They were the QTL at 411.7 Mb on 1A that was linked to yield and kernels spike-1 with all the favorable alleles increasing the traits from TAM 111; the QTL at 20.6 Mb on 2D that was linked to both yield and biomass yield and the QTL at 26 Mb on 4D that was linked to yield and harvest index with all favorable alleles from TAM 111; two additional QTL at 455.3 and 455.8 Mb on 4D that were linked to biomass grain yield and yield, dry biomass and single head dry weight, respectively, with all favorable alleles from TAM 112; the fifth QTL on 4D at 445.5 Mb that increased kernels spike-1 by 0.72 with favorable allele from TAM 112 while it increased dry biomass by 54.1 g m-2 based on the BLUE value from 14EP4 (Table 2 and S4 Table); the QTL at 19.6 Mb on 6D that was associated with yield and single head grain weight with favorable alleles from TAM 112; the QTL at 64.3, 66.3, and 68.3 Mb on 7D that were linked to yield, biomass yield, thousand kernel weight, kernels spike-1, and harvest index with the most favorable alleles increased yield and kernels spike-1 from TAM 111 while the favorable alleles increased thousand kernel weight and harvest index from TAM 112 (Table 2). The last two were not consistent QTL (Tables 1 and 2).

Major QTL allele effects comparisons

There were six QTL conditioning spikes m-2 and kernels spike-1 on chromosomes 1A, 4D, 6B and 7B (Table 1). Three QTL, Qspm.tamu.1A.356, Qkps.tamu.4D.446, and Qspm.tamu.6B.674 had favorable alleles from TAM 112 while QTL Qkps.tamu.1A.412, Qspm.tamu.4D.485, and Qkps.tamu.7B.648 had favorable alleles from TAM 111. The combinations of these 6 QTL were classified into 21 groups with two to ten RILs in each group based on 105 RILs with no missing data for those tightly linked SNPs (S4 Fig). Since QTL for SPM on 1A was only presented in 13EP4, QTL for KPS on 1A and for SPM on 4D were presented in 17CVI, QTL for KPS on 4D and SPM on 6B were presented in 12BD, QTL for KPS on 7B was presented in 17DMS where QTL for yield on 1A at 412 Mbp was presented, we summarized the means and standard deviations of KPS, SPM and yield in these four environments. However, some trait in certain environments did not show many differences, such as SPM and KPS in 12BD, and KPS in 17 CVI. All 21 groups were in four clusters, five with 6B QTL, five with both 6B and 7B QTL, six with 7B QTL, and five with neither 6B nor 7B QTL. Within each cluster, the number of QTL for SPM were arranged as increasing while the number of KPS QTL as decreasing. For SPM in 17 CVI, in the first five groups with 6B that have one, two and three QTL, the trends of increasing SPM mean values could be observed but the differences were not significant in 12BD and 17CVI; however, the increasing trends were not observed in 13EP4 (S4A–S4C Fig). For KPS in the same five groups, the trends of decreasing could be observed in 12BD and 17DMS but not in 17CVI (S4D–S4F Fig). For yield from all three environments, the increasing trends could be seen in 17CVI and 17DMS but not in 12BD (S4G–S4I Fig). In the second cluster with six groups having 7B KPS QTL added, we can see the increasing trends on most groups of SPM and decreasing trends of KPS in 12BD and 17CVI, as well as decreasing trends of yield in all three environments. In the third cluster of six groups with 7B KPS QTL, decreasing trends of KPS only were observed in all three environments, 12BD, 17CVI and 17DMS (S4D–S4F Fig). In the fourth cluster with neither 6B nor 7B QTL, decreasing trends were observed on KPS and yield for most groups in 12BD, 17CVI, and 17DMS (S4D–S4I Fig). In summary of the four QTL in 1A and 4D, compensation effects were observed where the groups with more SPM QTL had higher values of SPM and those groups with more KPS QTL had higher KPS values in most cases. However, this only happened in some environments. In order to observe the effect of 6B and 7B, the groups with both, one or none of these two QTL were compared. For SPM in 17CVI, the groups were rearranged (S4J–S4L Fig). It was very clear to tell that 6B QTL increased SPM and had the similar values as those of groups with both 6B and 7B QTL if groups already had 1A SPM and 4D KPS QTL or 1A and 4D SPM (S4J Fig). on the other hand, for KPS in 17DMS, 7B QTL increased KPS and had similar values as the groups with both 6B and 7B if the groups already had 1A KPS QTL, or 1A SPM and 4D KPS QTL, or 1A and 4D SPM QTL (S4K Fig). After the yield values from 17CVI were arranged in the same order of groups, yield values from groups with only 6B or 7B or both 6B and 7B were very similar (S4L Fig). The yield values from 17DMS had the similar trends as that of 17CVI. In this case, the compensation and complementary effects of SPM and KPS for yield could be observed for these six QTL on chromosome 1A, 4D, 6B and 7B with favorable alleles from both TAM 111 and TAM 112. However, no epistasis interactions were not found for these six QTL (S5 Table), which is good for breeding application.

Epistasis, epistasis-by-environment, and additive-by-environment interactions

Only those with overall LOD scores > 5.0 were summarized for the epistasis and additive-by-environment interactions (S5 Table). Among 375 interactions for yield, only 56 had overall LOD scores > = 10.0, but none of the epistasis and additive-by-environment interactions had LOD > 10.0 (S5 Fig). Among 28 interactions that increased yield by more than 10 g m-2, there were six additive-by-environment interactions at 17BI with favorable alleles from TAM 111 that increased yield from 10.4 to 17.5 g m-2; two additional additive-by-environment interactions at 14EP4 and 17CVI, respectively increased yield by 10.2 and 10.4 g m-2 with favorable alleles from TAM 112. Among 19 epistasis-by-environment interactions, 17 interactions at 17BI increased yield by 10.1 and 13.1 g m-2 with seven favorable alleles from TAM 112 and 10 favorable alleles form TAM 111 (S5 Table and S5 Fig).

Among 234 interactions for spikes m-2 with LOD > 5.0, only eight interactions had overall LOD score > = 10.0 and no epistasis LOD > = 10.0. Five epistasis-by-environment interactions occurred in ME2 (including12BD, 13EP4, 13EP5, and 17CVI) in which three of them had favorable alleles from TAM 112 and increased spikes m-2 by 10.8 and two increased spikes m-2 by 14.4 with alleles from TAM 111 (S5 Table and S5 Fig).

For thousand kernel weight (TKW), among 581 interactions with overall LOD > = 5.0, there were 123 with overall LOD > 10.0 and 26 with epistasis LOD > = 10.0. There were five epistasis that increased TKW by 0.4–0.7 g with four having favorable alleles from TAM 112. However, only two of the five had epistasis LOD > = 10.0. Two of the five occurred within mega-environment ME2 including 14EP4, 14EP5, 17BI, 17CVI, and 17DMS. Among the eight additional interactions that increased TKW by 0.4 to 0.7 g, four were epistasis-by-environment interactions with all favorable alleles from TAM 112 and three of the four occurred in 17BI while the four additive-by-environment interactions occurred in 17CVI and 11BD with three having favorable alleles from TAM 111 (S5 Table and S5 Fig).

Among 243 interactions with LOD > 5.0 for kernels spike-1, only four had overall LOD > = 10.0 but none of them could increase the trait by > 0.4. Among six interactions that increased kernels spike-1 by 0.4, four epistasis-by-environment interactions increased kernels spike-1 by 0.4–1.0 with two having favorable alleles from TAM 111. The one increased by 1.0 had favorable alleles from TAM 111 at 13EP5 while the same interaction increased 0.4 with favorable allele from TAM 112 in 17DMS (S5 Table and S5 Fig).

For harvest index, there were 240 interactions that had LOD > 5.0 but only one interaction had overall LOD > 10.0. Four additive-by-environment and six epistasis-by-environment interactions increased harvest index by 0.5–0.8% at 17CVI with eight having favorable alleles from TAM 111.

For biomass grain yield, 190 interactions had LOD > 5.0 but only one had LOD > 10.0. All 16 epistasis-by-environment interactions at 17BI increased biomass grain yield by 15.1 to 19.6 g m-2 with ten having favorable alleles from TAM 112 and six having favorable alleles from TAM 111 (S5 Table). For total dry biomass, no interaction had LOD > = 10.0. Only four interactions explained 10.3–12.2% of total phenotypic variations but none of epistasis and epistasis-by-environment interactions increased more than 10 g m-2 (S5 Table).

Discussion

Evaluation of yield and yield component in individual and mega-environments

Yield is a complex trait affected by genetic, environment and genetic-by-environment interactions. Management in crop growing conditions, such as drought or irrigated, can also interfere with grain yield. Therefore, yield trials from multiple years at multiple locations are crucial to provide data of yield and yield components under various weather and management conditions including dryland and irrigated, and further lead to more reliable genetic analysis for yield plasticity [9]. In this study, we used an alpha lattice experimental design to conduct the trials in five growing seasons and up to five locations, which provided diverse growing conditions to evaluate yield and yield-related traits and thus being able to detect effects due to genetic and genetic-by-environment interactions. Through combined ANOVA and heritability analyses, trait data showed genetic variance at a significance level with heritability higher than 0.05 were used for QTL analysis. In addition, Pearson’s product moment correlation was conducted among all traits, and most of the correlations are significant, which was further supported by co-localized QTL linked to yield and yield components to indicate the presence of pleiotropy in genomic regions modulating the quantitative traits [28], and the positive correlation thus suggests a possible linkage existing in coupling phase or presence of positive pleiotropic effects [28].

Mega-environments (MEs) have been initially defined by CIMMYT as similar biotic and abiotic stress, cropping system requirements, and environments conditions by a volume of production [29]. Besides individual environment QTL analysis from genome-wide scan in this study (S2 Fig), QTL analyses across all individual environments (S3 Fig), and within mega-environments were also conducted, which minimized environment effects within MEs (S4 Table). This also increased the accuracy to identify a potential major QTL under mega environment and they are very important for local adaptation.

Compensation effects of QTL for spikes m-2 and kernels spike-1

Since the data were collected in multiple environments and QTL were identified only in certain environments, the QTL effects were compared based on data from QTL presenting environments could be difficulty. Furthermore, the summary of the QTL effects only explained a portion of the total phenotypic variations. For example, three QTL on chromosome 1A, 4D and 7B only explained up to 47.2% total phenotypic variations and 24.2% additive phenotypic variations while the three QTL for SPM on chromosome 1A, 4D and 6B only explained 37.9% total and 8.9% additive phenotypic variations (Table 1). Compensation effects of the six QTL for SPM and KPS with favorable alleles from both TAM 111 and TAM 112 were compared (S4A–S4L Fig). Pyramiding of various favorable alleles from different QTL would increase the trait values including overall grain yield but may be only effective for certain environments and genetic backgrounds with various combinations of QTL.

Yield component QTL comparisons

Several QTL studies have been conducted in wheat to identify closely linked markers associated with yield and yield components which can be applicable in marker-assisted selection and therefore improve the efficiency of breeding efforts. However, it is well-known that most of yield related traits are highly dependent on the interactions of genotype by environment. Therefore, in this study, we presented a set of 36 unique consistent QTL identified to be associated with one trait but from at least two out of the four analyses from individual, across all individual, within and across defined mega environments. Another large QTL cluster located on chromosome 7D ranging from 64 Mb to 68 Mb was identified in several environment analysis (Table 1) and has a confounding pleiotropic effect (Table 2) but it is far away from the grain size gene TaGS3-D1 that is at 6.5–6.8 Mbp [30, 31]; however, two markers, wsnp_Ra_c6894_11980338 at 57.0 Mbp and Excalibur_c22419_460 at 67.0 Mbp, closely linked with a major QTL on 7D for flour yield [29–32], overlapped with the QTL cluster in this study; however, other three QTL on chromosomes 3D at 414.6–552.0 Mbp, 4B at 139.7–407.4 Mbp and 6B at 697.5–712.4 Mbp were different from what were found from this study based on the linked SNPs that were at 574 Mbp on 3D, 659 Mbp on 4B, and 674 Mb on 6B (Table 1 and Fig 1). The TaGW2 was located at 237.8–240.3 Mbp based on the IWGSC RefSeq V1.0 and V 2.0 [33–35], which is far away from those QTL mapped on chromosomes 6A, 6B and 6D (Table 1 and Fig 1).

Since there is a common parent, TAM 111, in this study and previously published QTL from CO 960293-2/TAM 111 (CT111) by Assanga et al. [19], comparisons of QTL locations based on the linked SNP s associated with the yield and yield components were conducted. The two QTL at 13.8 and 411.7 Mbp on chromosome 1A associated with KPS and SPM were different from those two QTL in CT111 that were located around 462.4 and 471.7 Mbp on 1A. The major QTL associated with increasing yield and thousand kernel weight with favorable alleles from TAM 111 at 56.8–65.1 Mbp on chromosome 2BS was not found in this study. The increasing yield QTL at 12.4 Mbp on chromosome 6A was different from the QTL in CT111 at 176.6–185.2 Mbp and TaGw2-6A at around 237.0 Mbp. The QTL at 673.8 Mbp associated with SPM was different form the one in CT111 at 557.4 Mb (Table 1 and Fig 1).

Comparing the physical location of QTL for yield and yield components based on linked SNPs in a population derived from ND 705 and PI 414566, the QTL at 337.7 Mbp on 1A, 14.8 Mbp in 1B, 531.0 and 603.4 Mbp on 3D, 13.4 Mbp on 4B, 499.1 Mbp on 4D, 46.8 Mbp on 6B and 131.7 Mbp on 6D were different from the QTL we identified in this study [36]. Through a QTL mapping in a population derived from “Attila” and “CDC Go”, QTL for yield and kernel weight at 584.1 Mbp on 1B, 454.7 and 496.4 Mbp on 6A, 26.6 and 56.0 Mbp on 6B [37] were different from those consistent QTL in this study (Table 1 and Fig 1).

Dissection of QTL by environment, epistasis and additive-by-environment interactions

Some QTL were very significant for the total LOD score but not for the additive effect LOD scores LOD(A) (Tables 1 and 2, S4 Table). For example, Qhi.tamu.4D.29 had LOD(A) of 4.6 among total LOD of 16.9 and the total explained phenotypic variations by the QTL additive effect was 2.5% compared with 13.0% explained by the additive-by-environment. The total additive effect for harvest index was only increased by 0.2% when it was analyzed across all individual environments (S4D Fig). On the other hand, when analyzed within ME1 including 14EP4, 17BI, 17 CVI, and 17DMS, all irrigated or with high rainfall in that year, the same QTL had LOD(A) of 13.3 from total LOD of 14.1, and variation explained by additive effects was 15.2% from total 19.0% and increased harvest index by 0.6% (S4 Table). The same QTL had additive effect that increased harvest index by 1.1% at 17CVI. This is the advantage of dissecting the additive effects from the additive-by-environment interactions to identify the major QTL with higher additive effects but less additive-by-environment interactions. Among the 75 epistasis and epistasis-by-environment interactions that had LOD > = 10 or the interaction effects increased the traits more than those of most major QTL (10.0 g m-2 for yield and biomass yield, 0.4 g for thousand kernel weight, 0.4 for kernels spike-1, 10 for spikes m-2, 0.5% for harvest index, 15 g m-2 for biomass yield), only six out of 87 significant QTL involved with the epistasis-by-environment interactions (Tables 1 and 2, and S5 Table). They were Qyld.tamu.1B.376, Qtkw.tamu.2D.531, Qhi.tamu.4D.29, Qtkw.tamu.4D.409, Qyld.tamu.6A.12, and Qkps.tamu.7B.19, which can be a warning for breeding selection. Since breeders can only fix the additive effects by selection, through these analyses, breeders can have a better idea for what QTL are worthy of consideration for selection in breeding practice.

Conclusion

In this study, the wheat 90K Infinium iSelect SNP array and whole genome ddRADseq were used in the construction of high-saturated genetic map for QTL mapping associated with yield and yield components collected from 11 environments across five years and five locations across Texas and New Mexico in the US Southern High Plains. QTL were analyzed using single trait in single environment, single trait across multiple environments, and single trait within and across mega-environments in which lines performed similarly. In addition to additive effects, the interactions of additive-by-environment, epistasis and epistasis-by-environment were dissected. Among 87 significant QTL for nine traits, 36 consistent QTL were identified with presence in at least two above-mentioned analyses and ten pleiotropic QTL were found associated with more than one trait. The eight consistent and pleiotropic QTL were located at 411.7 Mb on chromosome 1A, at 20.6, 26.0, 445.5, 455.3 and 455.8 Mb on chromosome 4D, at 19.6 Mb on 6D and at 64.3 Mb on 7D. They increased dry biomass by 12.1 g m-2, harvest index by 0.6%, thousand kernel weight by 0.9 g with favorable alleles from TAM 112 and increased biomass grain yield by 16.3 g m-2, kernels spike-1 by 0.7, and yield by 20.3 g m-2 with favorable alleles from TAM 111. Only six of 75 epistasis-by-environment interactions were involved with the major QTL. Major QTL with larger additive effects and less interaction effects were identified and their compensation effects of various QTL from both parents were discussed.

Supporting information

[Figure omitted. See PDF.]

S1 Fig.

A. Boxplot analysis of yield and yield component traits. B. Dendrogram and biplot analysis of yield and yield component traits to classify mega-environments for each trait.

https://doi.org/10.1371/journal.pone.0237293.s001

(PDF)

S2 Fig. LOD profile and additive effects of QTL detected for all traits in each of the 11 environments.

https://doi.org/10.1371/journal.pone.0237293.s002

(PDF)

S3 Fig. Whole genome significance LOD(A) and LOD(AbyE) profiles of quantitative trait loci for yield and yield components based on across all the environments for each trait.

https://doi.org/10.1371/journal.pone.0237293.s003

(PDF)

S4 Fig. QTL allele effects on spike m-2 and kernels spike-1 from various environemnts with the presence of significant QTL.

https://doi.org/10.1371/journal.pone.0237293.s004

(PDF)

S5 Fig. Whole genome significance profiles of epistasis at LOD > 10 for yield and its components.

https://doi.org/10.1371/journal.pone.0237293.s005

(PDF)

S1 Table. The combined ANOVA, heritability, and mean performance for all traits across environments.

https://doi.org/10.1371/journal.pone.0237293.s006

(XLSX)

S2 Table. Genetic correlation among grain yield and other traits based on overall means and means measured in individual environments.

https://doi.org/10.1371/journal.pone.0237293.s007

(XLSX)

S3 Table. Mapped SNPs on 21 chromosomes, their genetic and physical length.

https://doi.org/10.1371/journal.pone.0237293.s008

(XLSX)

S4 Table. Significant SNPs for yield and yield components in all analysis based on 1000 permutation threshold.

https://doi.org/10.1371/journal.pone.0237293.s009

(XLSX)

S5 Table. Epistases and epistasis by environment interactions for yield and yield components in all analysis.

https://doi.org/10.1371/journal.pone.0237293.s010

(XLSX)

Acknowledgments

We thank the technical support from Lisa Garza, Hangjin Yu, Cody Shachter, Srirama Reddy and Chor Tee Tan from the wheat genetic program and graduate students Silvano Assanga and Bharath Reddy from Department of Soil and Crop Science, TAMU.

Citation: Yang Y, Dhakal S, Chu C, Wang S, Xue Q, Rudd JC, et al. (2020) Genome wide identification of QTL associated with yield and yield components in two popular wheat cultivars TAM 111 and TAM 112. PLoS ONE 15(12): e0237293. https://doi.org/10.1371/journal.pone.0237293

References

1. Wu X, Chang X, and Jing R. Genetic insight into yield-associated traits of wheat grown in multiple rain-fed environments. PloS ONE 2012;7: 31249. pmid:22363596

2. Barton NH, and Keightley PD. Understanding quantitative genetic variation. Nature Rev. Genet. 2002;1: 11–21.

3. Goldringer I, Brabant P, and Gallais A. Estimation of additive and epistatic genetic variances for agronomic traits in a population of doubled-haploid lines of wheat. Heredity 1997;79.

4. Kumar N, Kulwal PL, Balyan HS, and Gupta PK. QTL mapping for yield and yield contributing traits in two mapping populations of bread wheat. Mol. Breed. 2007;19: 163–177.

5. Snape JW, Foulkes MJ, Simmonds J, Leverington M, Fish LJ, Wang Y, et al. Dissecting gene × environmental effects on wheat yields via QTL and physiological analysis. Euphytica 2007;154: 401–408.

6. Zhang K, Tian J, Zhao L, and Wang S. Mapping QTLs with epistatic effects and QTL× environment interactions for plant height using a doubled haploid population in cultivated wheat. J. Genet. Genomics 2008;35: 119–127. pmid:18407059

7. Reif JC, Maurer HP, Korzun V, Ebmeyer E, Miedaner T, and Würschum T. Mapping QTLs with main and epistatic effects underlying grain yield and heading time in soft winter wheat. Theor. Appl. Genet. 2011;123: 283–292. pmid:21476040

8. Patil RM, Tamhankar SA, Oak MD, Raut AL, Honrao BK, Rao VS, et al. Mapping of QTL for agronomic traits and kernel characters in durum wheat (Triticum durum Desf.). Euphytica 2013;190: 117–129.

9. Grogan SM, Anderson J, Baenziger PS, Frels K, Guttieri MJ, Haley SD, et al. Phenotypic Plasticity of Winter Wheat Heading Date and Grain Yield across the US Great Plains. Crop Science 2016;56: 2223–2236.

10. Torada A, Koike M, Mochida K, and Ogihara Y. SSR-based linkage map with new markers using an intraspecific population of common wheat. Theor. Appl. Genet. 2006;112: 1042–1051. pmid:16450184

11. Collard BC, and Mackill DJ. Marker-assisted selection: an approach for precision plant breeding in the twenty-first century. Phil. Trans. R. Soc. 2008;1491: 557–572. pmid:17715053

12. Jannink JL, Lorenz AJ, and Iwata H. Genomic selection in plant breeding: from theory to practice. Brief. Funct. Genomics. 2010;9: 166–177. pmid:20156985

13. Wicker T, Mayer KF, Gundlach H, Martis M, Steuernagel B, et al. Frequent gene movement and pseudogene evolution is common to the large and complex genomes of wheat, barley, and their relatives. The Plant Cell 2011;23: 1706–1718. pmid:21622801

14. Cavanagh CR, Chao S, Wang S, Huang BE, Stephen S, Kiani S, et al. Genome-wide comparative diversity uncovers multiple targets of selection for improvement in hexaploid wheat landraces and cultivars. Proc. Nat. Acad. Sci. 2013;20: 8057–8062. pmid:23630259

15. Wang S, Wong D, Forrest K, Allen A, Chao S, et al. Characterization of polyploid wheat genomic diversity using a high‐density 90 000 single nucleotide polymorphism array. Plant Biotechnol. J. 2014;12: 787–796. pmid:24646323

16. Liu S, Ocheya S, Dhakal S, Gu X, Tan CT, Yang Y, et al. Validation of chromosomal locations of 90K array SNP in US wheat. Crop Sci. 2016;56: 364–373.

17. Lazar MD, Worrall WD, Peterson GL, Fritz AK, Marshall D, Nelson LR, et al. Registration of ‘TAM 111’ Wheat Registration by CSSA. Crop Sci. 2004;44(1):355–6.

18. Rudd JC, Devkota RN, Baker JA, Peterson GL, Lazar MD, Bean B, et al. ‘TAM 112’ Wheat, Resistant to Greenbug and Wheat Curl Mite and Adapted to the Dryland Production System in the Southern High Plains. J. Plant Regis. 2014;8: 291–297.

19. Assanga SO, Fuentealba M, Zhang G, Tan C, Dhakal S, Rudd JC, et al. Mapping of quantitative trait loci for grain yield and its components in a US popular winter wheat TAM 111 using 90K SNPs. PLoS ONE. 2017;12: 0189669.

20. Liu S, Rudd JC, Bai G, Haley SD, Ibrahim AM, Xue Q, et al. Molecular markers linked to important genes in hard winter wheat. Crop Sci. 2014;54: 1304–1321.

21. Dhakal S, Tan CT, Anderson V, Yu H, Fuentealba MP, Rudd JC, et al. Mapping and KASP Marker Development for Wheat Curl Mite Resistance in ‘TAM 112’ Wheat Using Linkage and Association Analysis. Mol. Breed. 2018;38: 119.

22. Yang Y, Basnet BR, Ibrahim AM, Rudd JC, Chen X, et al. Developing KASP Markers on a Major Stripe Rust Resistance QTL in a Popular Wheat TAM 111 Using 90K Array and Genotyping-by-Sequencing SNPs. Crop Sci. 2018;59: 165–175.

23. Elshire RJ, Glaubitz Jeffrey C., Qi Sun, Poland Jesse A., Ken Kawamoto, Buckler Edward S., Mitchell SE. A Robust, Simple Genotyping-by-Sequencing (GBS) Approach for High Diversity Species. PLoS ONE. 2011;6(5): e19379. pmid:21573248

24. Poland J, Endelman J, Dawson J, Rutkoski J, Wu S, Manes Y, et al. Genomic Selection in Wheat Breeding using Genotyping-by-Sequencing. The Plant Genome. 2012;5(3):103–13.

25. Ballare KM, Pope NS, Castilla AR, Cusser S, Metz RP, Jha S. Utilizing field collected insects for next generation sequencing: Effects of sampling, storage, and DNA extraction methods. Ecology and Evolution. 2019;9(24):13690–705. pmid:31938475

26. Alvarado G, López M, Vargas M, Pacheco Á, Rodríguez F, Burgueño J, et al. 2015. META-R (Multi Environment Trail Analysis with R for Windows) Version 6.04. CIMMYT Research Data & Software Repository Network, V23. https://hdl.handle.net/11529/10201

27. Meng L, Li H, Zhang L and Wang J. QTL IciMapping: Integrated software for genetic linkage map construction and quantitative trait locus mapping in biparental populations. Crop J 2015;3: 269–283.

28. Mackay TFC, Stone EA, Ayroles JF. The genetics of quantitative traits: challenges and prospects. Nature Rev. Genet. 2009;10: 565–577. pmid:19584810

29. Rajaram S, Van Ginkel M, and Fischer RA. CIMMYT’s wheat breeding mega- environments (ME). In Proc. Intl. Wheat Genet. Symp. 1994;1101–1106.

30. Rasheed A, Wen W, Gao F, Zhai S, Jin H, Liu J, et al. Development and validation of KASP assays for genes underpinning key economic traits in bread wheat. Theor. Appl. Genet. 2016; 129(10): 18. pmid:27306516

31. Zhang Y, Liu J, Xia X and He Z. TaGS-D1, an ortholog of rice OsGS3, is associated with grain weight and grain length in common wheat. Mol. Breed. 2014; 34: 1097–1107.

32. Cabral A. L., Jordan M. C., Larson G, Somers D. J., Humphreys D. G., and McCartney C. A. Relationship between QTL for grain shape, grain weight, test weight, milling yield, and plant height in the spring wheat cross RL4452/‘AC Domain’. PLoS ONE, 2018; 13(1), e0190681. pmid:29357369

33. Zhang Y., Li D, Zhang D, Zhao X, Cao X, Dong L, et al. Analysis of the functions of TaGW2 homoeologs in wheat grain weight and protein content traits. The Plant Journal, 2018; 94(5): 857–866. pmid:29570880

34. Qin L, Zhao J, Li T, Hou J, Zhang X and Hao C. TaGW2, a Good Reflection of Wheat Polyploidization and Evolution. Frontiers in Plant Science 2017; 8(318). pmid:28326096

35. Su Z, Hao C, Wang L, Dong Y and Zhang X. Identification and development of a functional marker of TaGW2 associated with grain weight in bread wheat (Triticum aestivum L.). Theor. Appl. Genet. 2011; 122(1): 211–223. pmid:20838758

36. Kumar A., Mantovani E. E., Simsek S., Jain S., Elias E. M. and Mergoum M. Genome wide genetic dissection of wheat quality and yield related traits and their relationship with grain shape and size traits in an elite × non-adapted bread wheat cross. PloS ONE. 2019; 14(9): e0221826. pmid:31532783

37. Yang RC, Randhawa H and Spaner D. Mapping QTLs Controlling Agronomic Traits in the ‘Attila’ × ‘CDC Go’ Spring Wheat Population under Organic Management using 90K SNP Array. Crop Sci. 2017; 57(1): 365–377.