Genome-wide association study in historical and contemporary U.S. winter wheats identifies heightreducing loci

Sintayehu D. Da Priyanka Tyagi, Gina Brown-Guedirc, Mohsen Mohammadi*

aDepartment of Agronomy,Purdue University, West Lafayette,IN 47907,USA

bDepartment of Crop and Soil Sciences,North Carolina State University, Raleigh,NC 27695,USA

cU.S.Department of Agriculture,Agricultural Research Services,Southeast Area,Plant Science Research,Raleigh,NC 27695,USA

Keywords:GWAS Marker-trait associations Rht-B1 Rht-D1 Candidate genes Position-dependent procedures

A B S T R A C T

1. Introduction

Development of semi-dwarf and early flowering wheat cultivars was among the major achievements of the Green Revolution and led to large increases in grain yield [1,2].Introduction of semi-dwarfing genes improved lodging resistance, thereby increasing fertilizer response in cultivars and ultimately leading to increased kernel number and grain yield [1]. Two of the reduced-height (Rht) genes deployed worldwide are Rht-B1b (also known as Rht1) and Rht-D1b (also known as Rht2), which were derived from cultivar Norin 10 [2-4]. In Europe, Rht8, derived from Akakomugi, was also used widely. At least 20 Rht loci(https://wheat.pw.usda.gov/cgi-bin/GG3/browse.cgi?query=%2ARht%2A;class=gene;begin=1) have been reported to affect wheat plant height and are located on chromosomes 2A,2B,2D,4B,4D,5A,7B,and perhaps others.The dwarfing genes in wheat are either gibberellic acid-insensitive (GAI)or gibberellic acid-responsive (GAR). For instance, cultivars carrying the Rht-B1b and Rht-D1b mutant alleles do not respond to application of GA, whereas cultivars carrying GAR genes such as Rht8c and Rht13 in the absence of the Rht1 mutation respond to application of GA[1,5].

Agricultural drought has increased since the 1970s in most parts of the world [6], suggesting a critical need for cultivars with sufficient drought resilience to maintain or increase grain yield. Although Rht-B1b and Rht-D1b dwarfing alleles were integral components of the Green Revolution, they are associated with poor emergence and reduced early seedling growth if conditions are unfavorable [1,5,7]. Zanke et al. [8]reported that the wild-type alleles(conferring the tall phenotype)at Rht-B1 locus are associated with greater kernel weight.

At least 70% of current wheat cultivars worldwide carry Rht-B1,Rht-D1,and Rht8c[1,9].Forty-five percent of soft winter wheat cultivars carry Rht-D1b and 28% carry Rht-B1b [9]. Of 59 soft red winter wheat lines developed after 2000 at Purdue University,58 carried Rht-B1b and only one carried the Rht-D1b allele. In fact, the Rht-B1b and Rht-D1b alleles may already have been fixed in most breeding populations.Rebetzke et al.[10] demonstrated that GAR dwarfing genes such as Rht4,Rht12, and Rht13 increased kernel number and yield without compromising aerial biomass or coleoptile length in bread wheat. Identification of alternative loci for reduced height using allele mining and genetic studies of diverse germplasm would provide further useful genetic diversity for wheat improvement.

Genome-wide association studies (GWAS) provide opportunities to identify loci that control traits of interest. GWAS uses a collection of genetically diverse accessions that are genotyped using a dense molecular-marker platform followed by identification of genotype-phenotype relationships using marker-trait association(MTA)[11].The present study used a unique collection of germplasm developed before and after the Green Revolution, covering 202 years (1814-2015) of selection and breeding in the Midwestern region of the United States of America.The objectives of this study were to identify MTAs for plant height and possibly identify associated candidate genes nearby using the annotated wheat reference genome.

2. Materials and methods

2.1. Plant materials and phenotyping

The mapping panel of 260 historical and contemporary collections of soft red winter wheat accessions were planted for three years (2016-2018) at the Agronomy Center for Research and Education of Purdue University, West Lafayette,IN, USA. Accessions were planted in 1-m single-row plots spaced 25 cm apart. Plant height was recorded in all the three years. All 2017 records were discarded owing to moisture stress during that season; and the height records from 2016 and 2018 were used. The accessions are described in Table S1.

The field received 106 kg N ha-1of N fertilizer in both years. Accessions that were tall and prone to lodging were protected from lodging with metal poles and ropes until the end of the cropping season. The plant height datasets generated in 2016 and 2018 are referred to as PLH16 and PLH18, respectively. The combined dataset (LSMEANS) is referred to as PLH168. Changes in plant height over the breeding history were evaluated as trends over time by considering four groups classed by year of release: before 1920,1920-1960,1960-2000,and after 2000.

2.2. Genotyping

Genotyping was performed to generate single nucleotide polymorphism(SNP)markers using genotyping-by-sequencing(GBS).The details of the genotyping procedure,marker filtering,and Kompetitive allele specific PCR(KASP)assays are described in Ref. [12]. All SNP markers not assigned to any chromosome were removed. The remaining markers were filtered by minor allele frequency (MAF) ≥5% and missing values ≤30%. Missing genotypic data were imputed using the Linkage Disequilibrium K-Number Neighbor Imputation (LDKNNi) [13] algorithm implemented in TASSEL 5.0 [14]. In addition to the genome-wide SNP marker data,KASP markers for the known large-effect loci Rht-B1 and Rht-D1, were used, following Ellis et al. [15] and Rasheed et al.[16].

2.3. Genome-wide association study (GWAS)

GWAS was performed with the three phenotype datasets(PLH16, PLH18, and PLH168) and 38,695 markers including two KASP markers (Rht-B1 and Rht-D1) using TASSEL. A mixed linear model was fitted to perform full optimization and estimate population parameters for each SNP,resulting in the estimation of the genetic and error variances for each SNP by the Efficient Mixed-Model Association (EMMA) algorithm [17]. In GWAS, the first three principal components (PCs) accounted for a fraction of population structure, while a familial relatedness matrix estimated from marker data accounted for kinship.The association genetics model included markers and population structure as fixed effects and kinship and error as random effects. Manhattan plots were produced using the negative logarithm to base 10 of the P-values,abbreviated as –lg (P), using the R package qqman [18].Marker-trait associations with –lg (P) ＞4 (P-values ＜0.0001) were used for summary and those passing a false discovery rate (FDR) threshold of ＜0.001 were retained for further study. In a follow-up analysis using the combined data (PLH168), the Rht-B1 and Rht-D1 markers were incorporated as fixed effects to assess the effect of these large-effect loci on other potential associations. Genes in close proximity (±250 kb) to the major MTAs were identified using the IWGSC RefSeq v1.0 annotation v1.0,iwgsc_refseqv1.0_HighConf_2017Mar13.gff3.zip (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Annotations/v1.0/).

For identified loci, goodness of fit of a multi-locus linear regression of plant height on number of height-reducing alleles for the 16 MTAs was evaluated. The change in frequency of favorable alleles over the four year-of-release accession groups was also assessed.

3. Results

3.1. Phenotypic distribution and trend over breeding history

Plant height varied widely in individual years and across years(Table 1). Plant height in 2016 ranged from 64 cm to 155 cm,representing a 2.4-fold difference,with a mean of 122 cm.For 2018, the minimum and maximum plant heights were 73 cm and 137 cm, respectively, representing a 1.9-fold difference,with mean plant height of 108 cm. There was a mean plant height reduction of 14 cm in 2018 compared to 2016, which may be explained by the change of field location and seasonal condition.The 2016 field was mist-irrigated,whereas the 2018 crop was grown under rainfed conditions. The mean plant height in the combined data over the two years was 115 cm,with a minimum of 68 cm and maximum of 142 cm. A decreasing trend in plant height during the breeding and selection history of wheat,particularly after 1960,was evident(Fig. 1). This period coincides with the beginning of deployment of the Green Revolution genes Rht-B1 and Rht-D1[3,4].

The broad-sense heritability of height was high(H2= 0.82).This finding was in agreement with a relatively high correlation (r = 0.82, P ＜0.0001) between plant heights across the two years(Fig.S1).

3.2. Population structure

A 3D-plot of the first three principal components estimated from 38,695 markers (Fig. S2) shows three distinct clusters of accessions. Each cluster includes accessions released indifferent periods.Daba et al.[12]reported that the translocation from T.timopheevii representing by TaSus2-2B accounted for the clustering pattern. The stratification of the population may indicate a need to include population structure parameters in association analysis.

Table 1-Plant height (cm) of the 260 accessions in year 2016 (PLH16), year 2018 (PLH2018), and combined years(PLH168).

3.3. GWAS and characterizing the MTAs

For determining the number of PCs to use as population structure in GWAS,we assessed models including 0 to 5 PCs in the model to account for population structure. The P-P plots(Fig. S3) of the expected against observed P-values for the models showed deviations from the diagonal for the cases of 5 and 0 PCs in the model. The P-P plots generated by models with two, three, or four PCs showed nearly identical patterns and close to the diagonal.Therefore,we selected three PCs for presentation.

Over the three analyses, 60 SNP and two KASP markers were associated with plant height using values from two years (PLH16 and PLH18) and combined over the two years(PLH168) at -lg (P) ≥4.0, with 40 of them located on chromosome 2D in a stretch of about 2.4 Mb (Table S2). The correlation coefficients among these 40 markers ranged from 0.94 to 1.00 (Table S3), suggesting that these 40 SNP markers represent a single MTA locus. Therefore, we represented the MTAs on 2D by a single marker with high signal in all the three analyses.Table S4 lists the Chinese Spring alleles in the 40 marker-trait associations that represent this region on chromosome 2D.The 62 MTAs represented 16 loci distributed over 12 chromosomes:1A,2B,2D,3B,4B,4D,5A,5D,6A,6B,7A,and 7D(Fig.2A-C,Table 2).Two of the well-known genes(Rht-B1 and Rht-D1) affecting plant height were also among these MTA loci. Three of the MTAs (QPLH-2D, QPLH-4B.2, and QPLH-4D) were detected in all three analyses and remained significant after correction for multiple comparison using FDR (Padj≤0.001). Two other MTAs (QPLH-2B and QPLH-4B.1)were detected in two of the three analyses, but nonsignificant after FDR correction. Our discussion focuses on the three major MTAs.

The MTA with greatest association signal was recorded for QPLH-4D that represent Rht-D1, with –lg (P) values of 8.3 for PLH16, 12.4 for PLH18, and 12.7 for PLH168 (Table 2). The Rht-D1b allele (the height-reducing allele) was associated with a reduction in plant height of 15.9 cm in PLH18. On chromosome 4B, we detected an MTA QPLH-4B.2 that represents Rht-B1 with-lg(P)of 4.3 in PLH16,7.5 in PLH18,and 7.5 in PLH1618.The Rht-B1b allele at this locus led to a reduction of plant height by 11.2 cm in PLH18. The third strong signal was for QPLH-2D, with a -lg (P) of 5.4 for PLH16, 5.2 for PLH18, and 6.2 for PLH1618. The height-reducing allele at this MTA loci reduced plant height by 15.6 cm in PLH18.Details of the allelic forms in the mapping panel for the three major loci or MTAs detected in this study are presented in Table S1.

3.4. Cumulative effect and trend in favorable allele combinations

We assessed the cumulative effect of height-reducing allele combinations at the 16 MTAs.The number of height-reducing allele combinations in the panel ranged from 1 to 10, with 98.5%of the panel having 2-9 favorable alleles(Fig.3A).As the number of height-reducing alleles increased, plant height clearly decreased. For example, accessions with two heightreducing alleles (n = 11) had a mean plant height of 118 cm,whereas accessions with nine height-reducing alleles (n = 8)had a mean plant height of 95 cm,a difference of 23 cm.

It is expected that alleles targeted by breeding increase in frequency during breeding history. Assessing the change in frequency of alleles at a locus during selection and breeding could shed light on the history of trait selection in breeding programs. The frequency of height-reducing alleles at most loci increased over the four-year groups, particularly after 1960 (Fig. 3B). The exceptions were QPLH-2B and QPLH-4B.1,which showed a decreasing trend. Three of the high signal MTAs (QPLH-2D, QPLH-4B.2, and QPLH-4D) began to appear only after 1960.

Fig.1- Changes in plant height across four year-of-release accession groups. Plant height is based on Best linear unbiased estimate(BLUE)values over 2016 and 2018.

Fig.2-Manhattan plots representing–lg(P)values associated with each marker on the Y axis for year 2016(A),2018(B),and for the combined-year analysis(C).

Table 2-Identified marker-trait associations (MTAs) with their GWAS statistics for years 2016 (PLH16) and 2018 (PLH18),and the combined-year(PLH168)dataset.

3.5.Multiple regression of plant height(PLH168)on the MTAs

The multiple regression analysis of the 16 MTAs is shown in Fig. 4A. The two MTAs QPLH-4B.2 and QPLH-4D representing the widely-used height-reducing genes Rht-B1 and Rht-D1,respectively had the largest effect, followed by QPLH-2D,which also contributed greatly to height reduction. When predictability of plant height based on the 16 MTAs was assessed, the correlation between observed and predicted values was r = 0.83. However, when QPLH-2D, QPLH-4B.2, and QPLH-4D were excluded from the model, the correlation between the observed and predicted values dropped to r =0.67(Fig.4B).

3.6. Wheat genes located near MTAs

As expected, two of the MTAs (QPLH-4B.2 and QPLH-4D),representing Rht-B1 and Rht-D1, were located within 200 bp of genes TraesCS4B01G043100 and TraesCS4D01G040400,respectively, encoding GAI-like protein 1 (Table 3). The QPLH-2D locus could be of interest to breeders as an alternative to Rht-B1 and Rht-D1. The same locus could also be of interest to plant biologists interested in identifying candidate genes in this region. When we examined the wheat genome assembly near this MTA, we noticed a gene(TraesCS2D01G055700) that encodes GRAS transcription factor (GIBBERELLIC-ACID INSENSITIVE (GAI), REPRESSOR of GAI (RGA),and SCARECROW(SCR))located 105.8 kb from the MTA. BLASTX alignment of the TraesCS2D01G055700 coding sequence against the nonredundant protein database (http://blast.ncbi.nlm.nih.gov/) returned the refseq XP_020188264 protein accession in Aegilops tauschii subsp.tauschii encoding a 312 amino acid long scarecrow-like protein 3(SCL3)with 100%amino acid identity.

4. Discussion

Plant breeding programs worldwide intensively used the height-reducing genes Rht-B1 and Rht-D1, particularly after 1960. This period coincides with the era of the Green Revolution, when height-reducing genes were integral components of yield improvement [2]. The heightreducing genes helped to improve lodging resistance and increased assimilate partitioning to developing kernels,thereby increasing grain yield. Currently, over 70% of world wheat cultivars carry at least one of the three majoreffect genes Rht-B1,Rht-D1,and Rht8c[9,19].Guedira et al.[9]also reported that 61%of the cultivars developed in the U.S.Midwest(the states of Ohio,Indiana,Illinois,and Kentucky)carried Rht-B1b and 33% carried Rht-D1b. According to Tian et al. [19], Rht24 is the most frequent gene in Chinese elite wheat germplasm, with 76% of tested elite lines carrying Rht24. Our assessment of 59 Purdue breeding lines developed after 2000 revealed that 58 of them carried Rht-B1b and only one Rht-D1b.Thus,it is not surprising that plant height has decreased during the breeding and selection history of wheat, as shown in our trend analysis. The increase in the favorable allele frequency for widely used Rht genes(Rht-B1 and Rht-D1) after 1960 also supports the trend of height reduction over this period.In addition to these two loci,our study also showed a considerable increase of the favorableallele frequency after 1960 at almost all loci except QPLH-2B and QPLH-4B.1.This finding suggests that Rht-B1 and Rht-D1 were selected intentionally while the other loci were also selected unintentionally during the breeding and selection history of wheat.

Fig.3- The distribution of height-reducing alleles(left bar plot)along with cumulative effects of height-reducing alleles on plant height(the trend graph on the left bar plot)(A),and the trend in the frequency of height-reducing alleles at the identified MTAs across the four year-of-release groups(B).

Fig.4- Scatter plot of predicted plant height against observed plant height in the combined data(PLH168)using multi-locus regression using the 16 MTA loci as predictors(A)and excluding three high-signal MTAs(QPLH-2D,QPLH-4B.2,and QPLH-4D)from the regression model(B).

Optimizing plant height as a prime target in cereal breeding has reduced lodging and improved yield. Plant height remains an important practical consideration fordesigning and diversifying breeding populations in the future.Mapping studies for reduced height in wheat started as early as 1972, when Morris et al. [20] used monosomic analyses to map a gene for plant height to chromosome 4A. Since then,height-reducing genes have been mapped to all 21 wheat chromosomes[19].At least 24 Rht genes have been reported to date [19,21] (Table 4), and some of these genes are located in the same chromosomes as the MTAs we identified. These include Rht8 on 2D,Rht-B1 and Rht11 on 4B,Rht-D1 on 4D,Rht12 on 5AL,and Rht9 and Rht13 on 7B.Several QTL for plant height in wheat are also archived in T3/Wheat Toolbox (https://triticeaetoolbox.org/wheat/qtl/qtl_report.php; Table S5), with 43 QTL were shown to be significant at P ＜0.0001 and located on all chromosomes except 1A,2A,2D,3D,4D,5A,and 7B.Our study detected 16 genomic regions implicated in plant height,distributed on chromosomes 1A,2B,2D,3B,4B,4D,5A,5D,6A,6B,7A,and 7D.

Table 3-Some wheat genes with potential function in plant growth.

Some of our MTAs coincided with known genes reported previously. The KASP markers representing Rht-B1 and Rht-D1 were significantly associated with plant height, implying that our MTAs on 4B and 4D are represent the two widely-used genes.Bellucci et al. [22] reported two QTL on chromosome 2D, with one of these loci (QPht.cau-2D.2) marked by RAC875_c74_204,which is positioned at 22,352,745 bp. Gasperini et al. [23] also fine-mapped Rht8 between DG279 and xfdc53 on chromosome 2D. Applying Basic Local Alignment Search Tool (BLAST)analysis using the primer sequences for these markers, we located them at positions 22,202,306-22,202,325 bp and 23,024,744-23,024,763 bp, respectively (Fig. 5). This finding could imply that our MTA on chromosome 2D, QPLH-2D(positioned at 22,352,745 bp), is the same as QPht.cau-2D.2 as well as the height-reducing gene Rht8. Chai et al. [24] also detected two QTL (QPht/Sl.cau-2D.1 and QPht/Sl.cau-2D.2) that affect both plant height and spike length in the same genomic region as QPLH-2D.The GRAS transcription factor linked to QPLH-2D via TraesCS2D01G055700 was also located in the same region of 2D (22,246,914-22,248,492 bp), and was also reported previously by Chai et al.[24]as a possible candidate gene for Rht8.This protein was reported to have GA signal repression function as in the case of DELLA protein[25].TraesCS2D01G055700 could be a candidate gene for Rht8 and invites follow-up validation studies.Tian et al.[19]reported a height-reducing gene in wheat(Rht24),which was called a modern Green Revolution gene with 6.0-7.9 cm height-reduction effect,but with significant increase in kernel weight (2.0-3.4 g). This gene is among GA-sensitive height-reducing genes in wheat. In our study, we detected one MTA locus on 6A (QPLH-6A) at the position of 93,311,903 bp.However, one of the two markers(Xwmc256) linked to Rht24 as reported by Tian et al.[19]was located between 549,524,058 and 549,524,079 bp.This finding suggests that the MTA we identified on 6A is different from Rht24.

Several studies [1,5,7,26-28] have associated negative attributes to the widely used Rht genes(Rht-B1 and Rht-D1),including short coleoptile growth, reduced kernel weight, and decreased resistance to diseases. Alternative height-reducing genes may offer breeding solutions to reverse these negative effects. One such locus could be QPLH-2D.Approximately 100 kb away from QPLH-2D resides the gene TraesCS2D01G055700,which encodes a GRAS family transcription factor with 100%sequence identity with the Aegilops tauschii subsp. tauschii scarecrow-like protein 3(SCL3).An ortholog of this gene in Arabidopsis thaliana has been shown to be a positive regulator of GA signaling. The loss-offunction mutant scl3 showed a reduced-shoot-elongation phenotype [29,30]. Whether SCL3 is a candidate gene for Rht8 awaits further investigation.

Table 4-Previously reported reduced-height loci in wheat.

Fig.5-Approximate position of Rht8 determined on chromosome 2D according to the flanking markers reported by Gasperini et al.[23]and position of QPht.cau-2D.2 reported by Bellucci et al.[22].The positions of the markers were determined by Basic Local Alignment Search Tool(BLAST)analysis using the primer sequences of the markers in International Wheat Genome Sequencing Consortium(IWGSC,https://urgi.versailles.inra.fr/blast_iwgsc/blast.php).

In conclusion, our study suggests alternative heightreducing loci for further validation and possible breeding application.We also detected two well-known and widely used height-reducing genes (Rht-B1 and Rht-D1). We showed that combining association genetics with a position-dependent gene search strategy can suggest candidate genes for agronomically important MTAs for further study.

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2019.09.005.

Acknowledgments

The authors thank the National Small Grains Collection,USDA for providing seed of historical wheat accessions.Financial support from USDA Hatch Grant 1013073 to MM via Purdue College of Agriculture is greatly appreciated.