Bin-based genome-wide association studies reveal superior alleles for improvement of appearance quality using a 4-way MAGIC population in rice.

INTRODUCTION: The multiparental population provides us the chance to identify superior alleles controlling a trait for genetic improvement. Genome wide association studies at bin level (bin-GWAS) are expected to be more power in QTL mapping than GWAS at SNP level (SNP-GWAS).
OBJECTIVES: This study is to estimate genetic effects of QTL conferring grain appearance quality in rice by SNP-GWAS and bin-GWAS, compare their power in QTL mapping and identify the superior alleles of all detected QTL from 4 parents for genetic improvement.
METHODS: A 4-way MAGIC population and its four founders were cultivated in two environments to dissect the genetic basis of rice grain appearance quality. Both SNP-GWAS and bin-GWAS were conducted for QTL mapping. Multiple comparison among 4 parental bin/alleles was used to identify the superior alleles.
RESULTS: A total of 16 and 20 QTL associated with grain appearance quality were identified by SNP- and bin-GWAS, respectively. A minor chalkiness QTL qPGWC8.2/qDEC8 was assigned to a 30-kb genomic region, in which OsMH_08T0121900 is the potential candidate gene because its encoded protein, glucan endo-1,3-beta-glucosidase precursor is involved in the starch and sucrose metabolism pathway. The superior parental alleles for GS3, GL3.1, GW5, GW7, and Chalk5 and two QTLs were almost carried by the high-quality parents Cypress and Yuejingsimiao (YJSM), while the poor-quality parent Guichao-2 (GC2) always carried the inferior alleles. The top five recombinant inbred lines with the highest quality of grain shape and chalkiness traits all carried gene combinations of superior alleles.
CONCLUSIONS: Both SNP- and bin-GWAS methods are encouraged for joint QTL mapping with MAGIC population. qPGWC8.2/qDEC8 is a novel candidate gene strongly associated with chalkiness. The superior alleles of GS3, GW5, GL3.1, GW7, Chalk5 and qPGWC8.2 were identified, and the pyramiding of these superior alleles is helpful to improve rice appearance quality.

Introduction

Rice (Oryza sativa L.) is a major staple crop across Asian countries, feeding more than half of the world's population. Grain appearance has a considerable effect on the rice market demand, even though grain quality preferences vary among consumers worldwide. Uniform shape and a translucent endosperm are the best indicators for rice quality. Thus, grain shape and chalkiness are focal traits in rice genetics and breeding programs [1].

Conceptually, chalkiness, measured as the degree of endosperm chalkiness (DEC) and percentage of grains with chalkiness (PGWC), is a quantitative trait controlled by multiple genes. However, numerous quantitative trait loci (QTLs) have been detected using different types of populations according to the Gramene QTL database (). qPGWC7 has been fine mapped to a 44-kb region on chromosome 7 using a set of CSSLs, which was derived from a cross between PA64s and 9311 [2]. A total of 22 QTLs for appearance quality traits were identified with a population involving 66 chromosomal segment substitution lines (CSSLs) from the cross between Asominori (Jap.) and IR24 (Ind.) across eight environments, and nine QTLs were consistently detected [3]. Then, qPGWC8, a major QTL, was located at a 142-kb genomic region of chromosome 8 [4]. Three QTLs related to PGWC were detected on chromosomes 5, 8, and 10. They explained 50% of the phenotypic variation on F8 RILs derived from the cross between japonica rice varieties Koshihikari and C602 [5]. By utilizing high-throughput single nucleotide polymorphism (SNP) markers, ten common qPGWC and qDEC were identified in a CSSLs population with the japonica cultivar Nipponbare as donor parent and the indica cultivar ZS97 as recurrent parent [6]. A total of 19 QTLs associated with chalkiness on chromosomes 1, 4, 6, 7, 9 and 12, were identified in RILs of PA64s (Jap.) and 9311 (ind.) varieties, and qACE9 was fine mapped to a 22-kb genomic region on chromosome 9 [7]. Recently, qPCG1 was fine mapped to a 139-kb region on chromosome 1 using a residual heterozygous line [8]. Currently, only one major QTL (Chalk5) encoding a vacuolar H++-translocating pyrophosphatase has been cloned as a positive regulator on chromosome 5, and upregulation of Chalk5 expression increased chalkiness in the endosperm [9]. In addition to Chalk5, many mutants and transgenic lines have been used to identify genes involved in regulatory pathways of grain chalkiness formation including starch biosynthetic and metabolism pathways like SBE3, OsZIP58, GIF2 and OsBT1 [10], [11], [12], [13], seed storage protein biosynthesis like FLO2, GPA3, GLUP6 [14], [15], [16], and other chalkiness genes like OsAlaAT1, OsNFYB1 and Amy1A [17], [18], [19].

Grain length (GL), grain width (GW), and grain thickness (GT) reflected grain size or shape. A large number of QTLs control grain size. A total of 14 QTLs with significant effects on grain size have been cloned in rice [20]. Two of them (GW7 and GW6a) have a profound impact on both GL and GW; GW7 has opposite effects on GL and GW without reducing grain weight [21], [22], while GW6a had positive impacts on GL and GW and exhibits an extremely strong effect on grain weight [23], and the newly cloned gene WG7 increased the grain width [24]. The other 12 QTLs affect either grain length or grain width and resulted in varied grain weight; GW2, GS5, qSW5/GW5, and GW8 mainly controlled grain width [25], and GS2/GL2, OsLG3, qLGY3/OsLG3b, GS3, GL3.1/qGL3, GL4, TGW6, and GLW7 regulated grain length [25], [26], [27], [28], [29]. In addition to these major QTLs, hundreds of minor QTLs have been identified in the rice genome in the Gramene QTL database (). Likewise, SG3, a minor QTL, encoded an R2R3-MYB protein and negatively regulated GL [30].

Biparental population such as F2 and RILs populations usually restricts the number of detected QTLs and mapping resolution [31]. Additionally, linkage analysis and genome-wide association studies (GWAS) are efficient tools for the association mapping of phenotypic and genotypic data. GWAS takes advantage of natural variations in a germplasm collection to establish the association between traits and SNPs. Such studies utilize historically accumulated recombination to fine-map QTLs with high resolutions [7], [32], [33], [34], [35], [36], [37]. Since the cost of resequencing sharply decreased, GWASs have become an essential tool to dissect various agronomic traits because of their powerful mapping approach [38]. Although several advanced statistical models have been developed to control population structure effects, false positives are difficult to exclude because the noise of the population structure cannot be completely shielded yet [39].

Over the past decade, the MAGIC populations, a second-generation mapping resource, have raised novel opportunities in crops due to their complex pedigree structure. These populations offer great potential for both dissecting genomic structure and improving breeding populations [40], [41]. MAGIC populations harbor more allelic diversity than typical bi-parental mapping populations, and the multiple cycles of intercrossing provide more significant opportunities for recombination and, in turn, result in higher resolution for QTL location [42].

Recently, MAGIC populations have overcome the limitations of existing mapping populations, and have been recognized as some of the most potent mapping resources because of their diverse genetic multifounder contributions and weak population structure [43]. To date, scientists have established many MAGIC populations in cereals such as rice [44], [45], [46], [47], [48], barley [49], and maize [50].

In MAGIC populations, the primary input data for association mapping or linkage mapping are SNPs or parent probabilities. Each single-recombinant chromosome fragment (bin) in a progeny can be traced back to a specific parent, according to the pedigree. Thus, bin map development is a crucial step for executing multiple comparisons among parental alleles. A bin-mapping strategy can clarify the contributors of haplotypes and further display the breakpoints in the genome. MAGIC population is expected to have smaller haplotype blocks and higher mapping resolution than biparental population due to a higher recombination rate in MAGIC population. Thus, a MAGIC population can provide a chance to break the linkage drag between yield and yield-related traits [51].

In this study, a 4-way MAGIC population of 248 recombinant inbred lines (RILs) at the F6-7 generations and its parents were used to dissect the genetic basis of rice appearance quality. The main objective of the study was to identify loci associated with five-grain appearance quality traits with high resolution by combining SNP-GWAS and bin-GWAS. Beneficial parental alleles from these founders of the MAGIC population were identified for genetic improvement of rice appearance quality.

Materials and methods

Plant materials

A four-way MAGIC population of 248 F7 lines was constructed by crossing four founders (Table S1), including three indica (Xian) cultivars (IR34, Guichao-2 (GC2), Yuejingsimiao (YJSM)) and one japonica (Geng) cultivar (Cypress) [47], [52].

Field experiment, daily temperature recording and trait measurement

Twenty-five-days-old seedlings of the 4-way MAGIC population and its parents were transplanted into a one-row plot on 15th June 2015 and 2016. There were ten plants in each row with a planting density of 16.5*26.4 cm within and between rows. The field trial was carried out at Huazhong Agricultural University, Wuhan, China (114°21′ E, 30°28′ N), following a randomized complete block design with two replicates. Daily temperature was recorded in the field during the period from 1st June to 30th September every year using Ethernet temperature and humidity transmission recorder (Model Cos-03, Shandong Renke Measurement and Control Technology Co., Ltd., China).

The grain length (GL), width (GW), and length–width ratio (LWR) were captured using an image processing technology scanner (ScanMaker i800 Plus, Microtech Company, Hsinchu 30075, Taiwan) with the supplied software without enhancing the images and measured precisely using SmartGrain Software [53]. Three months after harvest, the physical and chemical properties of the seeds became stable. The tested samples were placed in a dry-ventilated area or an air-conditioned room for one week to keep the moisture content at approximately 14–16%. Fifty grams of seeds were dehulled to obtain 30 g of brown rice, which were then divided into three 10 g replicates and milled for 60 sec to reach the national standard for the first-class rice. Hence, after cooling, 10 g of seeds were sieved to a 1.5 mm diameter to remove the chaff and broken grains. Full rice grains were used for chalkiness evaluation. DEC, which was measured as the ratio of total chalky area to the total kernel area of all sampled grains, and PGWC were utilized to describe grain chalkiness [54]. PGWC and DEC were measured manually according to He and the National Standard of the People Republic of China (NSPRC 1999) [55].

Annotation of SNP variations and bin map construction

The sequencing depth for the four parents and each MAGIC RILs were 30× and 2×, respectively. The whole-genome variants were annotated with the general feature format (GFF) file of Minghui 63 as a reference genome () using SnpEff software to categorize the effects of variants on genome sequences [56]. Therefore, the bin map for the MAGIC population was constructed based on a sliding window approach. The exact method for bin development was described in our previous study, which determined the parental origin of each bin via identity by state and identity by descent between the MAGIC RILs and parents [47].

GWAS at the SNP level and bin level

A total of 843,505 high-quality SNPs were used for GWAS and to develop the bin map, and 5934 bins were developed for the genome. The SNP-based GWAS was conducted by FaST-LMM [57], and a kinship matrix was adopted to control for population structure, with a 3.7E-6 threshold obtained by improved Bonferroni methods [58]. The leading SNPs were obtained by the sliding window method, in which the window size was set as 1 Mb, and the step length was 1 SNP. The bin-based GWAS was conducted by the Random-B model in MagicQTL software [59], with a threshold of 3.0E-4 obtained by 1,000 permutations.

Genetic statistics

Two-way analysis of variance (ANOVA) was used to estimate genotype (G)-by-environment (E) interactions (G*E) for all investigated traits using ANOVA functions in R software (). ANOVA was used to calculate the phenotypic variation explained by each detected QTL with a linear model that included all peak markers [47]. The genetic effects of parental alleles at GS3, GL3.1, GW5, GW7, and Chalk5 were determined by multiple comparisons, which were performed with the Duncan test function in the agricolae package in R software [60]. GS3, GW5, and Chalk5 were amplified with designated primers for Sanger sequencing (Table S2). The PCR products were purified to carry out the sequencing reaction using the Big Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). Meanwhile, GL3.1 and GW7 were tested according to a 30x sequencing depth among the four parents of the MAGIC population.

Results

Variations in grain appearance quality traits

Parental performance in terms of grain appearance quality traits was presented in Table 1. GC2 had the highest values of DEC and PGWC and exhibited the smallest LWR. IR34 and Cypress had moderate chalkiness. On the contrary, YJSM was the highest-quality parent with the lowest values of DEC and PGWC in both years. All these parents exhibited significantly higher values for DEC and PGWC in 2016 than in 2015. In contrast, the MAGIC population exhibited wide continuous variations in grain appearance quality traits (Fig. 1). In general, the average of DEC and PGWC for both population and parents were more variable in 2016 than in 2015 (Table 1). However, the grain shape traits showed more differences in 2015 than in 2016, with equivalent average values in the two years.

Table 1

Phenotypic performance in terms of grain quality traits in the 4-way MAGIC parents and their heritability.

Traits	GC2	YJSM	Cypress	IR34	Parents (M ± SD)	Population (M ± SD)	Heritability %
DEC_2015%	11.5	1.5	0.8	4.3	4.5 ± 4.5	8.5 ± 9.0	83
DEC_2016%	34.4	4.4	34.1	24.5	24.3 ± 13.1	26.6 ± 15.6
PGWC2015 %	50.2	5.1	2.4	23.4	20.3 ± 20.4	26.4 ± 23.6	83
PGWC2016 %	92.1	14.5	69.4	70.5	61.6 ± 30.7	61.6 ± 25.9
GL2015 (mm)	4.94	6.09	6.19	6.26	5.87 ± 0.58	5.52 ± 0.72	89
GL2016 (mm)	4.74	5.65	5.75	6.19	5.58 ± 0.56	5.53 ± 0.49
GW2015 (mm)	2.57	1.97	1.92	2.14	2.15 ± 0.27	2.19 ± 0.19	96
GW2016 (mm)	2.39	1.85	1.91	1.95	2.03 ± 0.23	2.11 ± 0.18
LWR2015	1.92	3.10	3.24	2.93	2.80 ± 0.55	2.55 ± 0.41	94
LWR2016	1.98	3.07	3.03	3.19	2.82 ± 0.52	2.65 ± 0.35

DEC: degree of endosperm chalkiness, PGWC: percentage of grains with chalkiness, GL: grain length, GW: grain width, and LWR: length–width ratio. M: Mean, SD: standard deviation.

Fig. 1

Distribution of different appearance quality traits in two years in the 4-way MAGIC population.

The broad-sense heritability was 96% for GW, 94% for LWR, and 89% for GL, while chalkiness traits DEC and PGWC showed a moderate heritability of 83%. The traits were highly significantly affected by the environment except for GL. However, G*E was not significant for any of the studied traits (Table S3). DEC and PGWC were more significantly affected by the environment than were the grain shape traits.

Correlation among grain appearance quality traits

The phenotypic correlations between traits in 2015 were similar to those in 2016. The phenotypic correlations in a given character between the two years were highly significant (Fig. 2). Highly significant positive associations were observed between the chalkiness traits (DEC and PGWC) within a year and between years. Although a negative correlation was detected between GW and GL in the two years, it was not significant. As expected, GL and GW were positively and negatively correlated with LWR, respectively. Interestingly, DEC and PGWC in 2015 had significant positive correlation with GW, but they had no such significant correlation in 2016.

Fig. 2

Correlation heatmap between grain quality traits in the 4-way MAGIC population in 2015 and 2016. DEC: degree of endosperm chalkiness, PGWC: percentage of grains with chalkiness, GL: grain length, GW: grain width, and LWR: length–width ratio. The numbers 15 and 16 indicate the data from 2015 and 2016, respectively. The negatively correlated cells represented in blue, while the positively correlated cells represented in red. The darker color, the higher correlation was (P < 0.05; P < 0.01).

QTL detection by SNP-based GWAS

The genome-wide average LD decay in the MAGIC population was 2.5 Mb (Han et al. 2020). Nine and 12 QTLs were detected for the studied traits by the SNP-GWAS in 2015 and 2016 (Fig. 3, Fig. 4), respectively, and five QTLs were commonly detected (Table 2). A major QTL was mapped for DEC, explaining 13.1% of the phenotypic variance in 2015. Two QTLs for PGWC were distinguished each year, but no common QTLs were identified in either year. No QTLs were located near the known chalkiness genes.

Fig. 3

Manhattan (A, C, E, G) and Q-Q (B, D, F, H) plots showing associated SNP markers for the degree of chalkiness (DEC) and percentage of grains with chalkiness (PGWC) detected by a GWAS at SNP level in 2015 and 2016. The X-axis shows the chromosome number, and the Y-axis shows −log10 (p), while the horizontal line indicates the threshold p-value at significant level (p < 0.0001). Chalk5 is indicated in the plots.

Fig. 4

Manhattan (A, C, E, G, I, K) and Q-Q (B, D, F, H, J, L) plots showing associated SNP markers for grain length (GL), width (GW) and length–width ratio (LWR) detected by a GWAS at the SNP level in 2015 and 2016. The X-axis shows chromosome number and the Y-axis shows −log10 (p), while the horizontal line indicates the threshold p-value at significance level (p < 0.0001). GS3, GL3.1, GW5 and GW7 indicated in the plots.

Table 2

QTLs for appearance quality traits identified by the SNP-GWAS in both years.

Traits	Year	QTL name	Peak SNPs	Neighbor genes	p-value	PVE %
DEC	2015	qDEC3	3:6444105		9.3E-07	13.1
PGWC	2015	qPGWC2	2:32,812,602		2.8E-06	10.9
		qPGWC5.2	5:5,242,234	+2.10 Mb to Chalk5	6.8E-09	14.7
	2016	qPGWC8.1	8:3,697,671		1.6E-06	11.6
		qPGWC8.2	8:7,598,441		5.2E-07	4.7
GL	2015	qGL3.1	3:16,600,755	+133.2 kb to GS3	7.3E-08	19.3
		qGL3.4	3:26,307,574	+64.9 kb to GL3.1	2.5E-06	3.5
	2016	qGL3.1	3:16,489,110	+21.6 kb to GS3	2.8E-16	32.8
		qGL3.2	3:21,008,148		1.6E-09	1.2
		qGL3.3	3:23,451,981		3.0E-06	0.4
		qGL7	7:24,098,173	+368 kb to GW7	8.6E-09	7.3
GW	2015	qGW5	5:5,246,588	+16.0 kb to GW5	2.1E-18	39.7
		qGW7	7:23,684,827	−45.9 kb to GW7	1.4E-06	6.6
	2016	qGW5	5:5,228,933	+1.7 kb to GW5	7.6E-19	39.0
		qGW7	7:23,666,067	−64.7 kb to GW7	1.1E-06	4.6
LWR	2015	qLWR3.1	3:16,627,196	+161.3 kb to GS3	7.6E-08	18.4
		qLWR7	7:23,684,827		7.6E-09	9.1
	2016	qLWR3.1	3:16,569,106	+103.2 kb to GS3	6.7E-13	28.5
		qLWR3.2	3:21,599,202		1.1E-06	0.4
		qLWR5	5:5,228,933	+1.7 kb to GW5	7.0E-09	18.4
		qLWR7	7:23,684,990		9.5E-11	9.1

Peak SNP, the SNP showing a peak association signal, presented in the format of chromosome number followed by its physical position with the Minghui 63 genome as the reference. Neighbor genes represent the distance between peak SNPs and previously cloned genes. p-value, probability value; PVE, the percentage of phenotypic variance explained by QTL.

Five QTLs were detected for grain length. Among these QTLs, one major QTL, qGL3.1, was commonly detected in both years; four QTLs on chromosome 3 and one on chromosome 7 were detected in only one year. In 2015, the peak SNP in qGL3.1 placed 133.2 kb from GS3, and the peak SNP in the minor QTL qGL3.4 was 64.9 kb from GL3.1. In 2016, qGL3.1 mapped closer (21.6 kb) to GS3, with a larger contribution (32.8%) to grain length variation. The peak SNP in qGL7 was 368 kb from GW7. The grain width QTLs qGW5 and qGW7 were detected in both years. qGW5 was closely linked to GW5 with a distance of 16.0 kb and 1.7 kb, and explaining 39.7% and 39.0% of the phenotypic variance in 2015 and 2016, respectively. Similarly, qGW7 was located at positions 45.9 and 64.7 kb relative to GW7 in the two years. Five LWR QTLs were mapped to the designated locations of GW and GL QTLs. Of them, qLWR3.1 and qLWR7 were commonly detected in both years. Two major QTLs qLWR3.1 and qLWR5 were closely linked to GS3 and GW5, respectively.

QTL detection by bin-based GWAS

The 4-way MAGIC population genome was divided by 5934 bins (Table S4), with a bin length ranging from 10 to 970 kb and an average bin length of 56.3 kb. Most bins (85%) were shorter than 100 kb, while 47% of the bins were shorter than 50 kb [47]. The bin-based GWAS identified 20 QTLs for these investigated traits (Table 3). For Chalkiness traits, a total of eight QTLs were detected on chromosomes 3, 5, 6, and 8, but no QTLs were detected in both years. Accordingly, two main DEC QTLs, qDEC3.1, and qDEC6, were discovered. They were traced to B03_064, and B06_096 and explained 11.5% and 10.6% of the phenotypic variance, respectively. Similarly, two PGWC QTLs were mapped to Bin05_079 (qPGWC5) in 2015 and to Bin08_083 (qPGWC8.1) in 2016, which explained 13.4% and 12.4% of phenotypic variance, respectively. Interestingly, the minor QTLs qDEC8 and qPGWC8.2 were located in the same bin (Bin8_134) in 2016, which covered a 30-kb genome region.

Table 3

QTLs identified for appearance quality traits using high-density SNP bin-map.

Traits	QTLs	Peak bins	Interval (Mb)	Bin length	Neighbor genes	2015		2016
						p-value	PVE %	p-value	PVE %
DEC	qDEC3.1	B03_064	3.24–3.33	90 kb		2.5E-04	11.5
	qDEC3.2	B03_113	5.54–5.71	170 kb		2.1E-04	4.5
	qDEC3.3	B03_360	17.19–17.20	10 kb		2.6E-04	6.1
	qDEC6	B06_096	5.71–5.78	70 kb				3.7E-04	10.6
	qDEC8	B08_134	7.57–7.60	30 kb				3.6E-05	9.7
PGWC	qPGWC5.1	B05_079	4.59–4.64	50 kb	1.465 Mb to Chalk5	1.6E-03	13.4
	qPGWC8.1	B08_083	3.57–3.62	50 kb				6.8E-06	12.4
	qPGWC8.2	B08_134	7.57–7.60	30 kb				2.3E-07	8.4
GL	qGL3.1	B03_345	16.58–16.67	90 kb	+112.4 kb to GS3	4.9E-04	12.3	3.0E-12	35
	qGL3.2	B03_389	19.59–19.60	10 kb				2.1E-05	3.9
	qGL3.3	B03_482	26.04–26.08	40 kb	−202.7 kb to GL3.1	1.0E-03	4.0
	qGL7	B07_367	23.99–24.15	160 kb	+340 kb to GW7			2.1E-07	5.9
GW	qGW1.1	B01_024	1.12–1.15	30 kb		4.9E-04	8.3
	qGW1.2	B01_251	21.56–21.85	290 kb		5.2E-04	4.8
	qGW5	B05_094	5.19–5.23	40 kb	+0.6 kb to GW5	2.5E-18	30.4	7.1E-25	37.2
	qGW7	B07_361	23.36–23.53	170 kb	−285 kb to GW7	3.9E-06	7.3	1.2E-05	5.1
LWR	qLWR3.1	B03_347	16.73–16.76	30 kb	+262.5 kb to GS3			9.3E-08	25.8
	qLWR3.2	B03_389	19.59–19.60	10 kb				7.4E-05	2.8
	qLWR5	B05_093	5.18–5.19	10 kb	−40.6 kb to GW5			4.4E-07	12.7
	qLWR7	B07_361	23.36–23.53	170 kb	−285 kb to GW7	8.9E-06	5.9	5.7E-12	7.9

Peak bin, the bin showing a peak association signal. The physical positions of known genes were recorded according to the Minghui 63 reference genome. Neighbor genes represent the distance between peak bins and previously cloned genes. p-value, probability value; PVE, the percentage of phenotypic variance explained by the QTL.

For grain length, one major QTL, qGL3.1 (Bin03_345), was located 112.4 kb from GS3 and explained 12.3% and 35.0% of the phenotypic variance in 2015 and 2016, respectively. One minor QTL (qGL3.2) was detected in the bin B03_482 in 2015, which was 203 kb from GL3.1. Additionally, two minor QTLs (qGL3.2 and qGL7) were discovered in 2016. qGL7 was mapped to the bin 340 kb from GW7. For grain width, four QTLs were detected in the two years. Among these QTLs, two (qGW5 and qGW7) were strongly expressed in both years. The major QTL qGW5 located in B05_094 was adjacent to the bin harboring GW5, and explained 30.4% and 37.2% of the phenotypic variance in 2015 and 2016, respectively. qGW7, a minor QTL, was located 285 kb to GW7. Moreover, two minor QTLs (qGW1.1 and qGW1.2) were mapped on chromosome 1. For LWR QTLs, qLWR3.1 was assigned to B03_347, which was 262.5 kb from GS3; qLWR5 was mapped to B05_093 with a distance of 40.6 kb from GW5; and the minor QTL qLWR7 located in B07_361 was detected in both years.

Comparison of genetic effects among parental alleles at four-grain shape genes

Four QTLs identified in this population were closely linked to GS3, GL3.1, GW5, and GW7 (Table 3, Table 4). Hence, we compared the parental genetic effects of these genes on the evaluated grain shape traits.

Table 4

Parental allele effects of four QTLs on grain shape in two environments.

QTLs/gene/Bin	Hap	Freq	GL15 (mm)	GL16 (mm)	GW15 (mm)	GW16 (mm)	LWR15	LWR16
qGL3.1/GS3/B03_343	GC2	89	5.2 ± 0.6b	5.2 ± 0.4b	2.2 ± 0.2 a	2.1 ± 0.2 a	2.4 ± 0.4b	2.4 ± 0.3b
	YJSM	41	5.7 ± 0.6 a	5.6 ± 0.5 a	2.1 ± 0.2 a	2.0 ± 0.2 a	2.7 ± 0.4 a	2.8 ± 0.4 a
	Cypress	23	5.7 ± 0.6 a	5.8 ± 0.3 a	2.2 ± 0.2 a	2.1 ± 0.2 a	2.6 ± 0.3 a	2.7 ± 0.3 a
	IR34	81	5.8 ± 0.7 a	5.8 ± 0.4 a	2.2 ± 0.2 a	2.1 ± 0.2 a	2.7 ± 0.4 a	2.8 ± 0.3 a
qGL3.3/GL3.1/B03_485	GC2	71	5.1 ± 0.7b	5.3 ± 0.5b	2.2 ± 0.2 a	2.1 ± 0.2 a	2.4 ± 0.4b	2.5 ± 0.3b
	YJSM	55	5.7 ± 0.7 a	5.6 ± 0.4 a	2.1 ± 0.2b	2.0 ± 0.2b	2.7 ± 0.4 a	2.8 ± 0.3 a
	Cypress	44	5.6 ± 0.5 a	5.5 ± 0.6 a	2.2 ± 0.2 a	2.2 ± 0.2 a	2.5 ± 0.3 a	2.6 ± 0.3b
	IR34	75	5.7 ± 0.7 a	5.7 ± 0.4 a	2.2 ± 0.2 a	2.1 ± 0.2 a	2.6 ± 0.4 a	2.8 ± 0.3 a
qGW5/GW5/B05_095	GC2	70	5.5 ± 0.7 a	5.5 ± 0.5 a	2.3 ± 0.2 a	2.3 ± 0.2 a	2.4 ± 0.3b	2.5 ± 0.3b
	YJSM	73	5.5 ± 0.8 a	5.5 ± 0.6 a	2.1 ± 0.2b	2.0 ± 0.1b	2.6 ± 0.4 a	2.7 ± 0.4 a
	Cypress	49	5.6 ± 0.8 a	5.6 ± 0.4 a	2.1 ± 0.1b	2.1 ± 0.1b	2.6 ± 0.4 a	2.7 ± 0.3 a
	IR34	55	5.6 ± 0.6 a	5.5 ± 0.5 a	2.1 ± 0.2b	2.0 ± 0.2b	2.7 ± 0.4 a	2.7 ± 0.4 a
qGW7/GW7/B07_365	GC2	57	5.4 ± 0.7b	5.4 ± 0.4b	2.2 ± 0.2 a	2.1 ± 0.2 a	2.5 ± 0.4b	2.6 ± 0.3b
	YJSM	74	5.6 ± 0.6 ab	5.4 ± 0.4b	2.2 ± 0.2 a	2.1 ± 0.2 a	2.5 ± 0.3b	2.6 ± 0.3b
	Cypress	43	5.8 ± 0.7 a	5.8 ± 0.5 a	2.1 ± 0.2b	2.0 ± 0.2b	2.8 ± 0.4 a	2.9 ± 0.4 a
	IR34	56	5.4 ± 0.8b	5.5 ± 0.6b	2.2 ± 0.2 a	2.1 ± 0.2 a	2.5 ± 0.5b	2.6 ± 0.4b

QTLs/gene/bin, QTLs identified in this study/linked known genes/the bin known gene located in.

In the MAGIC population, 89, 41, 23, and 81 RILs carried parental bins (B03_345 harboring GS3) of GC2, YJSM, Cypress, and IR34, respectively. GS3 had pleiotropic effects on GL and LWR (Table 4). Comparison sequencing of GS3 (Fig. S1a, b) showed that YJSM, IR34 and Cypress carried frameshift mutation (C165A-Hap1) in GS3-exon2 that resulted in a premature stop codon [61]. This haplotype (GS3-Hap1) was associated with longer grains, while GS3-Hap2 presented in GC2 was linked to a shorter grain length (Table S5). The phenotypic performance of GS3 alleles was visualized as 0.4–0.6 mm differences in grain length between Hap1 and Hap2.

GL3.1 increased not only GL but also GW in the MAGIC population. Seventy-one RILs carrying GL3.1GC2 had the shortest GLs of 5.12 and 5.25 mm, which were 0.4–0.6 mm shorter than those of lines carrying other parental GL3.1 allele (5.5–5.73 mm). Moreover, the YJSM allele decreased GW (Table 4). The variation in sequencing data showed that IR34 and YJSM shared the same GL3.1 allele (Hap2) and Cypress carried Hap3, which contained a typical 1-bp insertion at 3′-UTR (Table S6), while GC2 carried Hap1 with a 1-bp deletion.

For GW5, the lines with IR34, YJSM and Cypress alleles exhibited no considerable differences in GW, but they had significantly narrower grains (0.21–0.26 mm) than the 70 RILs carrying GW5GC2 allele (Table 4). Moreover, these alleles except GW5GC2 increased LWR. The variation in sequencing data showed that the GW5GC2 allele had one different amino acid (Pro63Ser) (Fig. S2a, b) and a 4-bp insertion in the site of 1268 bp downstream, which was not identified in the other three parents' alleles (Table S7). This insertion of the GW5GC2 allele likely resulted in a functional change in GW5.

The 36 MAGIC lines carrying the GW7 allele not only had a smaller average GW but also had a longer GL by 0.3–0.4 mm (and thus a larger LWR) than the lines carrying other parental GW7 alleles. However, no significant differences among the three parental GW7 alleles (except the Cypress allele) were observed across the years (Table 4). Comparative sequencing analysis showed that GC2, YJSM, and IR34 carried the same GW7 allele (Hap1), and Cypress carried a different allele of Hap2 (Table S8).

Comparison of genetic effects on DEC and PGWC among parental alleles

Although qPGWC5.2 was only detected by the SNP-GWAS and qPGWC5.1 was only identified by the bin-GWAS, they were located 2.11 Mb and 1.46 Mb from Chalk5, which was located in B05_048. Therefore, we suggest that Chalk5 underlies qPGWC5 because only a major gene for chalkiness has been frequently reported on chromosome 5 [62]. Then the MAGIC RILs were divided into four parental types according to the origins of B05_048. Seventy-three RILs carrying Chalk5GC2 exhibited a higher PGWC than those carrying the other three parental alleles (Table 5).

Table 5

Allele effects of Chalk5, qPGWC8.1 and qPGWC8.2 on chalkiness traits in two environments.

QTL/gene/bin	Hap	Freq.	DEC 2015%	DEC 2016%	PGWC 2015%	PGWC 2016%
qPGWC5/Chalk5/b05_048	GC2	73	11.0 ± 11.2 a	28.4 ± 15.8 a	35.1 ± 28.8 a	66.6 ± 24.7 a
	YJSM	48	9.0 ± 9.2 ab	28.1 ± 12.3 a	24.5 ± 20.4b	54.9 ± 22.0b
	Cypress	48	6.4 ± 6.2b	28.5 ± 16.7 a	20.6 ± 17.5b	64.6 ± 26.2 ab
	IR34	74	6.9 ± 7.2b	22.4 ± 16.2b	23.2 ± 22.2b	52.1 ± 27.7b
qPGWC8.1/b08_083	GC2	88	10.0 ± 9.6 a	29.8 ± 14.2 a	30.9 ± 26.0 a	68.5 ± 22.7 a
	YJSM	41	8.7 ± 8.4 ab	29.9 ± 16.3 a	27.8 ± 24.1 ab	67.1 ± 24.0 a
	Cypress	43	5.9 ± 5.6b	17.1 ± 16.2b	20.3 ± 17.4b	42.3 ± 29.6b
	IR34	64	7.5 ± 9.5b	25.7 ± 14.2 a	22.1 ± 20.8b	60.7 ± 22.6 a
qPGWC8.2/b08_134	GC2	74	10.6 ± 10.8 a	31.6 ± 14.6 a	31.9 ± 27.5 a	71.6 ± 20.6 a
	YJSM	45	8.1 ± 7.3 ab	25.0 ± 14.8b	25.7 ± 22.2 ab	59.4 ± 25.6b
	Cypress	45	5.0 ± 4.2b	16.2 ± 14.9c	17.8 ± 14.2b	41.7 ± 28.7c
	IR34	71	8.6 ± 9.6 a	28.9 ± 14.4 ab	26.3 ± 23.8 ab	65.7 ± 21.7 ab

QTLs/gene/bin, QTLs identified in this study/linked known genes/the bin known gene located in.

Assuming that alleles decreasing the DEC and PGWC are superior alleles, Chalk5IR34, qPGWC8.1Cypress and qPGWC8.2Cypress are superior alleles, while Chalk5GC2, qPGWC8.1GC2 and qPGWC8.2GC2 are inferior alleles.

Meanwhile, 11 SNPs and one 9-bp InDel were identified in the coding sequence of Chalk5 among the four parents, and all the SNPs caused amino acid changes (Table S9). According to the 12 polymorphic sites, three haplotypes were identified in the four parents. GC2, which carried Chalk5-Hap1, had two different amino acids (Ser281Ile, Arg282Gly) compared with those in the other parents, which carried Serine and Arginine in the same location (Fig. S3a, b); YJSM and IR34 carried the same Chalk5-Hap2. In addition, there are ten polymorphic sites between the Cypress haplotype (Hap3) and Hap2.

Both the SNP-GWAS and bin-GWAS detected qPGWC8.1 and qPGWC8.2 only in the hot season of 2016. The comparative analysis of parental alleles for both QTLs revealed that the RILs carrying Cypress allele had significantly lower chalkiness percentage than those lines carrying other parental alleles. Moreover, for qPGWC8.2, the Cypress allele had about 18% lower PGWC and 9% lower DEC than the allele carried by the highest-quality parent, YJSM, in 2016.

Validation of QTLs in the lines with extreme phenotypes of chalkiness and grain shape

Supposing the alleles that increased the LWR are superior alleles and the ones that decreased the LWR are inferior alleles, the lines with long and slender grains would be expected to contain combinations of superior alleles and vice versa. To validate this expectation, we selected the five lines with the best grain shape (long and slender) and five lines with short, round grains and identified their GS3, GL3.1, GW5 and GW7 gene combinations (Table S10). Indeed, the five best lines with an LWR greater than 3.2 carried the superior alleles GS3IR34, GW5YJSM, GL3.1IR34, and GW7Cypress (Table 4). The alleles together increased GL, decreased GW, and in turn led to long and slender grains. In contrast, the five lines with an LWR of less than 2.3 carried the inferior alleles GS3, GW5, GL3.1 and GW7GC2 or GW7. The alleles together decreased GL, increased GW, and in turn led to short and round grains.

Accordingly, we compared the five best lines and the five worst lines, according to chalkiness (Table S11). The best lines had very limited chalkiness (DEC less than 6% and 8%, and PGWC less than 18% and 20% in 2015 and 2016, respectively). The best lines carried the superior alleles qPGWC8.1Cypress, qPGWC8.2Cypress, and Chalk5IR34 or Chalk5YJSM, whereas the worst lines displayed a significant degree of chalkiness (DEC greater than 14% and 18%, and PGWC greater than 65% and 86% in 2015 and 2016, respectively). All the worst lines carried the inferior alleles from parent GC2 at the three targeted QTLs.

Discussion

Superior alleles for the improvement of grain appearance quality identified by multiple comparisons at bin-level

The identification of superior alleles in nature is beneficial to genetic development in rice. In this study, we identified the favorable alleles for four-grain shape genes by multiple comparisons (Table 4). The GS3YJSM, GS3Cypress, and GS3IR34 alleles are nonfunctional and contribute to long grains, while GS3GC2 is a functional allele leading to a short grain. Therefore, the nonfunctional GS3 alleles carried by the three parents other than GC2 resulted in long, slender grains and therefore were regarded as the superior alleles for grain shape. This result was consistent to previous studies, showing that the loss of function of GS3 donated by a 178-aa truncation in the C-terminus of the predicted protein led to slender grains [63]. Accordingly, the GL3.1 superior allele from YJSM, Cypress and IR34 had no effect on GW [64].

For GW5, the functional allele donated by GC2 increased GW but decreased GL and thus resulted in short, round grains. GW7 increased LWR by increasing GL but decreasing GW, which is in accordance with the findings of previous studies [21], [22], [65]. Therefore, GW7 is the superior allele, and all other parental alleles are inferior, leading to short, round grains. These superior alleles were validated by comparing the top five RILs with long, slender grains and five RILs with short, round grains, which carried 4-gene combinations of superior alleles and inferior alleles, respectively (Table S10). As expected, the parent Cypress with the highest grain shape quality harbored the superior alleles for all four genes, and the worst-performance parent GC2 carried the inferior alleles for all four genes. Interestingly, the high-quality parent YJSM carried an inferior GW7 allele, indicating the potential to improve its grain shape.

For chalkiness-related QTLs, the poor-quality parent GC2 contributed inferior alleles at three detected QTLs (Table 5). Cypress alleles seemed to be superior at qPGWC8.1 and qPGWC8.2, while IR34 donated the superior allele of Chalk5. Without exception, the top five RILs with the highest and lowest chalkiness content carried the gene combinations of superior alleles and inferior alleles, respectively (Table S11). All these superior alleles are beneficial for breeders to improve rice appearance quality by marker-aided selection.

High temperature has major effects on chalkiness and minor impacts on grain shape

The percentage of chalkiness in rice grain is an index that determines appearance quality. Chalky grains contained a lower density of starch granules than vitreous grains. Chalkiness-free rice without white core, white back, and white belly can significantly improve milling and cooking quality [66]. The traits for grain shape exhibited no significant differences between the two years, indicating that grain shape is of highly heritable and relatively stable across environments [67]. However, a significant difference in chalkiness was detected between the two years. The phenotypic values of parents and population were much greater (3-folds) in 2016 than those in 2015 (Table 1). The main reason was that the average temperature in the grain filling stage from 1st July to 20th August was 3.4℃ higher in 2016 than in 2015 (Fig. S4), which indicated that chalkiness-related traits were significantly affected by the environment. At the cellular level, moderate heat stress (34 °C) resulted in precocious endosperm cellularization. Hence, the seeds under high heat stress (42 °C) failed to cellularize [68]. High temperatures during the grain filling stage caused a lower density of starch granules, leading to a high PGWC and DEC [69], [70]. Tests of the japonica rice cultivar “Akitakomachi” under high- and low-temperature conditions in both open field and artificial greenhouse revealed results similar to those mentioned above [71]. Additionally, plants grown at high temperatures produced a larger number of immature grains and less translucent grains than those grown at low temperatures [71]. Consequently, both moderate and elevated temperatures negatively affected grain yield and grain quality in rice.

The inbred and hybrid japonica and indica rice varieties in northern China, with better climatic conditions for rice cultivation, showed better grain filling with a low PGWC and DEC than those in southern Chinese regions with high temperatures [72]. Consequently, to enhance rice quality in the production system, the timing of the grain filling stage should coincide with periods of warm temperatures. It is somewhat surprising that no common QTLs for chalkiness-related traits were identified between the years. Perhaps these QTLs detected in two years have different responses to heat stress. Moreover, high temperature led to a significant variation in chalkiness content in the parents (Table 1) and MAGIC population (Fig. 1 and Table 5), which was helpful for QTL mapping.

Bin-based QTLs detection has more power in QTL mapping

GWAS at the SNP level was conducted in several MAGIC populations [41], [44], [48], [49], [50], [73]. Nevertheless, the advantage of the MAGIC population in QTL mapping was not fully realized because genetic pedigree information was not utilized. The international rice research institute (IRRI) used the SNP-based GWAS for QTL detection with a 16-way MAGIC global population. Its results strongly indicates that the higher recombination rate in MAGIC population than the biparental population leads to a higher resolution of QTL detection in the MAGIC population [48]. A comparison of QTL mapping results for heading date showed that bin-GWASs are more powerful in QTL mapping than SNP-GWAS with the same 4-way MAGIC population used in this study [47]. It is noted that the four parents could be classified into only two types according to SNP genotypes. However, the four parents could be classified into as many as four types, allowing the MAGIC lines to be divided into four parental groups and multiple comparisons in any tested bins. Therefore, it is expected that bin-GWASs are more powerfully for identifying significant differences among parents. For example, there were three haplotypes among the four parents at GL3.1, GW5, and Chalk5. Also, the MAGIC lines were therefore divided into three groups for QTL analysis according to their parental haplotypes, which allowed for a fine-resolution comparison between any paired parental bins.

Here, 16 and 20 QTLs were identified for these five traits related to grain appearance quality in the SNP-GWAS and bin-GWAS, respectively (Table 2, Table 3). The major QTLs, such as those linked to GS3, GW5 and GW7, were commonly identified by both methods, but the bin-GWAS identified more minor QTLs than SNP-GWAS. Five QTLs for DEC were identified in the bin-GWAS, whereas only one QTL was detected in the SNP-GWAS. For example, the QTL qDEC8 detected only by the bin-GWAS was mapped to the same region as a previously identified QTL [1].

In the last decade, four chalkiness QTLs of qPGWC-7 [2], qACE9 [7], qPGWC8 [4] and qPGC1 [8] were fine-mapped in rice. In this investigation, the detected chalkiness QTLs were located in chromosomes 2, 3, 5, 6 and 8. Among these QTLs, qPGWC8 was placed into a 7.57–7.60 Mb region on chromosome 8, which was closely linked to the interval containing a previously mapped chalkiness QTL (qPGWC8; 1.5–7.5 Mb) [1]. Most likely, both QTLs were the same gene underlying chalkiness. Interestingly, this 30-kb genomic region includes four predicted genes. Among these predicted genes, OsMH_08T0121900 (http://rice.hzau.edu.cn/rice/) encoded glucan endo-1,3-beta-glucosidase, which is annotated as a member of the glycoside hydrolase family 17 (PF00332). Glycoside hydrolase is involved in carbohydrate metabolic process (http://pfam.xfam.org/family/pf00332). This data agrees with the result obtained from the transcriptomic analysis of 15 DAF caryopses of a high chalkiness NIL with its parental lines (CSSL50-1 and Asominori) for the assessment of grain endosperm chalkiness [74], which confirmed that glycosyl hydrolases family 1, 16, and 17 were involved in starch, sucrose and carbohydrate metabolic process and interacted with rice endosperm development [74]. This predictive gene is likely the potential candidate for qDEC8, which will be validated in future work. This information indicated the reliability of minor QTLs identified by bin-GWAS.

It is noted that no methods have precisely mapped Chalk5. One possible reason is that the chalkiness traits were environmentally dependent (Table 1, Table 6). In the MAGIC population, the lines even carrying the same parental Chalk5 allele had extremely different chalkiness phenotypes because they had drastically different heading dates [47]. To improve the precision of such kinds of QTL mapping, an advanced statistical model is still required to block the environmental effect. On the other hand, advanced experimental designs minimizing the environmental effect by synchronizing the heading date of the MAGIC population are more important [75]. There was no clear distinction in the resolution of QTL mapping between the two methods because some QTLs were mapped with a higher resolution by the SNP-GWAS, and others were mapped with a higher resolution by the bin-GWAS (Table 5, Table 6). Therefore, both methods are suggested for QTL mapping with MAGIC population.

Table 6

Allele effects of Chalk5, qPGWC8.1 and qPGWC8.2 on chalkiness traits in two environments

QTL/gene/bin	Hap	Freq.	DEC 2015 %	DEC 2016 %	PGWC 2015 %	PGWC 2016 %
qPGWC5/Chalk5/b05_048	GC2	73	11.0 ± 11.2 a	28.4 ± 15.8 a	35.1 ± 28.8 a	66.6 ± 24.7 a
	YJSM	48	9.0 ± 9.2 ab	28.1 ± 12.3 a	24.5 ± 20.4 b	54.9 ± 22.0 b
	Cypress	48	6.4 ± 6.2 b	28.5 ± 16.7 a	20.6 ± 17.5 b	64.6 ± 26.2 ab
	IR34	74	6.9 ± 7.2 b	22.4 ± 16.2 b	23.2 ± 22.2 b	52.1 ± 27.7 b
qPGWC8.1/b08_083	GC2	88	10.0 ± 9.6 a	29.8 ± 14.2 a	30.9 ± 26.0 a	68.5 ± 22.7 a
	YJSM	41	8.7 ± 8.4 ab	29.9 ± 16.3 a	27.8 ± 24.1 ab	67.1 ± 24.0 a
	Cypress	43	5.9 ± 5.6 b	17.1 ± 16.2 b	20.3 ± 17.4 b	42.3 ± 29.6 b
	IR34	64	7.5 ± 9.5 b	25.7 ± 14.2 a	22.1 ± 20.8 b	60.7 ± 22.6 a
qPGWC8.2/b08_134	GC2	74	10.6 ± 10.8 a	31.6 ± 14.6 a	31.9 ± 27.5 a	71.6 ± 20.6 a
	YJSM	45	8.1 ± 7.3 ab	25.0 ± 14.8 b	25.7 ±22.2 ab	59.4 ± 25.6 b
	Cypress	45	5.0 ± 4.2 b	16.2 ± 14.9 c	17.8 ±14.2 b	41.7 ± 28.7 c
	IR34	71	8.6 ± 9.6 a	28.9 ± 14.4 ab	26.3 ± 23.8 ab	65.7 ± 21.7 ab

QTLs/gene/bin, QTLs identified in this study/linked known genes/the bin known gene located in.

Assuming that alleles decreasing the DEC and PGWC are superior alleles, Chalk5IR34, qPGWC8.1Cypress and qPGWC8.2Cypress are superior alleles, while Chalk5GC2, qPGWC8.1GC2 and qPGWC8.2GC2 are inferior alleles.

Conclusion

This 4-way MAGIC population is one of the latest mapping resources for exploring the grain appearance quality through SNP- and bin-based GWAS methods. The combination of the two methods proved to be a more powerful tool for identifying significant candidate QTLs/genes. This investigation discovered five known genes and 15 novel QTLs for grain shape and grain chalkiness traits across the environments. qPGWC8.2/qDEC8, which was verified in 2016, carried a novel candidate gene strongly associated with chalkiness traits. This MAGIC population presented the superior and inferior alleles of GS3, GW5, GL3.1, GW7, Chalk5 and qPGWC8.2, which were carried by the MAGIC four founders of the MAGIC population. Additionally, the pyramiding of superior alleles of detected genes/QTLs can help rice breeders to develop new rice varieties with high grain quality. The validation of the novel candidate gene will be a goal in future work.

Ethical statement

This article does not contain any studies with human or animal subjects.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.