Plant materials
The Hard Winter Wheat Association Mapping Panel (HWWAMP) is a population of 258 hard red and 41 hard white winter common wheat inbred lines assembled by public and private breeding programs reflecting the wide precipitation and temperature gradients across the primarily dryland wheat cropping systems of the Great Plains region of the United States (Guttieri et al. 2015). The lines were previously categorized into northern (Montana, North Dakota, South Dakota, n = 39), central (Colorado, Kansas, Nebraska, n = 119), and southern (Oklahoma, Texas, n = 105) classes, based on the geographical region within the Great Plains where the lines were developed (Grogan et al. 2016b). A fourth category of “other” (n = 36) was used for lines that did not fit one of these three regional categories.
The CO940610/’Platte’ hard white winter wheat recombinant inbred line (COP-RIL) population was previously described (Dao et al. 2017). Remnant F5:6 seed from this population was advanced to produce F5:7 and F5:8 seed used in field experiments. CO940610 (GSTR 10702; pedigree - KS87H22/MW09) is an experimental hard white winter wheat line developed by the Colorado State University wheat breeding program and carries Ppd-B1b, Vrn-D3a, Wapo-A1a, and QSn.csu-6Ba alleles. ‘Platte’ (PI 596297; pedigree - ‘Tesia 79’/Chat ‘S’//‘Abilene’) is a hard white winter wheat cultivar developed by HybriTech Seed International and carries Ppd-B1a, Vrn-D3b, Wapo-A1b, and QSn.csu-6Bb alleles.
Eight HIFs were derived from COP-RIL F3:4 remnant seed by screening 528 individual plants with Kompetitve Allele Specific PCR (KASP) marker BobWhite_c22638_135 (Table S1). The 35 individual plants heterozygous for this marker were genotyped with three other KASP markers (KS0619_760354, KS0619_760357, and KS0619_760427) spanning the QSn.csu-6B locus. Individuals heterozygous for all four loci were then genotyped using KASP markers to distinguish WAPO-A1, VRN-D3 and PPD-B1 alleles (Table S1). From the 35 individuals, eight HIFs at the F4:5 generation were derived, each with a different combination of fixed PPD-B1, WAPO-A1 and VRN-D3 alleles (Table S2). From each HIF, 96 individual F5 plants were genotyped for the QSn.csu-6B locus with the same four KASP markers to identify single plants fixed for either the QSn.csu-6Ba or QSn.csu-6Bb allele. F5:6 seeds were derived from these individuals and used in field experiments.
Field experiments
Field experiments for the HWWAMP were planted in Greeley, Colorado (40.447, 104.636 W; elevation 1,425 m; soil type Nunn clay loam and Olney fine sandy loam) under irrigated and dryland conditions in the fall of 2011, referred to as “Greeley 2012 irrigated” and “Greeley 2012 dryland”, respectively. The population was planted again in Fort Collins, Colorado (40.648 N, 104.993 W; elevation 1,558 m; soil type Nunn clay loam) under irrigated and dryland conditions in the fall of 2012, referred to as “Fort Collins 2013 irrigated” and “Fort Collins 2013 dryland”, respectively. Plots approximately 3.6m2 were arranged in a row-column design with two check varieties ‘Settler CL’ (PI 655242; pedigree – ‘Wesley’ sib//’Millennium’ sib/’Above’ sib) and ‘Hatcher’ (PI 638512; pedigree - ’Yuma’/PI 372129//’TAM-200’/3/4*Yuma/4/KS91H184/’Vista’). The experimental genotypes were unreplicated for each environment.
Field experiments for the biparental populations were planted in Fort Collins, Colorado. The F5:6, F5:7, and F5:8 COP-RIL populations were planted under irrigated conditions in the fall of 2016, 2018 and 2019, referred to as “Fort Collins 2017”, “Fort Collins 2019”, and “Fort Collins 2020”, respectively. Each line was planted in a two-row plot 0.92 m long with 23 cm spacing between rows, 28 cm spacing between plots, and a planting density of approximately 2,500,000 seeds ha-1. In Fort Collins 2017, two replicates were planted in a randomized complete block design. In Fort Collins 2019 and Fort Collins 2020, one replicate was planted where the RILs were randomized. Eight F5:6 HIFs were planted in Fort Collins 2020 under the same conditions as the RILs. Each HIF was planted in three replicates using a randomized complete block design with a family block nested within each replicate block. Each replicate block consisted of all eight family blocks and each family block consisted of 18 lines with nine lines fixed for QSn.csu-6Ba and nine lines fixed for QSn.csu-6Bb.
Phenotyping
Kernel width was measured in the HWWAMP in four environments; Greeley 2012 dryland, Greeley 2012 irrigated, Fort Collins 2013 dryland, and Fort Collins 2013 irrigated. Grain length, width, and weight were measured for the COP-RIL population in Fort Collins 2017. Each sample consisted of ten spikes randomly selected from each plot. The grain was threshed using a wheat head thresher (Precision Machine Co., Inc, Lincoln, NE, USA, part # WHTA0100001). Broken grains were discarded and approximately 400 grains were spread evenly across a bed scanner (Microtek International Inc, Hsinchu City, Taiwan, part # MRS-3200A3L). Scanned images were collected at a vertical and horizontal resolution of 300 dots per inch and analyzed using the software GrainScan (Whan et al. 2014) to calculate mean length, width and grain number. Sample weight was divided by grain number to calculate mean grain weight. Grain weight, heading date, and spikelet number data were previously collected for the HWWAMP in the same four environments used for kernel width (Grogan et al. 2016a), and these phenotypic data were downloaded from the T3/Wheat database (Blake et al. 2016).
Spikelet number was measured in each environment by manually counting spikelets beginning from the first rachis node to the terminal spikelet. Only heads with an intact terminal and basal spikelet were counted. For Fort Collins 2017, ten spikes were randomly selected from each plot and grain harvested from those spikes were used to measure grain length, width and weight. Twenty spikes were selected for Fort Collins 2019 and Fort Collins 2020.
Heading date was recorded in Julian days from January 1 when approximately 50 % of the spikes in a plot had fully emerged from the flag leaf sheath. Heading date was recorded for each line of the COP-RIL and HIF populations in Fort Collins 2020. All raw phenotypic data is provided in Supplemental data 1.
Genotyping
Lines from the HWWAMP population were previously genotyped using the 90K SNP array (Wang et al. 2014) and data were downloaded from the T3/Wheat database (Blake et al. 2016). Markers with a minor allele frequency of greater than 0.05 and missing in less than 20 % of genotypes were retained for downstream analysis (Supplemental data 1). Peak markers BobWhite_c22638_135 and IWA5913 from the 90K SNP array data were used to differentiate between the QSn.csu-6B and WAPO-A1 alleles, respectively. Marker data for PPD-B1 was generated using a PCR-based assay (Beales et al. 2007) and were described for the HWWAMP previously (Grogan et al. 2016b). VRN-D3 was genotyped using a cleaved amplified polymorphic sequence marker previously described (Chen et al. 2010).
DNA was extracted using a modified Cetyl Trimethyl Ammonium Bromide (CTAB) extraction protocol (Doyle and Doyle 1990). DNA samples were diluted in water to a final concentration between 50 and 100 ng/µL prior to genotyping. DNA concentration was measured via a spectrophotometer (Thermo Scientific, Waltham, MA, USA, NanoDrop 1000). KASP reactions were performed in a total volume of 10 µL comprising 5 µL KASP V4.0 2X Master Mix (LGC Genomics, LLC, Beverley, MA, USA, Cat. # KBS-1016-017-US), 0.14 µL of 54 µM KASP primer mix, and approximately 100 ng of DNA template. The PCR conditions used consisted of: 94°C for 15 min for hot-start Taq activation; ten cycles of two-step PCR with a 94°C denaturing step for 20 s and an annealing-elongation step for 60 s starting at 61°C and decreasing by 0.6°C each cycle until reaching 55°C; 26 further cycles of two-step PCR with a 94°C denaturing step for 20 s and a 55°C annealing-elongation step for 60 s. PCR was carried out on an Applied Biosystems QuantStudio Q3 qPCR machine (Applied Biosystems, Waltham, MA, USA, Cat. # A28136) and analyzed using QuantStudio Design and Analysis software (version 1.4). KASP primer sequences for the VRN-D3 locus were designed from a previously described cleaved amplified polymorphic sequence marker and provided by the United States Department of Agriculture, Agricultural Research Service, Hard Winter Wheat Genetics Research Unit (Manhattan, KS, USA) (Table S1). KASP primer sequences for the WAPO-A1 locus were previously published (Table S1). KASP primer sequences for the QSn.csu-6B alleles were designed using the online primer design tool PolyMarker (Ramirez-Gonzalez et al. 2015). The PPD-B1 locus was genotyped using a presence-absence marker adapted from the previously described PCR-based assay by using the same KASP PCR conditions, but with 42 µM of KASP primer mix with only one forward primer including a VIC fluorophore tail sequence (Table S1) (Beales et al. 2007). Primers were synthesized by Integrated DNA Technologies (San Diego, CA, USA). All genotypic data is provided in Supplemental data 1.
Statistical analysis
Grain width best linear unbiased estimates (BLUEs) were calculated for the HWWAMP with the lmer function of the afex R package version 0.28-1 (Singmann et al. 2020) using a mixed model where environment and the row-column location of each plot was treated as a random effect and genotype was treated as a fixed effect. The BLUEs model also accounted for the nested structure of plot location within each environment. Genome-wide association mapping was conducted upon 16,058 markers and 298 genotypes using a mixed model method previously described (Yu et al. 2006) and R version 4.0.2 (Team R Development Core 2020) with the rrBLUP version 4.6.1 package (Endelman 2011). In the mixed model, population structure was treated as a fixed effect using an additive relationship matrix calculated from the markers, and genotype was treated as a random effect. Marker effect was treated as a fixed effect following the restricted maximum likelihood method previously described (Endelman 2011). QTL were identified as significant markers (-Log10 (P) >3) within 10.0 cM of one another, using the genetic map previously published for the 90K SNP array (Supplemental data 1) (Wang et al. 2014).
Spikelet number, heading date, and kernel weight BLUEs for the HWWAMP were calculated with a linear model considering both genotype and environment as fixed effects due to data from the T3/Wheat database being corrected prior to this analysis. BLUEs calculated for the COP-RIL population used the same method, however the Fort Collins 2017 data was averaged across two replicates prior to analysis to avoid confounding estimates due to the unbalanced nested structure of two replications in Fort Collins 2017; no replication was included in Fort Collins 2019 and Fort Collins 2020.
Spikelet number and kernel weight data for the HWWAMP were analyzed via one-way analysis of variance (ANOVA) using R version 4.0.2 (Team R Development Core 2020) with the package emmeans version 1.5.0 (Lenth 2019). Only 249 lines that had genotype and phenotype data for all four loci of interest (WAPO-A1, VRN-D3, PPD-B1 and QSn.csu-6B) and four environments (Greeley 2012 dryland, Greeley 2012 irrigated, Fort Collins 2013 dryland, and Fort Collins 2013 irrigated) were used to calculate the effect sizes and P-values.
Pearson correlations were estimated using the Pearson’s Product-Moment Correlation Test from the stats R package (R version 4.0.2, stats package version 4.0.2). HWWAMP correlations for each environment were calculated between spikelet number, heading date and kernel weight for each environment and BLUEs across all environments. The Fort Collins 2017 COP-RIL correlations were calculated among kernel length, kernel width, kernel weight, and spikelet number. The Fort Collins 2020 COP-RIL correlations were calculated between heading date and spikelet number.
One-way ANOVA was conducted for all kernel length, kernel width, kernel weight, heading date, and spikelet number observations of the COP-RIL population with the loci of interest as the dependent variable using R version 4.0.2 (Team R Development Core 2020) with the package emmeans version 1.5.0 (Lenth 2019). Prior to analysis, COP-RIL individuals that were segregating at the WAPO-A1, VRN-D3, PPD-B1 and QSn.csu-6B loci were removed. To identify interactions affecting spikelet number, a full linear model was constructed accounting for the four-way interaction between the WAPO-A1, VRN-D3, PPD-B1 and QSn.csu-6B loci and environment. Akaike information criterion was calculated for all possible model subset combinations of the full model using the MuMin package in R (R version 4.0.2, MuMIn package version 1.43.17) (Barton 2020). The model with the lowest Akaike information criterion was the full four-way interaction model, which was used for subsequent analyses on interactions. Significance of the model terms was calculated via the Tukey-adjusted type II ANOVA method.
The HIF spikelet number data was analyzed via a linear model accounting for the allelic state of the QSn.csu-6B locus by family and the replicate effect using R version 4.0.2 with the package emmeans version 1.5.0. Since the WAPO-A1, VRN-D3, and PPD-B1 alleles were nested within each family, these variables were not included in the model. Using the backwards selection approach, a model that accounted for the effect of the QSn.csu-6B locus, family, and replicate was identified, which was used for all downstream analysis.